Firn model intercomparison experiment (FirnMICE)

Jessica M.D. Lundin, C. Max Stevens, Robert Arthern, Christo Buizert, Anais
Orsi, Stefan R.M. Ligtenberg, Sebastian Simonsen, Evan Cummings, Richard
Essery, Will Leahy, et al.

► To cite this version:

Jessica M.D. Lundin, C. Max Stevens, Robert Arthern, Christo Buizert, Anais Orsi, et al. Firn model intercomparison experiment (FirnMICE). Journal of Glaciology, 2017, 63 (239), pp.401-422. 10.1017/jog.2016.114. hal-03226722

HAL Id: hal-03226722

https://hal.science/hal-03226722

Submitted on 16 May 2021

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

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

<figure><img src="image_1.png" /></figure>
Firn Model Intercomparison Experiment (FirnMICE)

JESSICA M.D. LUNDIN,<fnref n="1" /> C. MAX STEVENS,<fnref n="1" /> ROBERT ARTHERN,<fnref n="2" /> CHRISTO BUIZERT,<fnref n="3" />

ANAIS ORSI,<fnref n="4" /> STEFAN R.M. LIGTENBERG,<fnref n="5" /> SEBASTIAN B. SIMONSEN,<fnref n="6" />

EVAN CUMMINGS,<fnref n="7" /> RICHARD ESSERY,<fnref n="8" /> WILL LEAHY,<fnref n="1" /> PAUL HARRIS,<fnref n="1" />

MICHIEL M. HELSEN,<fnref n="5" /> EDWIN D. WADDINGTON<fnref n="1" />

ABSTRACT. Evolution of cold dry snow and firn plays important roles in glaciology; however, the phys-
ical formulation of a densification law is still an active research topic. We forced eight firn-densification
models and one seasonal-snow model in six different experiments by imposing step changes in tempera-
ture and accumulation-rate boundary conditions; all of the boundary conditions were chosen to simulate
firn densification in cold, dry environments. While the intended application of the participating models
varies, they are describing the same physical system and should in principle yield the same solutions. The
firn models all produce plausible depth-density profiles, but the model outputs in both steady state and
transient modes differ for quantities that are of interest in ice core and altimetry research. These differ-
ences demonstrate that firn-densification models are incorrectly or incompletely representing physical
processes. We quantitatively characterize the differences among the results from the various models.
For example, we find depth-integrated porosity is unlikely to be inferred with confidence from a firn
model to better than 2 m in steady state at a specific site with known accumulation rate and temperature.
Firn Model Intercomparison Experiment can provide a benchmark of results for future models, provide a
basis to quantify model uncertainties and guide future directions of firn-densification modeling.

KEYWORDS: ice-sheet modelling, mass balance, polar firn

# 1. INTRODUCTION

Snow that falls in the cold, dry interiors of polar ice sheets slowly transforms into ice through an intermediate stage called firn. Firn density, *ρ*, increases with depth and time due largely to overburden stress from accumulation of new snow. The densification of firn is often divided into three zones. In zone 1 (*ρ* < 550 kg m⁻<fnref n="3" />), firn becomes denser due to grain-boundary sliding and grain growth (Alley, 1987). The density of 550 kg m⁻<fnref n="3" /> (porosity 40%) marks the transition to zone 2 and corresponds to the maximum packing density of uniform spheres (Cuffey and Paterson, 2010, p. 21). In zone 2 (550 kg m⁻<fnref n="3" /> < *ρ* < 830 kg m⁻<fnref n="3" />), densification occurs primarily through sintering. During sintering, grain deformation occurs at the interface between grains. In this process the distance between the midpoints of neighboring grains is reduced, thereby densifying the firn. The transition to zone 3 occurs at the bubble close-off (BCO) density of ~830 kg m⁻<fnref n="3" />. Bubbles of air are trapped and are isolated from the overlying atmosphere. In zone 3, densification occurs through bubble compression until the density of ice is reached at ~917 kg m⁻<fnref n="3" />. These processes mentioned are the dominant ones; however, other processes may be important. Detailed observations of firn density profiles show that the transitions between the different zones are not as abrupt and clear as implied here (Hörhold

and others, 2011), suggesting that the transition between
these dominant modes of densification occurs gradually. The
discussion here focuses on the dry snow zone, where summer-
time melt occurs only very rarely and meltwater percolation and
refreezing do not need to be considered.

Understanding the processes of firn evolution including densification is important for several applications in glaciology. For ice core interpretation, firn-densification models are needed to establish a consistent chronology for the ice and for the gases trapped within it (Schwander and others, 1997; Goujon and others, 2003). Bubbles of air trapped in ice provide a direct measurement of past atmospheric composition, and stable isotopes in the ice provide a proxy for temperature and/or moisture source. The age of gas trapped in bubbles differs from that of the surrounding ice; this gas-age/ice-age difference (Δage) arises because gases diffuse rapidly through the porous and permeable firn (Schwander and Stauffer, 1984); for example, modern air is being trapped today in ice layers that can be hundreds to thousands of years old. Firn-densification models are needed to find the age of the firn at the lock-in depth (i.e. depth at which air can no longer exchange with the free atmosphere) in order to calculate Δage (e.g. see Buizert and Severinghaus, 2016).

Satellite and airborne altimetry using radar or LiDAR allow
centimeter-scale-resolution measurements of surface elevations

<footnote n="1">¹Department of Earth and Space Sciences, University of Washington, Seattle, WA, USA</footnote>

<footnote n="2">²British Antarctic Survey, Natural Environment Research Council, Cambridge, UK</footnote>

<footnote n="3">³College of Earth, Ocean and Atmospheric Sciences, Oregon State University, Corvallis, OR, USA</footnote>

<footnote n="4">⁴Laboratoire des Sciences du Climat et de l'Environnement, France</footnote>

<footnote n="5">⁵Institute for Marine and Atmospheric Research Utrecht (IMAU), Utrecht, The Netherlands</footnote>

<footnote n="6">⁶DTU Space, Technical University of Denmark, Kgs. Lyngby, Denmark</footnote>

<footnote n="7">⁷Department of Computer Science, University of Montana, Missoula, MT, USA</footnote>

<footnote n="8">⁸School of Geosciences, University of Edinburgh, Edinburgh, UK</footnote>

<footnote>Correspondence: Edwin D. Waddington <edw@uw.edu></footnote>
of glaciers and ice sheets, which can be used to measure volume changes (e.g. Wingham and others, 1998; Zwally and others, 2005; Helm and others, 2014). A firn-densification model must be used to convert measured volume changes to mass changes. Firn models are also used to determine density-depth profiles in ice shelves and mountain glaciers to correctly model their rheological properties when firn comprises a large fraction of their thickness (Gagliardini and Meyssonnier, 1997; Lüthi and Funk, 2000; Zwinger and others, 2007).

Uncertainties in firn-densification physics introduce uncertainties into each of these firn-model applications. Uncertainty from firn modeling is the largest contributor of uncertainty in mass-balance estimates from altimetry methods, where the thickness change from the densification can be greater than the change in thickness due to mass-balance changes (Helsen and others, 2008; Shepherd and others, 2012). In ice core research, fundamental lead-lag analysis between temperature changes and changes in atmospheric-gas concentrations requires a precise understanding of the Δage. Often the uncertainty in Δage estimation is larger than the lag between temperature change and atmospheric CO₂ concentration change (Parrenin and others, 2012).

In steady state at a given site (constant climate, no change in the firn density-depth profile and no ice-sheet thickness change) the firn-densification rate at any depth is also constant through time; the mass of the snow added by accumulation each year equals the mass of ice removed from the bottom of the firn column by downward ice flow. In reality, firn columns in ice sheets and glaciers are seldom, if ever, in the steady state envisioned in Sorge's Law (e.g. Bader, 1954; Cuffey and Paterson, 2010, p. 19). For example, the difference between gas age and ice age is highly variable through space and time because firn densification depends strongly on accumulation rate and temperature, both of which vary spatially and temporally. Mass-balance studies of ice sheets can have significant uncertainty in estimating the mass of polar firn due to spatial and temporal variations in firn densification (Pritchard and others, 2012).

The uncertainties in firn-densification models arise from the fact that no accurate, physically-based, widely applicable and widely accepted constitutive law yet exists for snow densification. Firn modeling has largely been an empirical process, calibrating models to present-day observed density profiles (e.g. Herron and Langway, 1980; Spencer and others, 2001). Some models determine effective viscosity, or firn stiffness, in a constitutive formulation that relates known overburden stress to measured strain densification rates (e.g. Arthern and others, 2010; Morris and Wingham, 2014). These models take firn densification one step closer to a fully physics-based set of equations that can model transient behavior. The effective viscosity is generally calibrated empirically to firn texture such as grain size (e.g. Arthern and others, 2010), to thermal history (e.g. Morris and Wingham, 2014), or to chemical impurities that control texture (e.g. Freitag and others, 2013). Arnaud and others (2000) described a model based on physics of grain sliding and plastic deformation, but its structural parameters are still empirical.

Three questions arise about firn models:

(1) When models are based on similar physical concepts, do their results agree? (inter-model variance);

(2) When models are calibrated under steady-state conditions, do they produce similar histories of change under transient conditions? (extension to transient physics);

(3) Do the models reliably represent 'reality', within and beyond their calibration range? (validation).

In this work, we address questions (1) and (2). Question (3) cannot be addressed fully until (1) and (2) are resolved, and we leave (3) as a topic for future work by the glaciological community.

Here, we present the results of the Firn Model InterComparison Experiment (FirnMICE). We compared the steady-state and transient behavior of eight firn-densification models and one seasonal-snow model by forcing all models with the same suite of temperature and accumulation-rate boundary conditions. By running FirnMICE we seek to provide a fundamental understanding in the variability among firn models and to provide direction for future firn research.

Our goal is to compare results and differences among results that would be obtained if a nonmodeler were to ask a modeler peer to solve a posed question about firn, with each peer expert using his or her model with its own governing equations in the manner that he or she determined was most appropriate for the question. We are not assessing numerical implementations of governing equations in the individual models. Because each of these models has previously been published in peer-reviewed literature, we accept here that basic concerns about numerics have been adequately addressed by authors and peer reviewers. Comparing results from different numerical implementations of a specific set of model equations is a different question, worthy of future research.

This paper is not a review paper; we do not provide here an exhaustive review of firn observations and modeling work. FirnMICE was not designed to select a favored firn model among numerous options. At the current stage of firn-model development, an appropriate firn model for any particular study must be chosen based on the complexity required for the analysis and the temporal and/or spatial domain and climate conditions. The ice core and altimetry scientific communities often use different firn-densification models. The spatial domain and resolution are the same for these two applications, but the time steps are often different: firn-model applications for ice core studies use sparse data to model long-term (multidecadal and longer) changes in the firn, whereas firn-model applications for altimetry studies use daily or monthly weather data and outputs from regional-climate models to investigate seasonal and interannual changes in the firn. However, firn evolution is ultimately based on physical laws, and a model implementing those physics should be able to predict firn evolution on any timescale regardless of the model's application. A goal of the FirnMICE project is to identify firn-model commonalities and discrepancies and to encourage connections among firn-research applications in order to improve the physical descriptions and parameterizations.

As an extension of this idea, we recognize that seasonal-snow models also compute densification rates of ice/air mixtures, but often emphasize different processes (such as melting and percolation) and are optimized for much shorter temporal and spatial steps. This can lead to significantly different responses when they are applied to firn. However, a long-term challenge would be development of a unified evolution model that could be applied to both seasonal snow and firn.
## 2. EXPERIMENT DESIGN

### 2.1. Participating models

We compared eight firn-densification models, designated as follows: HLA and HLS (Herron and Langway, 1980), BAR (Barnola and others, 1991), GOU (Goujon and others, 2003), ART (Arthern and others, 2010), LIG (Ligtenberg and others, 2011), SIM (Simonsen and others, 2013) and CMG (Cummings and others, 2013). The Herron and Langway (1980) model has two solutions presented: an analytical solution (Herron and Langway, 1980, Eqns (7) and (10)) and a stress-based solution (Herron and Langway, 1980, Eqn (4c)). Features of the models and details about their implementation are shown in Table 1. Specific details about the models and their original applications are described in the Appendix.

The models have similarities; for example, all of the models use temperature and accumulation rate as the primary drivers of firn densification. Most of the participating models differentiate between zone 1 and zone 2 densification, and all of the models have an Arrhenius temperature dependence. However, the models differ in numerical implementation and description of firn-densification physics. The activation energies used in the Arrhenius term vary among models. The ART model includes a grain-size factor.

All of the models are empirical to some degree. The HLA, HLS, LIG, CMG, GOU, SIM, HLS and BAR models each have between four and seven adjustable parameters that have been tuned to match density-depth profiles from polar firn. The particular range of climatic conditions to which each model has been tuned is unique to that model (e.g. LIG was tuned using cores from 47 Antarctic sites; HLA was tuned with seven cores from Greenland and ten cores from Antarctica), although there is overlap in the sites used to tune (e.g. both LIG and HLA are tuned using data from South Pole). The ART model has only two adjustable parameters tuned to match strain rates measured in boreholes at sites that are relatively warm and have relatively high accumulation rate.

Firn-densification models are often categorized as either being ice core age models or altimetry mass-balance models, and we recognize the original application of each model in Table 1. The models may also be categorized by the physical description of the densification rate $Dp/Dt$ in a Lagrangian (material-following) coordinate system and by the assumptions built into the model physics. In Table 1, we describe how models treat stress in their densification-rate equations: some models implicitly represent stress through the accumulation rate, some models include the stress directly in their evolution equations and some parameterize the stress using the mean accumulation rate over the lifetime of a parcel of firn. The latter method is effectively the same as using the stress. None explicitly incorporates conservation of mass, momentum and energy. If the stress is parameterized by the immediate accumulation rate, the densification rate of a given parcel will have an instantaneous change when the accumulation rate is changed. If stress is included directly or parameterized by the mean accumulation rate over the lifetime of a parcel of firn, the densification rate of a parcel of firn will gradually adjust to the changing load history that accompanies an accumulation-rate change.

Groups of firn models share common ancestries. The Herron and Langway (1980) model serves as a benchmark

firn-densification model. Herron and Langway (1980) compiled data of density as a function of depth from firn cores from Antarctica and Greenland. They assumed that the depth/density relationship was invariant in time and derived densification-rate equations, tuning parameters to fit the measured density profiles. In the absence of sufficient strain-rate data from varied sites spanning a broad range of climatic conditions, most modeling efforts have continued to use the steady-state assumption. In the present study, the Barnola and others (1991) model explicitly uses the Herron and Langway (1980) equation for zone-1 densification. Arthern and others (2010) measured firn strain rates to derive a densification equation (the model included in this study), but also invoked a steady-state assumption to find updated coefficients for the Herron and Langway (1980) model. The Ligtenberg and others (2011) and Simonsen and others (2013) models are both extensions of the Arthern and others (2010) steady-state model, where the coefficients are updated to apply to a larger range of polar climates. Cummings and others (2013) use the Ligtenberg and others (2011) densification equation with an enthalpy formulation.

The Goujon and others (2003) model draws upon the physics described by Alley (1987) and Arnaud and others (2000). The models that are compared in this work do not constitute an exhaustive list of firn models. Notable exceptions include the models described by Alley (1987), Spencer and others (2001), Li and Zwally (2004), Helsen and others (2008), Li and Zwally (2011) and Kuipers Munneke and others (2015).

The ESS seasonal-snow model (Essery and others, 2013) has been included to see how well a physics-based seasonal-snow model can perform when applied to long timescales and surface conditions far outside its calibration range. The model incorporates a surface energy balance and stress-based densification, and itself represents the results of a wide-ranging intercomparison of seasonal-snow models.

### 2.2. Synthetic boundary conditions

We forced the models using synthetic temperature and accumulation-rate boundary conditions. There are a number of advantages to using this approach for this type of research question. It isolates the differences between empirical tuning and physics among models by examining steady state and transient model responses. The models show steady-state differences that are attributed to differences in empirical tuning and transient differences that are attributed to the way that the models describe firn-densification physics. The introduction of a climatic step change isolates effects of temperature and accumulation rate on each model's densification physics.

Future work could include incorporating temperature and accumulation rates for a particular site from reanalysis or ice-core reconstructions, but there are complications with using reconstructed data for the purpose of understanding fundamental firn physics. Firn models use empirical tuning to account for an imperfect description of physics. If an inter-comparison experiment were designed to force the set of firn models for reconstructed boundary conditions for the last 2000 years, the modeled depth-density profile at present could be compared with observed depth-density profiles for a particular site. Models empirically tuned to
<table><caption>**Table 1.** The participating models vary in densification physics, numerics, and application. See Appendix for more detail</caption><thead><tr><th>Characteristics</th><th colspan="5">Models</th></tr><tr><th>Model name (abbreviation used in figures)</th><th>Herron and Langway (1980): Analytic solution (HLA)</th><th>Herron and Langway (1980): Stress-based solution (HLS)</th><th>Barnola and others (1991) (BAR)</th><th>Goujon and others (2003) (GOU)</th><th>Essery, et al. 2013 (ESY)</th></tr></thead><tbody><tr><td>Numerical method</td><td>Analytical solution</td><td>Finite difference, Lagrangian densification, Eulerian heat transport</td><td>Finite difference, Lagrangian densification, Eulerian heat transport</td><td>Finite difference; Eulerian, with Lagrangian steady-state spin-up</td><td>Finite Difference, Eulerian</td></tr><tr><td>Grid</td><td>Evenly-spaced grid points (0.01 m) to 1000 m depth</td><td>0.25–1000 m depth</td><td>0.25–1000 m depth</td><td>0.5 m until 150 m; 5 m until 470 m; 50 m until the bottom (3465 m)</td><td>1000 m depth, 1007 nodes, layer thickness increasing exponentially from 0.01 m at surface to 1 m for depths of 1 m and below</td></tr><tr><td>Time step</td><td>N/A</td><td>0.2 years, explicit densification, Crank-Nicholson heat</td><td>0.2 years, explicit densification, Crank-Nicholson heat</td><td>1 year; explicit</td><td>1 day, Implicit for Heat Conduction, Explicit for Water Drainage</td></tr><tr><td>Spin-up/ Initialization</td><td>None</td><td>Initialized with H&amp;L steady-state solution; 10 000 years; 15 000 years for experiment 4</td><td>10 000 years; 15 000 years for experiment 4</td><td>Initialized with a steady state profile of temperature and density; then ran for 10 000 years</td><td>1000 years with constant surface temperature and accumulation</td></tr><tr><td>Heat flow</td><td>None</td><td>Diffusion-advection</td><td>Diffusion-advection</td><td>Diffusion-advection</td><td>Diffusion</td></tr><tr><td>Densification physics: load</td><td>Accumulation-rate parameterized</td><td>Stress</td><td>Stress</td><td>Stress</td><td>Stress</td></tr><tr><td>Densification physics: unique features</td><td>Benchmark model</td><td>Benchmark model</td><td>n = 3 exponent on stress</td><td>n = 3 exponent on stress; includes microphysical properties: particle contact area, coordination number, grain-boundary viscosity</td><td>non-Arrhenius temperature dependence</td></tr><tr><td>Texture evolution</td><td>None</td><td>None</td><td>None</td><td>None</td><td>Grain Size (used only for albedo parameterization)</td></tr><tr><td>Data tuned to</td><td>Seven sites in Greenland and ten sites in Antarctica</td><td>Seven sites in Greenland and ten sites in Antarctica</td><td>‘Several Antarctic and Greenland density profiles with temperatures ranging from –14°C to –57°C and accumulation rates from 2.2 to 65 g cm<sup>-2</sup> a<sup>-1</sup>’</td><td>Vostok, GISP2</td><td>Kojima, 1967</td></tr><tr><td>Tunable parameters</td><td>Six empirical parameters tuned to density-depth profiles</td><td>Six empirical parameters tuned to density-depth profiles</td><td>Seven empirical parameters tuned to density-depth profiles</td><td>Four empirical parameters tuned to density-depth profiles</td><td>3 empirical parameters tuned to density-depth profiles</td></tr><tr><td>Original model application</td><td>Glaciological interest</td><td>Glaciological interest</td><td>Ice core Δage</td><td>Ice core Δage</td><td>Surface Boundary Conditions for Atmospheric Model (seasonal snow model)</td></tr><tr><td>Reference</td><td>Herron and Langway (1980)</td><td>Herron and Langway (1980)</td><td>Barnola and others (1991)</td><td>Arnaud and others (2000), Goujon and others (2003)</td><td>Kojima (1967); Pitman and others (1992); Best and others (2011); Essery and others (2013)</td></tr></tbody></table>
<table><caption>**Table 1. Continued.**</caption><thead><tr><th>Characteristics</th><th colspan="4">Models</th></tr><tr><th>Model name (abbreviation used in figures)</th><th>Arthern and others (2010) (ART)</th><th>Ligtenberg and others (2011) (LIG)</th><th>Simonsen and others (2013) (SIM)</th><th>Cummings and others (2013) (CMG)</th></tr></thead><tbody><tr><td>Numerical method</td><td>Finite volume; irregular Lagrangian</td><td>Finite volume; Lagrangian</td><td>Finite difference, modular: Lagrangian or Eulerian</td><td>Finite element, Lagrangian</td></tr><tr><td>Grid</td><td>Approximate domain depth is 1000 m; Number of layers is fixed before the start of the simulation by dividing the approximate domain depth by the smallest ice-equivalent accumulation expected within a time step. The thickness of new layers is controlled by the mass accumulation rate, time step and density. One layer is removed from the base for each new layer added at the surface</td><td>100 layers; layer thickness 0.1 m</td><td>Determined by layer thickness, which is a function of the time step and accumulation. Layers are removed at a depth of 1000 m</td><td>Depth 1000 m, 1470 nodes; variable grid size with largest at the base (10 m) and smallest at the surface (0.078125 m)</td></tr><tr><td>Time step</td><td>0.1 years; implicit heat, explicit density</td><td>1 h; implicit heat, explicit density</td><td>Monthly time steps; modular: explicit or implicit</td><td>0.05; semi-implicit (Liniger) densification and velocity, Crank-Nicholson Enthalpy</td></tr><tr><td>Spin-up/ Initialization</td><td>(1) Run isothermal for enough time to advect all layers through firn; then, (2) run with temperature evolution for ≥1000 years</td><td>5000 years: complete firn column is refreshed with new snow</td><td>Initialized with analytical solution of H&amp;L; dynamic firn-densification model is run until steady state</td><td>5000 years, or enough time for new surface layer to reach 500 m depth</td></tr><tr><td>Heat flow</td><td>Diffusion-advection</td><td>Diffusion-advection</td><td>Diffusion-advection</td><td>Diffusion-advection</td></tr><tr><td>Densification physics: load</td><td>Stress</td><td>Mean-accumulation-rate parameterized</td><td>Mean-accumulation-rate parameterized</td><td>Mean-accumulation-rate parameterized</td></tr><tr><td>Densification physics: unique features</td><td>Viscosity of the ice matrix surrounding pores is assumed to be controlled by lattice diffusion creep (Coble, 1970); grain-size dependent; calibrated to measured dp/dt</td><td>Dependent on natural log of accumulation rate</td><td>Tuned using radar layers</td><td>Nabarro-Herring creep around pores; normal grain growth; enthalpy formulation</td></tr><tr><td>Texture evolution</td><td>Grain growth</td><td>None</td><td>None</td><td>Grain growth (Arthern and others, 2010, Eqn (4))</td></tr><tr><td>Data tuned to</td><td>Vertical strain rates recorded in boreholes at two sites on the Antarctic Peninsula and one site on Berkner Island</td><td>Founded on H&amp;L and Arthern and others (2010) and tuned these with observations used in Van den Broeke (2008)</td><td>Founded on H&amp;L (1980) and Arthern and others (2010), with additional tuning from two Greenland ice core sites (NGRIP and Flade Isblink) and a radar profile at the EGIG-line in Greenland</td><td>Founded on Arthern and others (2010) and Ligtenberg and others (2011)</td></tr><tr><td>Tunable parameters</td><td>Two empirical parameters tuned to measured borehold strain rates</td><td>Four empirical parameters tuned to density-depth profiles</td><td>Four empirical parameters tuned to density-depth profiles</td><td>Four empirical parameters tuned to density-depth profiles</td></tr><tr><td>Original model application</td><td>Ice-sheet mass balance</td><td>Ice-sheet mass balance</td><td>Ice-sheet mass balance</td><td>Ice sheet mass balance: Ice lenses and water transport in ablation zones</td></tr><tr><td>References</td><td>Arthern and others (2010)</td><td>Ligtenberg and others (2011)</td><td>Simonsen and others (2013)</td><td>Cummings and others (2013)</td></tr></tbody></table>
particular sites or regions would likely outperform other models (e.g. firn models tuned for Summit, Greenland are not necessarily expected to fit East-Antarctica data well), and the experiment might demonstrate which models have been empirically tuned for the site rather than which has the best representation of firn-densification physics. An experiment using site-specific reconstructed accumulation rates and temperatures must also consider the uncertainties in the reconstructed boundary conditions. If there is a model-output misfit to observed density profiles, is it an error in the reconstructed boundary conditions or a poor choice in the model physics? For this study we focused on investigating how differences among the models, including physics and empirical tuning, affect the firn density and age.

## 2.3. Methods

We compared the models in a series of six experiments by forcing the models with temperature and accumulation-rate step changes. FirnMICE participants completed model runs and submitted model output for the suite of experiments. The models were spun-up to an initial steady state using initial temperature and accumulation-rate boundary conditions specific to each experiment. Each FirnMICE participant was allowed to choose how specifically to spin-up his or her model; Table 1 includes a description of how participating models were initialized. After the spin-up, the model run began ($t = 0$ years). A step change in either temperature (Experiments 1–3; 5 K increase) or accumulation rate (Experiments 4–6; 0.05 m a⁻¹ increase) was prescribed at $t = 100$ years into the 2000-year model run. Figure 1 shows the initial conditions and step changes. The accumulation-rate and temperature boundary conditions are within the range of values observed in polar regions. The six

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

<figure><img src="image_3.png" /><figcaption>Fig. 1. The accumulation-rate and temperature boundary conditions specified for the six experiments. In Experiments 1–3 there was a positive 5 K temperature step change while accumulation rate remained constant at 0.1 m a⁻¹ (ice equiv.). In Experiments 4–6 there was a positive 0.05 m a⁻¹ accumulation rate step change while the temperature remained constant at -30°C.</figcaption></figure>

experiments were designed to (1) probe 12 steady-state model solutions, as the initial and final steady-state conditions for each of the experiments are unique, and (2) investigate the transient response of each model to the climate perturbation. The following constraints were also prescribed:

(1) The surface density was held constant at $\rho_0 = 360 \text{ kg m}^{-3}$.

(2) The Neumann (gradient) temperature boundary condition at the bottom of the firn was set to 0 K m⁻¹, ~1000 m below the surface. This large model domain (~10 times the firn column thickness) was chosen to account for the large thermal mass of the ice sheet.

(3) At high latitudes where surface temperature is driven by sunlight, the annual temperature cycle often shows a coreless winter, because the lowest portion of the sinusoidal insolation signal is effectively truncated in periods of total darkness. This effect can be further enhanced by higher windiness in winter; this can disrupt the boundary-layer temperature inversion more frequently in winter (Thompson, 1969).

Participants could choose to implement a seasonal temperature cycle $T_{seas}$ (K) that varied around the mean annual temperature. The recommended seasonal cycle (supplementary material in Orsi and others, 2012) was

$$T_{seas} = A(\cos(2\pi t) + 0.3\cos(4\pi t)), \quad (1)$$

where the amplitude of the seasonal cycle $A$ was given as 10 K, $t$ is time in years, and the second term generates a coreless winter. The optional addition of this cycle was included to allow participants to preserve authenticity of the models as they are typically used.

The initial grain size was not prescribed, nor was the parameterization for thermal conductivity. Participants were free to use their preferred values in their models.

## 2.4. Model output and derived quantities BCO and depth-integrated porosity (DIP)

Model outputs for the six experiments were firn age, density and temperature as functions of depth. From these outputs, it is possible to calculate the depth and age of lock-in and BCO as well as the DIP, which is the total amount of air in the firn column. The lock-in age is needed to calculate Δage in the past. Δage is the lock-in age less the age of the enclosed gas; however, the latter is very small (usually on the order of tens of years). For our synthetic modeling experiments in which we investigate differences among models, all participants calculated the depth and age of the 815 kg m⁻³ horizon, which we call the BCO horizon, following Barnola and others; if the results were to be used for gas-mobility studies, we would consider the firn age and depth at the lock-in density. DIP is used in studies of ice-sheet mass balance as a correction to altimetry measurements (e. g. Sandberg Sørensen and others, 2011). The DIP is the porosity $\phi$ integrated over depth $z$ from the surface to depth $z_i$ where ice density $\rho_i$ is reached:

$$DIP = \int_{0}^{z_i} \phi(z) dz = \int_{0}^{z_i} \frac{\rho_i - \rho(z)}{\rho_i} dz. \quad (2)$$
# 3. MODELING RESULTS

## 3.1. Steady-state firn-model variability

Because the six experiments start at steady state and end very near a new steady state (see Section 3.2), the DIP and BCO values from times $t = 0$ and $t = 2000$ years provide steady-state model results for 12 temperature and accumulation-rate pairs. Figure 2 shows steady-state (or near steady-state) values of DIP, BCO depth and BCO age for those 12 steady climates. Figures 2a-c show DIP, BCO depth and BCO age as a function of temperatures from -50°C to -25°C for a constant accumulation rate of $0.1 \text{ m a}^{-1}$. Figures 2d-f show DIP, BCO depth and BCO age as a function of accumulation rates from $0.02 \text{ m a}^{-1}$ to $0.30 \text{ m a}^{-1}$ for a constant temperature of -30°C. The standard deviations (SDs) provide a measure of the variation among the models for each calculation. The general pattern of variability includes relatively higher SD at the lower and upper bounds of the temperatures and accumulation rates in the experiment, because not all of the models have used data

from those extremes in their calibrations. Figure 2f shows the highest variability in BCO age at low accumulation rates, for which the SD is ~200 years at $0.02 \text{ m a}^{-1}$. This is attributable in large part to the ART model, which was calibrated only with warmer, higher-accumulation conditions in the Antarctic Peninsula and was not intended to apply to deeper, very-cold firn. The highest SD for depth is ~8 m at -20°C in Figures 2b, e.

The models' similar approach to tuning (to measured depth-density profiles) may explain some of the agreement in their steady-state responses. The ART model's lower sensitivity to accumulation rate may reflect differences in the number of adjustable parameters, the target used for tuning (warmer, higher accumulation sites), the representation of the underlying physical mechanisms, or in a combination of the three. Until the sensitivity of density-depth profiles to accumulation rate and temperature can be reproduced by a physically-based model without any adjustable parameters, there will be some uncertainty as to the physical explanation for this sensitivity. Further, this will translate into uncertainty

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** The FirnMICE models were at steady state at $t = 0$ after spin-up and near steady state at $t = 2000$ years (see discussion of transients in Section 3.2, providing 12 steady-state realizations of DIP and BCO). An accumulation rate of $0.1 \text{ m a}^{-1}$ and a range of temperatures produced the steady-state values for DIP and BCO shown in a-c. A temperature of -30°C and a range of accumulation rates produced the steady-state results for DIP and BCO shown in d-f. The SD among the eight models is shown in the lower panel for each DIP and BCO figure. See Section 2.1 and Appendix A for explanation of model legend.</figcaption></figure>
in the accuracy of the models when they are applied to cli-
matic conditions different from the tuning datasets, or to
modeling the response to temporal changes in climate. The
transient FirnMICE simulations are designed to quantify this
uncertainty in more detail.

Excluding the accumulation rate = 0.02 m a⁻¹ data point, 1 SD, shown in Figure 2, is ~2 m for DIP, <8 m for BCO-depth, and <60 years for BCO-age. The range of steady-state values given by FirnMICE models can be interpreted as the limit to which firn models can be tuned, such that this tuning can be extrapolated to other sites in a similar temperature and accumulation range. For example, DIP is unlikely to be inferred with confidence from a firn model to better than 2 m in steady state at a specific site with known accumulation rate and temperature. These limits are due to calibration differences, differences in the thermal coefficients in the advection-diffusion equation and other factors beyond temperature and accumulation rate that are not included in the models.

3.1.1. *Sensitivity to temperature*

Each of the FirnMICE models incorporates sensitivity to tem-
perature through an Arrhenius exponential factor with a con-
stant activation energy. Their activation energies for the
second stage range from 17.6 kJ mol⁻¹ (ART, LIG, in a
steady-temperature regime) to 60 kJ mol⁻¹ (BAR, GOU).
This Arrhenius relationship causes a nonlinear sensitivity to
temperature: the (absolute value of) slopes of steady-state
DIP, BCO depth and BCO age with respect to temperature
increase with decreasing temperature (Figs 2a-c), i.e. the
models predict a larger difference in steady-state DIP and
BCO values for a given difference in steady-state temperature
at colder steady-state temperatures. The models with lower
activation energies (HLA, ART, SIM, LIG) appear to have a
slightly higher sensitivity to temperature. From the form of
the model equations (see Appendix), however, we might
expect the differing activation energies among the models
to have a major influence on the temperature sensitivities.
However, that influence is likely counteracted by other
model factors, for example differing sensitivities to accumu-
lation rate.

It is not clear whether colder firn should have an increased sensitivity to temperature. Capron and others (2013) compared firn thickness in Antarctica predicted by the Goujon and others (2003) firn-densification model to ice core $\delta^{15}$N data, a proxy for firn thickness in the past. They found that the firn model failed to predict the observed $\delta^{15}$N variations during the last deglaciation and suggested that the sensitivity to temperature in firn models may be overestimated and sensitivity to accumulation rate may be underestimated. Zwally and Li (2002) suggested that activation energy is itself a function of temperature. The effective activation energy for firn densification also may depend on two or more activation energies for separate processes (e.g. see Arthern and others, 2010), which need to be determined independently.

3.1.2. *Sensitivity to accumulation rate* $\dot{b}$

FirnMICE models show an increase in DIP and BCO depth
with increasing accumulation rate $\dot{b}$. A higher $\dot{b}$ yields
younger firn at any given depth due to faster downward
advection. This younger firn has had less time to compact
and thus has lower density (higher porosity); as a result, the
close-off density is not achieved until a greater depth.

Because densification is driven by overburden load $\sigma$, densification rate at a given depth will be lower when the average density above that depth is lower due to higher accumulation rate. This reduced densification rate at a given depth, together with reduced time to compact before reaching that depth, further increases the sensitivity of model DIP and BCO depth to accumulation rate $\dot{b}$. This increased sensitivity is enhanced further in models with higher stress exponent. However, the BCO age decreases with increasing accumulation rate for all models, because the greater stress $\sigma$ due to greater overburden load at any age tends to produce faster densification so that close-off density is reached sooner.

These generalizations hold even though the models differ widely in their representation of the dependence of the densification rate on accumulation rate. Although ART qualitatively follows these trends, it is an outlier in (Figs 2d-e), showing less sensitivity to accumulation rate. ART is based on Nabarro-Herring creep, which has a linear stress dependence, and therefore might be expected to have less sensitivity to $\dot{b}$ than GOU and BAR, which are based on dislocation creep, with a cubic dependence on $\sigma$. Furthermore, densification is also opposed by larger grain size. Like stress $\sigma$, grain size $r^2$ also grows linearly with time, but from a nonzero starting point $r_0^2$. This further reduces the sensitivity to $\dot{b}$. As a result, ART is almost insensitive to the accumulation rate.

LIG used a steady-accumulation formulation of ART that includes a linear dependence on $\dot{b}$ (see Appendix), and introduced a dependence on $-\ln(\dot{b})$. CMG follows the same equation as LIG and has a very similar response. GOU and BAR are based on dislocation creep, with a dependence on $\sigma^3$, and interestingly give similar results to LIG and CMG. HLS and SIM are the most sensitive to the accumulation rate, and for these models the densification rate is proportional to $\sqrt{\dot{b}}$.

3.1.3. *Challenges constraining temperature and accumulation-rate sensitivities*

It is difficult to accurately constrain the stress dependence
and activation energy with modern field data because in
practice, accumulation and temperature are almost always
correlated (accumulation rate increases with temperature)
and increasing temperature and increasing accumulation
rate have opposite effects on firn-densification rate.
Improved constraints on activation energy and stress depend-
ence may come from controlled experiments on deformation
of firn in a laboratory, where the temperature and loading
history can be isolated. Progress has been made on this
problem for snow using micro-computed tomography (e.g.
Schleef and Löwe, 2013). However, under typical stress
levels found in nature, the long timescales for firn densifica-
tion may prove to be challenging for laboratory experiments.
Additionally, setting the rate of incremental loading would be
very difficult. Therefore, field investigations will also be
needed at targeted suites of sites where accumulation
varies significantly while temperature does not, and vice
versa. For example, on a 20-km south-north transect across
Taylor Dome, temperature changes by only 3°C
(Waddington and Morse, 1994), while accumulation rate
varies by a factor of four (Morse and others, 1999). That
site could isolate accumulation-rate effects. Identifying
other sites or suites of sites with similar accumulation rates
but substantially different temperatures could help isolate
temperature effects.

**3.2. Transient firn-model variability**

3.2.1. *Density variability*

Figure 3 shows each model's predicted depth-density profile at t = 2000 years and the mean density profile of all the models $\bar{\rho}(z, t = 2000)$ for Experiment 1. In Figure 4, we characterize the differences among the models' predicted depth-density profiles using the $SD_p(z, t)$ of the models' predicted densities $\rho(z, t)$ at selected times t. Because $\bar{\rho}(z, t)$ increases monotonically with depth z, we use $\bar{\rho}(z, t)$ as a proxy for depth to produce the same independent variable for all six experiments. The pattern of density variability is qualitatively the same through time for the six experiments, shown at 0, 150, 250 and 2000 years. The density variability is small at the surface where the models are constrained by the prescribed initial condition and at depth where they all approach ice density.

A local minimum in density variability among the models
also occurs at ~600 kg m⁻³, which is near the density com-
monly accepted as the transition from zone 1 to zone
2. This minimum is a result of the variance in the models'
densification rates in zone 1 and zone 2, as can be seen in
Figure 3. Compared with the other models, the GOU
model predicts a higher densification rate as a function of
depth in zone 1 and a transition to zone 2 densification at
a lower density. As a result, SDₚ(z, t) increases in zone 1 as
GOU diverges from the other models. Below the zone 1/
zone 2 transition, SDₚ(z, t) decreases as GOU predicts a

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

**Fig. 3.** Depth-density profiles and the mean depth-density profile $\bar{\rho}(z,t)$ for Experiment 1 at the end ($t = 2000$ years) of the model run. GOU predicts a faster densification rate in zone 1, and its transition to zone 2 happens at a shallower depth and lower density than the other models (see Appendix). This behavior creates the first maximum in standard deviation $SD_p(z,t)$ seen in **Figure 4**. Deeper in the firn, the models' predicted depth-density profiles diverge, leading to the second $SD_p(z,t)$ maximum seen in **Figure 4**. The figure is representative of all the experiments. See Section 2.1 and Appendix A for explanation of model legend.

lower densification rate as a function of depth in the upper
part of zone 2 and its predicted depth-density profile
approaches the mean depth-density profile. Deeper in the
firn, variance in the models' predicted densification rates
causes the modeled depth-density profiles to diverge, and
$SD_p(z, t)$ increases to a maximum between the densities of
$750$ and $850 \text{ kg m}^{-3}$ before the models converge to the ice
density at depth. Although Figure 3 only shows results from
Experiment 1, it is representative of all the experiments:
GOU consistently predicts a higher densification rate in
zone 1, but no single model consistently predicts higher or
lower densification rates in zone 2.

3.2.2. *Temperature variability*

Temperature profiles of the model results for Experiment 1 are shown at times 0, 150 and 2000 years in Figure 5. At initialization, the models have (mostly) reached steady state, with nearly constant and uniform temperature of -50°C. At time t = 150 years, the temperature profiles have been responding for 50 years to the 5°C increase in temperature. The Herron and Langway (1980) analytic model immediately takes the new temperature of -45°C and the CMG model is slowest to respond, due to differences in the parameterization of the diffusion-advection equation. The time-dependent temperature models differ at 150 and 2000 years, reflecting different tuning in heat-equation terms, including the thermal conductivity, density, and specific-heat terms of the heat equation. At 2000 years the temperature results show that a new thermal steady state has not yet been reached; however, we can treat density profiles as approximately steady state, and BCO and DIP, two measures of interest, undergo only small differences beyond 2000 years.

The time to reach thermal equilibrium after a surface-
temperature change (Experiments 1, 2 and 3) is longer than
the 2000 years of the model run. This (diffusive) timescale
is a result of the large thermal mass of the underlying ice
sheet, which introduces a strong memory effect. We can esti-
mate the diffusive e-folding timescale $\tau$ to 200 m depth for
the temperature perturbation by nondimensionalizing the
heat equation: $\tau \approx z^2/\kappa$. Approximating the thermal diffusivity
of firn $\kappa = 8 \times 10^{-7} \text{ m}^2 \text{s}^{-1}$ (Giese and Hawley, 2015):
$\tau \approx (200 \text{ m})^2/(8 \times 10^{-7} \text{ m}^2 \text{/s}) \approx 1600$ years. The advective
heat-transport timescale is of similar magnitude: it can be
approximated as $\tau \approx z/\dot{b} = 200 \text{ m}/0.1 \text{ ma}^{-1} = 2000$ years.

**3.3. DIP results**

The DIP (Eqn (2)), which is the total amount of air in the firn
column, is the key model result needed for studies of ice-
sheet mass balance in order to correct observations of
surface-elevation change for changes in firn-air content. In
Figure 6, DIP is shown through time for the six experiments.
The DIP varies among models in both the steady-state and
transient behavior. First, we note that following spin-up (i.e.
during the first 100 years), the DIPs calculated by the
various models vary by a range of typically 4–8 m of air in
the firn column. These differences among the models
represent a large uncertainty in the total amount of air in
the firn column, which is typically between 20 and 40 m.

The transient behavior (i.e. after the temperature or accu-
mulation rate perturbation 100 years into the run) generally
results in a decrease in DIP with an increase in temperature
(Figs 6a-c), and an increase in DIP with an increase in
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** The variation of standard deviation SD<sub>p</sub>(z) among the models at depth z where mean density is ρ̄(z) is shown at times t = 0, 150, 250 and 2000 years. Using ρ̄(z) as a depth proxy allows us to show results from all six experiments on a common independent variable. The depth pattern of density variability among the models is maintained through time for the six experiments, with a minimum in variation at ~600 kg m<sup>-3</sup> and maxima at 450 kg m<sup>-3</sup> and at 750–850 kg m<sup>-3</sup>.</figcaption></figure>

accumulation rate (Figs 6d–f), as noted in Section 3.1. After the transients have run their course, the models generally tend to settle back into the same relative order of DIP; however, the amount of change varies among the models, shown in Figure 7. These differences among models indicate the range of uncertainty, which should be taken into account when calculating the mass change in an ice sheet in response to climate changes.

The majority of the models tend to relax monotonically toward their new DIP states over centennial timescales. This is to be expected for models that depend only on climate states represented by surface temperature *T* and accumulation rate *b*, because their evolution depends explicitly only on firn temperature and overburden load. However, in recognition that firn structural properties affect densification rate, models based on Arthern and others

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Temperature profiles from the firn models in Experiment 1 are shown at 0, 150 and 2000 years. The results have not reached steady state at 2000 years. Differences among the models arise both from differing rates of densification and differences in the firn-thermal-diffusivity parameterization. See Section 2.1 and Appendix A for explanation of model legend.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** The DIP for six experiments from participating models. The time is shown to 1000 years to highlight variation after the step change at 100 years. (a–c): For all models in steady state (t = 0–100 years), DIP is smaller when the mean annual temperature is greater, and the step-change temperature increase causes a decrease in DIP. (d–f): DIP is larger when the accumulation rate is higher, and the accumulation-rate increase causes an increase in DIP. The magnitude of the change is larger when the climate perturbation is proportionately larger (e.g. Experiment 4 is a 250% increase in accumulation rate, while Experiment 6 is just a 20% increase). See Section 2.1 and Appendix A for explanation of model legend.</figcaption></figure>

(2010) either explicitly (ART) or implicitly (LIG, SIM, CMG) incorporate the additional physical process of grain growth. These additional processes can influence densification rates over shorter timescales because aggregates of larger grains compact more slowly in those models. This effect can modify or even reverse the changes in DIP over time, producing histories that are more complex. Output from LIG shows an initial increase in DIP following the step increase in temperature (Figs 6a–c). However, DIP ultimately decreases, following the general trend of the other models. This pattern could be explained by initial stiffening of the near-surface firn due to larger snow grains or enhanced grain growth initially exceeding the ultimate softening due to warmer temperatures throughout the firn column. This is not the case for the Simonsen and others (2013) model (SIM), although they share the same temperature dependency with the same coefficients in their Arrhenius term (see Appendix). SIM does, however, include an additional Arrhenius term in its zone-2 densification equation, which could explain the different transient responses to temperature.

Unlike in the other models, which all show only a monotonically increase in DIP driven by enhanced burial rate of lower-density firn following the step increase in accumulation rate (Figs 6d–f), the initial increase in DIP in the ART and SIM models is then followed by a gradual decrease. This 'overshoot' behavior occurs due to the different

timescales of grain growth and advection. The increase in accumulation rate causes the firn to thicken as more air is incorporated (see Fig. 10); the new climate effectively adds snow to the surface more rapidly, and that snow has less time to consolidate before reaching any given depth; just as in the simpler models, this factor initially increases the DIP. The increased accumulation rate also increases the downward advection rate of a parcel of firn relative to the transient surface, meaning that a firn parcel at a given depth is younger than a parcel at the same depth before the accumulation increase. Grain growth is a time-dependent firn-stiffening mechanism, and so the younger firn is softer. Eventually, the increased softness of the firn becomes the dominant process as the older, stiffer firn is replaced by the younger, softer firn. The softer firn column then compacts more rapidly at any depth when compared with times before the climate step, and the DIP decreases again. However, as in the other models, DIP ultimately still reaches a greater value than it had initially.

Figure 8 shows the time derivative of the DIP, i.e. the slopes of the curves in Figure 6, illustrating the rate at which firn responds to the climate changes with each model. The models generally predict that DIP changes quickly after the climate step change, but the magnitudes of the predicted rates of change vary. The rate of DIP change decreases through time, but not necessarily monotonically.
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** As in Figure 6, but each model has its initial steady-state value subtracted from the time series to highlight the variation in model-predicted DIP change that results from the climate step change. See Section 2.1 and Appendix A for explanation of model legend.</figcaption></figure>

## 3.4. BCO results

BCO age and depth are key model results for ice core paleoclimate studies. The BCO age and depth are shown for the six synthetic experiments in Figures 9 and 10. The BCO age generally decreases with the increase in temperature and accumulation rate. The BCO depth decreases with an increase in temperature due to the increase in densification rate; however, the BCO depth increases with an increasing accumulation rate since the increase in advection rate outweighs the increase in densification due to the overburden load.

As with the model-predicted DIP, the transient BCO responses vary considerably among the models. In the temperature-step-change experiments (1–3) there is disagreement in the sign of the BCO age response immediately following the step change (Figs 9a–c). Additionally, ART shows an initial decrease in BCO depth, followed by an abrupt increase and then again a decrease towards the new steady state (Figs 10a–c), as the grain-size and viscosity continue to adjust after the step change. The time the models take to reach equilibrium with a new accumulation rate (experiments 4–6) is relatively short and approximately equal to the BCO age. This is the (advective) timescale needed to refresh the firn column with snow deposited under the new climatic conditions. Once the entire firn column has been refreshed down to the BCO, no memory remains of the previous climatic state. Note that the DIP response time (Fig. 6) is longer, because it involves density changes in the entire firn column and not just in the upper firn above the BCO depth. The magnitudes of the predicted changes in the BCO metrics vary as well, similarly to the DIP changes shown in Figure 7.

In Experiment 4, the ART and SIM models show sudden discontinuous jumps in BCO age (Fig. 9d) and depth (Fig. 10d).

These models incorporate grain-size dependence such that presence of large grains slows densification. Snow deposited after the accumulation-rate step change (an increase by a factor of 3.5) can reach any depth faster than snow deposited before the step increase, and at that depth the snow has a history of smaller grains and lower effective viscosity. It has therefore experienced significantly more densification and achieved higher density. For some time following the step change, the firn deposited after the step change can achieve a density higher than the older firn below it that was deposited before the step change. During that time span, BCO is still achieved only in the pre-step-change firn, but a density reversal can develop in the density-depth profile (i.e. higher-density firn shallower than lower-density firn). When the first layer of post-step-change firn reaches the close-off density, it creates a second, shallower BCO transition and air exchange between the atmosphere and the deeper pre-step-change firn is abruptly blocked. Until the last pre-step-change firn achieves BCO, there are two BCO transitions, a deep BCO in the pre-step-change firn and a shallow BCO in the post-change firn. However, only the shallower BCO is tracked, because it is the important one for air exchange and bubble occlusion.

# 4. DISCUSSION

## 4.1 Model variability and extension to transient physics

We address questions (1), (2) and (3) posed in Section 1.

(1) When models are based on similar physical concepts, do their results agree? (model variance) Agreement can be
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** The time derivative of the depth-integrated porosity (dDIP/dt) for participating models in the six experiments. (a–c): Experiments 1–3, step increase in temperature. (d–f): Experiments 4–6, step increase in accumulation rate. Most models exhibit monotonic relaxation to the new steady state. Models with additional physics can show more complex adjustment patterns in the initial decades to centuries. In (a–c): see LIG. In (d–f): see ART. See Section 2.1 and Appendix A for explanation of model legend.</figcaption></figure>

assessed with measures such as density variation with depth and age, and with measures such as DIP and BCO, which are integrated through the firn column. We recognize that experiments 1–6 explore subtle differences in model variation, where a small change in density $ρ(z)$ can produce a large change in DIP and BCO. Figure 2 shows variability in BCO and DIP among models for steady-state temperature and accumulation rate. Figure 4 shows that the models are more tightly clustered at some densities than others. Aside from Experiment 4 (which is a climate likely outside of the models’ calibration ranges, and was designed to push the models to their limits), the models generally show similar SDs, $SD_p(z)$ at similar mean densities $\bar{\rho}(z)$ in zone 1 ($SD_p(z) \sim 30 \text{ kgm}^{-3}$ at $\bar{\rho}(z) \sim 450 \text{ kgm}^{-3}$) and zone 2 ($SD_p(z) \sim 10 - 20 \text{ kgm}^{-3}$ at $\bar{\rho}(z) \sim 700 - 800 \text{ kgm}^{-3}$). The peak SD is then ~7% of the mean.

Figures 9 and 10 show significant steady-state variation at $t = 0$ and $t = 2000$. Transient variation is most notably shown in time series for DIP and BCO, in Figures 6, 9 and 10. Here we see disagreement in the shape of the response signal and response timescale for step changes in surface conditions. Figure 8 shows how this disagreement among physical models could even lead to conflicting estimates of the sign of mass change in an ice sheet. General agreement among the models includes that increased temperature or accumulation rate increases the densification rate; however, the large variation in the DIP and BCO has the

ability to change interpretations for mass balance and paleoclimate analysis from ice core data.

(2) When models are calibrated under steady-state conditions, do they produce similar histories of change under transient conditions? (extension to transient physics); The Herron and Langway (1980) steady-state equation is extended to transient dynamics in different ways. A variety of empirical tuning options can be used to fit model results to measured (or assumed) steady-state profiles, including adjusting coefficients and exponents in different mathematical terms that are loosely tied to physical mechanisms. We see that these models may effectively overfit steady-state results, where they disagree in the generalization to transient cases.

(3) Do the models reliably represent reality, within and beyond their calibration range? (validation) We leave the topic of validation of models to future work. We have purposely avoided awarding a single model the title of 'best' since most are tuned to particular sites and in choosing any single site we would favor an arbitrary 'best' model. We encourage work developing physics-based models that are supported by in-situ field studies targeting formulation of constitutive relationships.

## 4.2. Parameter dependence and nonuniqueness

In Section 1 we discussed the difficulty of separating the influences of site accumulation rate and temperature when
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** The BCO age of the six experiments for the participating model. (a–c): The BCO age is younger with increased temperature and (d–f) increased accumulation rate. See Section 2.1 and Appendix A for explanation of model legend.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** The BCO depth for the six experiments for participating models. (a–c): The BCO depth decreases with increased temperature and (d–f): increases with increased accumulation rate. See Section 2.1 and Appendix A for explanation of model legend.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Snow models and firn models both compute the densification of ice/air mixtures, but are calibrated under very different conditions in which different processes can be dominant. Here we compare depth-density profiles from the Essery and others (2013) snow model (ESS, blue) and from the Herron and Langway (1980) model (HLA, red) for the six firn experiments. Dashed curves are initial steady states, and solid curves are final states.</figcaption></figure>

empirically tuning models to steady-state firn-profile data. Inferring model parameters from a dataset is an inverse problem (e.g. Aster and others, 2005). When models are calibrated with incomplete data over a limited range, nonuniqueness is difficult to avoid, in that multiple combinations some model parameters can all allow output from a model to match a dataset at an acceptable tolerance level.

When transients are considered, additional confounding interactions can arise because model parameters are not independent. Additional trade-offs are possible. For example different values of thermal conductivity can produce different transient temperature profiles, but could still match transient data from depth-density profiles with a different value of the activation energy, which can stiffen or soften the firn to compensate.

As in any inverse problem, the way to reduce nonuniqueness is to incorporate additional model parameters that cause the model to respond in more specific ways to specific forcings or to specific data and to acquire additional data that can discriminate more effectively among the forcing influences. This is additional motivation for future physics-based models that can specifically incorporate more processes.

## 4.3. Firn models and snow models

Firn-evolution models clearly share some elements with models of seasonal-snow evolution, since both deal with

densification of ice/air mixtures. There is a large literature on snow models (see Essery and others (2013) and references therein), because they have been applied widely in global climate models. Because seasonal snow changes rapidly in response to changing weather and seasonal conditions, snow models are already necessarily strongly based on physical processes, with some empirical parameters. In an inter-comparison study of snow models, Essery and others (2013) found an unexpectedly wide variance in their responses to common forcing based on data from a site in the French Alps, and presented recommendations for reducing that variation.

Because we anticipate that future firn models will become more directly process-based, we have included the process-based Essery snow model (ESS) in the six experiments. Seasonal snow by definition lasts <12 months, is often on the order of 1 m deep, and seldom develops densities approaching the ice density, while firn profiles develop over centuries and to depths of typically 100 m. Therefore, the six FirnMICE experiments present a severe extrapolation challenge for snow models. Figure 11 shows the initial and final depth-density profiles for the Essery and others (2013) model (ESS) and the Herron and Langway (1980) model (HLA). The snow model produces profiles that densify too rapidly with depth, because physical processes that inhibit densification at high densities are incomplete in or missing from snow models; these processes are not needed in the typical calibration ranges of snow models.
Development of process-based firn models is not as advanced as the development of snow models, and we can be confident that current firn models would do a poor job of simulating seasonal snow cover. In the future, a model that can accurately simulate ice/air mixtures over the complete spectrum of time, depth and density found in nature will need to incorporate the physics currently in snow models or even empirical fits to simplified models of those processes, while adding additional processes that become dominant at higher densities, longer times and colder temperatures.

## 4.4. Moving beyond the current generation of empirical firn-densification models

In reality, firn-densification rate is affected by the microstructural characteristics and processes occurring at the grain scale, including density layering, hoar, nonuniformity of grain shapes and sizes, impurity loading and crystal orientation and bonding, yet the current generation of densification models has a limited representation of these processes and most ignore microstructure. Most of the firn models included in this work are empirical, with constitutive relations based on simplified, idealized physics, tuned to fit depth-density profiles that are assumed to be in steady state. An important goal in firn research is to improve our understanding of the relevant firn physics on the scale of individual grains, with the aim of building a 'next-generation' densification model; however, a full model description of all the relevant physics at the microstructure level will necessarily have more free parameters than can be constrained using field observations. When simulating firn evolution on large spatial or temporal scales, these data are unavailable. In the interior of the Antarctic continent surface precipitation is transformed into ice over a period on the order of several thousand years. Observational programs on these timescales are not possible; therefore, some degree of empiricism is likely to be unavoidable for practical applications, particularly when studying firn with small densification rate or with large spatial/temporal scales.

The hope is that improved understanding of grain-scale densification physics will lead to more accurate constitutive relationships. Arthern and others (2010) measured in-situ densification rates in Antarctic firn and used those data to calibrate a firn model that includes grain growth and temperature seasonality. Our comparisons highlight the need for time-dependent observational constraints that move beyond temperature and density profiles that are assumed to be in steady state. Constraints on transient behavior of depth-density profiles and microstructural profiles could come from modern observations and also from ice core records of abrupt climatic variations; an example is the WAIS Divide core, where a pronounced accumulation anomaly ~12 ka BP is accompanied by a response in $δ^{15}$N, a proxy for firn thickness in the past(WAIS Divide Project Members, 2015).

Over recent decades, other strategies have been developed for dealing with model shortcomings. First, in reconstructing past Δage in Greenland ice cores, constraints are available during abrupt Dansgaard-Oeschger (D-O) climatic variations. The D-O events are recorded both in the ice phase ($δ^{18}$O of water) and in the gas phase (CH$_4$, thermal fractionation in $δ^{15}$N of N$_2$); comparing the two provides a

direct Δage estimate that is independent of firn models (Leuenberger and others, 1999; Rasmussen and others, 2013). Moreover, the present-day Δage can be estimated very accurately using firn air sampling experiments (Schwander and others, 1993; Battle and others, 2011; Buizert and others, 2012). The firn model is then forced to fit the data-based constraints, thereby minimizing potential model biases.

Another approach to firn modeling involves firn model ensembles of model physics and coefficients (Buizert and others, 2014, 2015). Parrenin and others (2012) used estimates of firn-column thickness derived from $δ^{15}$N data together with an ice-flow model to reconstruct ice core Δdepth (the difference in the depths of ice and trapped air of equivalent age), rather than Δage.

The three firn models that have historically been used for Δage reconstructions are the Barnola and others (1991) model, the Goujon and others (2003) model and the dynamical Herron and Langway (1980) model (Table 1). None of these three models include microstructure evolution, as is used in the Arthern-based models (ART, LIG, SIM). Comparison of the transient responses (Figs 9 and 10) suggests that when grain evolution is included in the model it does impact short-term model behavior. As such, an investigation of the effects of grain evolution on Δage reconstructions is warranted, particularly in Greenland where transient effects are expected to be most important given the evidence of past rapid climate changes (Severinghaus and Brook, 1999).

The application of densification models to ice-sheet altimetry (e.g. Arthern and Wingham, 1998; Zwally and Li, 2002) is a more recent development than Δage reconstruction (Barnola and others, 1991). Notable recent developments in the application to altimetry include empirical corrections that allow the models to be used over a wider range of climatic conditions (e.g. Ligtenberg and others, 2011) and inclusion of additional physical processes, such as meltwater percolation and refreezing (e.g. Reeh, 2008; Ligtenberg and others, 2011). New experimental techniques have allowed the measurement of strain rates in the upper layers of the ice sheets (Arthern and others, 2010; Morris and Wingham, 2011; Hawley and Waddington, 2011). These methods provide data that can be directly related to densification rates without assuming that a firn column is in steady state, and they allow a more direct exploration of the effects of temperature on the rate of densification. For the purposes of correcting satellite altimetry on continental scales, models of firn densification are usually driven by climatic information derived from regional climate models (Ligtenberg and others, 2011) or from remote sensing (Zwally and others, 2005).

Differences in the evolution of DIP are significant for mass-balance modeling from altimetry, where firn densification is the largest source of uncertainty (Shepherd and others, 2012). The disagreement among the models in the sign of the immediate firn response to a climate change could cause misleading inferences of ice-sheet mass changes derived from altimetry data. In real-world applications, surface elevation change $dh/dt$ predictions from firn models furthermore depend strongly on the quality of input data (temperature, accumulation rate), which is not investigated here. Altimetric applications of firn models have not usually explored the effect of using different models of firn densification on their results, or different forcing data, but the results of this study show that significant differences can arise
depending upon the choice of model. This indicates that ensembles that bracket a range of different models, a range of different parameter choices and a range of forcing data might prove just as useful in altimetric applications as they have in studies of Δage reconstruction.

The strategies identified here for moving beyond the shortcomings of the current generation of firn-densification models (i.e. developing a next-generation physically-informed model, and increased use of model ensembles and data constraints) will have to be pursued in parallel. Developing the next-generation model requires time-dependent strain-rate observations and in the meantime important scientific questions in ice-sheet altimetry and ice core science must be addressed. We believe that the FirnMICE project both provides a motivation to increase efforts in firn-model development, as well as informs the glaciological community about the existing model shortcomings in steady-state and transient behavior.

## 5. CONCLUSIONS

The FirnMICE experiments consist of step changes in boundary conditions; however, this work is concerned with the variation among model output in a broader sense than any particular combination of boundary conditions or numerical implementations of governing equations. Our goal is to explore the range of results that a researcher could expect when collaborating with a range of modeler colleagues, each of whom runs a respective model at its optimum capabilities for a given problem.

Empirically-tuned models included in this study can match firn-density and strain-rate data reasonably well in the region(s) for which they were calibrated, but we have shown that the models do not agree with each other well in a steady state using the metrics of DIP, BCO depth and BCO age, which are in most cases the desired quantities for which the models were developed. The temperature sensitivities of the models are similar because there is consensus among the models that firn densification has an Arrhenius temperature dependence, and differences in activation energy have been compensated by differences in other empirical factors. The accumulation-rate sensitivities show more variance, indicating that a similar consensus on the relation between stress and strain rate has not been achieved.

The transient behavior of the models also shows variance; the general trends of change in the firn are consistent after the step changes, but the models disagree on the rates of change and the magnitudes of short-term and long-term changes. Examining the range of transient results here does not necessarily give insight into transient firn evolution: most of the models are calibrated using depth-density profiles that are assumed to be steady-state, and transients are not considered in calibration. It may not be appropriate to extend these steady-state-based models to transient evolution.

We have avoided choosing a 'best' firn model; the best-firn model for one particular application may be different than the best model for another application. We used the mean of the nine models as a metric to compare the models, but the mean of the set of models is unlikely to be the 'best' or 'true' solution. We used synthetic boundary conditions specifically to identify the variation among the models rather than testing which model can best match the data.

This work underscores the need to develop a 'next-generation' firn model that includes physically-based constitutive

relationships rather than an empirical tuning to match measured firn-density profiles. As ice core and altimetry data continue to be sampled at higher temporal resolution, firn models need to improve to be able to predict past Δage and to convert surface elevation change to mass change at the accuracies demanded by those data. Strain-rate data offer the potential to constrain transient responses in densification, and concurrent measurements of firn microstructure and chemistry offer an opportunity to better understand the constitutive relationship. These data are currently sparse, but progress is being made (e.g. Arthern and others, 2010; Freitag and others, 2013; Morris and Wingham, 2014; Proksch and others, 2015) to expand these datasets and to better understand firn evolution.

## ACKNOWLEDGEMENTS

We are grateful to Jeff Severinghaus and Ed Brook for inspiring the idea for this intercomparison. We acknowledge Michiel van den Broeke's leadership and contributions to IMAU firn modeling. This work was supported by National Science Foundation Partnerships in International Research and Education (PIRE) grant 0968391.

## CONTRIBUTION STATEMENT

J. Lundin designed the experiment, produced model output and wrote the manuscript. C. M. Stevens wrote and edited the manuscript, processed data and produced figures. C. Buizert and A. Orsi designed the experiment, produced model output and contributed feedback to the manuscript. R. Arthern, S. Ligtenberg, S. Simonsen, R. Essery and E. Cummings produced model output and contributed feedback to the manuscript. W. Leahy processed the data and produced figures. P. Harris developed model code. M. Helsen contributed feedback to the manuscript. E. Waddington wrote and edited the manuscript.

## REFERENCES

* Alley RB (1987) Firn densification by grain-boundary sliding: a first model. *J. Phys.*, **48**(C1), C1–C1249

* Arnaud L, Barnola J-M and Duval P (2000) Physical modeling of the densification of snow/firn and ice in the upper part of polar ice sheets. International Symposium on Physics of Ice Core Records. Shikotsukohan, Hokkaido, Japan, September 14-17, 1998. *Phys. Ice Core Rec.*, 285–305

* Arthern RJ and Wingham DJ (1998) The natural fluctuations of firn densification and their effect on the geodetic determination of ice sheet mass balance. *Clim. Change*, **40**(3–4), 605–624

* Arthern RJ, Vaughan DG, Rankin AM, Mulvaney R and Thomas ER (2010) In situ measurements of Antarctic snow compaction compared with predictions of models. *J. Geophys. Res.*, **115**(F3)

* Aster RC, Borchers B and Thurber C (2005) *Parameter estimation and inverse problems*. Elsevier, Amsterdam

* Bader H (1954) Sorge's Law of densification of snow on high polar glaciers. *J. Glaciol.*, **2**(15), 319–322

* Barnola J-M, Pimienta P, Raynaud D and Korotkevich YS (1991) CO₂-climate relationship as deduced from the Vostok ice core: a re-examination based on new measurements and on a re-evaluation of the air dating. *Tellus B*, **43**(2), 83–90

* Battle MO and 8 others (2011) Controls on the movement and composition of firn air at the West Antarctic Ice Sheet Divide. *Atmos. Chem. Phys. Discuss.*, **11**(6), 18633–18675
Best MJ and 16 others (2011) The Joint UK Land Environment Simulator (JULES), model description – Part 1: Energy and water fluxes. *Geoscientific Model Development*, **4**(3), 677–699

Buizert C and Severinghaus JP (2016) Dispersion in deep polar firn driven by synoptic-scale surface pressure variability. *Cryosphere*, **10**(5), 2099–2111

Buizert C and 25 others (2012) Gas transport in firn: multiple-tracer characterisation and model intercomparison for NEEM, Northern Greenland. *Atmos. Chem. Phys.*, **12**(9), 4259–4277

Buizert C and 9 others (2014) Greenland temperature response to climate forcing during the last deglaciation. *Science*, **345** (6201), 1177–1180

Buizert C and 9 others (2015) The WAIS Divide deep ice core WD2014 chronology–Part 1: Methane synchronization (68–31 ka BP) and the gas age–ice age difference. *Clim. Past*, **11**(2), 153–173

Capron E and 12 others (2013) Glacial-interglacial dynamics of Antarctic firn columns: comparison between simulations and ice core air-δ¹⁵N measurements. *Clim. Past*, **9**(3), 983–999

Coble Robert L (1970) Diffusion models for hot pressing with surface energy and pressure effects as driving forces. *J. Appl. Phys.*, **41**, 4798

Cuffey KM and Paterson WSB (2010) *The Physics of Glaciers*, 4th edn. Butterworth-Heinemann/Elsevier, Burlington, MA

Cummings E, Johnson J and Brinkerhoff D (2013) Development of a finite element firn densification model for converting volume changes to mass changes. *ArXiv e-prints* (https://arxiv.org/abs/1308.6616)

Essery R, Morin S, Lejeune Y and Ménard CB (2013) A comparison of 1701 snow models using observations from an alpine site. *Adv. Water Resour.*, **55**, 131–148

Freitag J, Kipfstuhl S, Laepple T and Wilhelms F (2013) Impurity-controlled densification: a new model for stratified polar firn. *J. Glaciol.*, **59**(218), 7

Gagliardini O and Meyssonnier J (1997) Flow simulation of a firn-covered cold glacier. *Ann. Glaciol.*, **24**, 242–248

Giese AL and Hawley RL (2015) Reconstructing thermal properties of firn at Summit, Greenland, from a temperature profile time series. *J. Glaciol.*, **61**(227), 503–510

Goujon C, Barnola J-M and Ritz C (2003) Modeling the densification of polar firn including heat diffusion: application to close-off characteristics and gas isotopic fractionation for Antarctica and Greenland sites. *J. Geophys. Res.*, **108**(D24)

Gow AJ, Meese DA and Bialas RW (2004) Accumulation variability, density profiles and crystal growth trends in ITASE firn and ice cores from West Antarctica. *Ann. Glaciol.*, **39**, 101–109

Hawley RL and Waddington ED (2011) In situ measurements of firn compaction profiles using borehole optical stratigraphy. *J. Glaciol.*, **57**(202), 289–294

Helm V, Humbert A and Miller H (2014) Elevation and elevation change of Greenland and Antarctica derived from CryoSat-2. *Cryosphere*, **8**(4), 1539–1559

Helsen MM and 7 others (2008) Elevation changes in Antarctica mainly determined by accumulation variability. *Science*, **320** (5883), 1626–1629

Herron M and Langway CC (1980) Firn densification: an empirical model. *J. Glaciol.*, **25**(93), 373–385

Hörhold MW, Kipfstuhl S, Wilhelms F, Freitag J and Frenzel A (2011) The densification of layered polar firn. *J. Geophys. Res.*, **116**(F1)

Kojima K (1967) Physics of Snow and Ice, Inst. Low Temp. Sci., Hokkaido Univ., Sapporo, Japan, vol. 1 Part 2, chap. Densification of seasonal snow cover, 929–952

Kuipers Munneke P and 9 others (2015) Elevation change of the Greenland Ice Sheet due to surface mass balance and firn processes, 1960–2014. *Cryosphere*, **9**(6), 2009–2025

Leuenberger MC, Lang C and Schwander J (1999) Delta 15N measurements as a calibration tool for the paleothermometer and gas-ice age differences: a case study for the 8200 BP event on GRIP ice. *J. Geophys. Res.: Atmos.* (1984–2012), **104**(D18), 22163–22170

Li J and Zwally HJ (2004) Modeling the density variation in the shallow firn layer. *Ann. Glaciol.*, **38**(1), 309–313

Li J and Zwally HJ (2011) Modeling of firn compaction for estimating ice-sheet mass change from observed ice-sheet elevation change. *Ann. Glaciol.*, **52**(59), 1–7

Ligtenberg SRM, Helsen MM and van den Broeke MR (2011) An improved semi-empirical model for the densification of Antarctic firn. *Cryosphere*, **5**(4), 809–819

Lüthi M and Funk M (2000) Dating ice cores from a high Alpine glacier with a flow model for cold firn. *Ann. Glaciol.*, **31**(1), 69–79

Morris EM and Wingham DJ (2011) The effect of fluctuations in surface density, accumulation and compaction on elevation change rates along the EGIG line, central Greenland. *J. Glaciol.*, **57**(203), 416–430

Morris EM and Wingham DJ (2014) Densification of polar snow: measurements, modelling and implications for altimetry. *J. Geophys. Res.: Earth Surf.*, **119**(2), 349–365

Morse DL and 7 others (1999) Accumulation rate measurements at Taylor Dome, East Antarctica: techniques and strategies for mass balance measurements in polar environments. *Geogr. Ann.: Ser. A Phys. Geogr.*, **81**(4), 683–694

Orsi AJ, Cornuelle BD and Severinghaus JP (2012) Little Ice Age cold interval in West Antarctica: evidence from borehole temperature at the West Antarctic Ice Sheet (WAIS) Divide. *Geophys. Res. Lett.*, **39**(9)

Parrenin F and 8 others (2012) On the gas-ice depth difference (Δdepth) along the EPICA Dome C ice core. *Clim. Past*, **8**(4), 1239–1255

Pitman AJ, Liang ZL, Cogley JG and Henderson-Sellers A (1992) Description of bare essentials of surface transfer for the Bureau of Meteorology Research Centre AGCM BMRC Research Report 32. Bureau of Meteorology, Melbourne, Australia

Pritchard HD and 5 others (2012) Antarctic ice-sheet loss driven by basal melting of ice shelves. *Nature*, **484**(7395), 502–505

Proksch M, Löwe H and Schneebeli M (2015) Density, specific surface area, and correlation length of snow measured by high-resolution penetrometry. *J. Geophys. Res.: Earth Surf.*, **120**(2), 346–362

Rasmussen SO and 9 others (2013) A first chronology for the North Greenland Eemian Ice Drilling (NEEM) ice core. *Clim. Past*, **9**(6), 2713–2730

Reeh N (2008) A nonsteady-state firn-densification model for the percolation zone of a glacier. *J. Geophys. Res.: Earth Surf.*, **113** (F3)

Robin G de Q (1958) Seismic shooting and related investigations *Norwegian-British-Swedish Antarctic Expedition*, 1949-52. *Scientific Results 5*.

Sandberg Sørensen L and 7 others (2011) Mass balance of the Greenland ice sheet (2003–2008) from ICESat data–the impact of interpolation, sampling and firn density. *Cryosphere*, **5**, 173–186

Schleef S and Löwe H (2013) X-ray microtomography analysis of isothermal densification of new snow under external mechanical stress. *J. Glaciol.*, **59**(214), 233–243.

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

Schwander J and 6 others (1993) The age of the air in the firn and the ice at Summit, Greenland. *J. Geophys. Res.: Atmos.*, **98**(D2), 2831–2838

Schwander J and 5 others (1997) Age scale of the air in the Summit ice: implication for glacial-interglacial temperature change. *J. Geophys. Res.*, **102**(D16), 19483–19493

Severinghaus JP and Brook EJ (1999) Abrupt climate change at the end of the last glacial period inferred from trapped air in Polar Ice. *Science*, **286**(5441), 930–934

Shepherd A and 9 others (2012) A reconciled estimate of ice-sheet mass balance. *Science*, **338**(6111), 1183–1189

Simonsen SB and 5 others (2013) Assessing a multilayered dynamic firn-compaction model for Greenland with ASIRAS radar measurements. *J. Glaciol.*, **59**(215), 545–558
Spencer MK, Alley RB and Creyts TT (2001) Preliminary firn-densification model with 38-site dataset. *J. Glaciol.*, **47**(159), 671–676

Thompson DC (1969) The coreless winter at Scott Base, Antarctica. *Q. J. R. Meteorol. Soc.*, **95**(404), 404–407

Van den Broeke MR (2008) Depth and density of the Antarctic firn layer. *Arct. Antarct. Alp. Res.*, **40**(2), 432–438

Waddington ED and Morse DL (1994) Spatial variations of local climate at Taylor Dome, Antarctica: implications for paleoclimate from ice cores. *Ann. Glaciol.*, **20**(1), 219–225

WAIS Divide Project Members (2015) Precise interpolar phasing of abrupt climate change during the last ice age. Nature, 520 (7549), 661–665

Wingham DJ, Ridout AJ, Scharroo R, Arthern RJ and Shum CK (1998) Antarctic elevation change from 1992 to 1996. *Science*, **282** (5388), 456–458

Zwally HJ and Li J (2002) Seasonal and interannual variations of firn densification and ice-sheet surface elevation at the Greenland summit. *J. Glaciol.*, **48**(161), 199–207

Zwally HJ and 7 others (2005) Mass changes of the Greenland and Antarctic ice sheets and shelves and contributions to sea-level rise: 1992–2002. *J. Glaciol.*, **51**(175), 509–527

Zwinger T, Greve R, Gagliardini O, Shiraiwa T and Lyly M (2007) A full Stokes-flow thermo-mechanical model for firn and ice applied to the Gorshkov crater glacier, Kamchatka. *Ann. Glaciol.*, **45**(1), 29–37

# APPENDIX A
## PARTICIPATING MODELS

In combination with Table 1, we provide further detail about the participating firn models. Symbols representing some variables have been altered from the referenced work to allow comparison of the constitutive relationships and empirical tuning present in the densification equations.

### A.1. HLA, HLS – (Herron and Langway, 1980)

The Herron and Langway (1980) model is a benchmark firn-densification model and is the basis for a significant proportion of later models. The model was developed to be a ‘widely applicable model of the first and second stages of polar snow densification’ (Herron and Langway, 1980). In the Herron and Langway model, the mathematical description of the rate of densification changes at critical densities and divides the firn into three densification zones. Herron and Langway (1980) derived empirical equations based on the assumption attributed to Robin (1958) that the change in density is proportional to the change in stress that results from new accumulation. That assumption leads to a linear relationship between the natural log of the density ratio $ρ/(ρ_i - ρ)$ and depth z. By assuming steady state, the model parameterizes the annual increase in overburden stress using annual accumulation rate, and the densification rate is assumed to have an Arrhenius-type temperature dependence. The model was tuned using seven depth-density profiles from Greenland and ten profiles from Antarctica. We describe three ways to formulate the Herron and Langway (1980) model: (1) analytically, (2) dynamically and (3) stress-based. We recognize that the original Herron and Langway (1980) equations have inconsistencies in the formulation of the units; we do not attempt to rectify those here except in the cases where corrections are needed to make the equations work.

#### A.1.1. HLA – analytic model

The analytic model is a steady-state model that determines the density, depth and age relationship based on the surface density $ρ_0$, accumulation rate $\dot{b}$ (m w.e. a⁻¹) and temperature T.

The density at depth z is

$$ \rho(z) = \rho_i \frac{\frac{\rho_0}{\rho_i - \rho_0} \exp\left(\frac{\rho_i}{\rho_w} k_0 z\right)}{1 + \frac{\rho_0}{\rho_i - \rho_0} \exp\left(\frac{\rho_i}{\rho_w} k_0 z\right)}, \quad (A1) \\ \rho \le 550 \text{ kg m}^{-3} $$

and

$$ \rho(z) = \rho_i \frac{\frac{\rho_{550}}{\rho_i - \rho_{550}} \exp\left[\frac{\rho_i}{\rho_w} \frac{k_1}{\dot{b}_w} (z - z_{550})\right]}{1 + \frac{\rho_{550}}{\rho_i - \rho_{550}} \exp\left[\frac{\rho_i}{\rho_w} \frac{k_1}{\dot{b}_w} (z - z_{550})\right]}, \quad (A2) \\ \rho > 550 \text{ kg m}^{-3}. $$

We have added the $1/ρ_w$ factor in the exponential terms for unit consistency. The analytic age of the firn is given by

$$ \text{Age}(\rho) = \frac{1}{k_0 \dot{b}_w^a} \ln \left[ \frac{\rho_i - \rho_0}{\rho_i - \rho} \right], \quad \rho \le 550 \text{ kg m}^{-3} \quad (A3) $$

and

$$ \text{Age}(\rho) = \frac{1}{k_1 \dot{b}_w^b} \ln \left[ \frac{\rho_i - \rho_{550}}{\rho_i - \rho} \right] + \text{Age}(550), \quad (A4) \\ \rho > 550 \text{ kg m}^{-3}. $$

The powers a and b are site-specific constants determined from the slope of the line formed by $\ln(ρ/(ρ_i - ρ))$ plotted as a function of $\dot{b}_w$. The mean and SD of a and b are respectively, 1.1 ± 0.2 and 0.5 ± 0.2. Based on that result, Herron and Langway set a = 1 and b = 0.5. The coefficients $k_0$ (units m⁻¹) and $k_1$ (units m⁻¹ a⁻¹) in the Herron and Langway model incorporate an Arrhenius rate law including temperature T in Kelvin and the gas constant R, 8.314 J mol⁻¹ K⁻¹ and are given by:

$$ k_0 = 11 \exp\left[-\frac{10160}{RT}\right] \quad (A5) $$

$$ k_1 = 575 \exp\left[-\frac{21400}{RT}\right], \quad (A6) $$

where 10160 and 21400 are the effective activation energies for the two zones (units J mol⁻¹).

#### A.1.2. Dynamic model

The dynamic Herron and Langway (1980) model assumes the rate of densification changes at critical densities. Results from the dynamical model are not included in the FirnMICE project; however, the equations are included
here since multiple future models build on these physics.

$$ \frac{D\rho}{Dt} = c_0(\rho_i - \rho) \quad \rho \le 550 \text{ kg m}^{-3} \qquad (A7) $$

$$ \frac{D\rho}{Dt} = c_1(\rho_i - \rho) \quad \rho > 550 \text{ kg m}^{-3} \qquad (A8) $$

with

$$ c_0 = \dot{b}_w^a k_0 \qquad (A9) $$

$$ c_1 = \dot{b}_w^b k_1 \qquad (A10) $$

where $a = 1$ and $b = 1/2$.

### A.1.3. HLS – stress-based model

The Herron and Langway (1980) model includes a suggested modification from reviewer Sigfus J. Johnsen, incorporating a substitution for the overburden stress $\sigma$ for the second densification zone:

$$ \frac{D\rho}{Dt} = \frac{k_1^2 (\sigma_\rho - \sigma_{550})(\rho_i - \rho)}{\rho_w g \ln[(\rho_i - 550)/(\rho_i - \rho)]} \quad (A11) $$
$$ \rho > 550 \text{ kg m}^{-3}. $$

A value of 42.6 kJ mol⁻¹ was recommended for the activation energy in Eqn (A6). We note the omission of the factor $\rho_w g$ in the denominator of Eqn (4c) of Herron and Langway (1980), and a typo in which $\rho$ was represented by $\rho$ on the left side of the equation.

## A.2. ART – Arthern and others (2010)

The Arthern and others (2010) model was originally developed for ice-sheet mass balance, but could also be applied to ice core Δage. The model used in FirnMICE is the coupled densification, heat-transfer and grain-growth model described in Appendix B of Arthern and others (2010). The dependence of the densification rate on temperature, snow load, density and grain-size is derived from a sintering theory. This theory applies to sintering by lattice diffusion creep, a grain-size dependent process. The grain-size profile is evolved dynamically by solving an equation for normal grain growth. Temperature is evolved by solving the heat equation; it influences both creep rates and grain growth.

The coupled densification, grain growth and heat-transfer system of equations (Arthern and others, 2010, Eqns (B1) and (B2)) includes densification by Nabarro-Herring creep incorporating the grain-size radius $r$,

$$ \frac{D\rho}{Dt} = k_c(\rho_i - \rho) \exp\left[-\frac{E_c}{RT}\right] \frac{\sigma}{r^2}, \qquad (A12) $$

using a grain-growth model dependent on the temperature (Gow and others, 2004),

$$ \frac{D(r^2)}{Dt} = k_g \exp\left[-\frac{E_g}{RT}\right]. \qquad (A13) $$

The model has two adjustable parameters $k_{c1}$ and $k_{c2}$ which are the coefficients for lattice-diffusion creep applied to low and high density snow respectively; values for these were

selected to give the best agreement with strain-rates measured in boreholes on the Antarctic Peninsula and Berkner Island. The coefficient $k_g$ for grain growth is $1.3 \times 10^{-7} \text{ m}^2 \text{ s}^{-1}$. The model has not been tuned to match any density-depth profiles from ice-cores. Even for simple step changes in climate forcing, the interactions between grain growth, heat conduction, firn advection, and compaction can give a complicated response.

Arthern and others (2010) also used a steady-accumulation assumption to formulate a simplified version of their model. The Nabarro-Herring creep and grain-growth physics were coupled to find rate coefficients $c_0$ and $c_1$ for the firn-densification Eqns (A7) and (A8):

$$ c_0^{(ART)} = 0.07 \dot{b} g \exp\left[-\frac{E_c}{RT} + \frac{E_g}{RT_{av}}\right], \qquad (A14) $$

and

$$ c_1^{(ART)} = 0.03 \dot{b} g \exp\left[-\frac{E_c}{RT} + \frac{E_g}{RT_{av}}\right]. \qquad (A15) $$

The activation energies for creep $E_c$ and grain-growth $E_g$ are 60 kJ mol⁻¹ and 42.4 kJ mol⁻¹. Here $\dot{b}$ has units of kg m⁻² a⁻¹ and the leading numerical coefficients have units of Pa⁻¹.

## A.3. ESS – Essery and others (2013)

The Essery and others (2013) model determines the seasonal snow-model surface boundary condition for atmospheric modeling and is therefore expected to be an outlier among the participating models, due to model tuning for mountain snow packs outside of polar regions. The densification equation,

$$ \frac{D\rho}{Dt} = \rho \frac{\sigma}{v_0} \exp\left[\frac{k_T}{T_m} - \frac{k_T}{T} - \frac{\rho}{\rho_0}\right], \qquad (A16) $$

has no limitation of compaction beyond firn densities and differs in the exponential temperature dependence. The terms include viscosity $v_0$ (kg m⁻¹ s⁻¹), reference density $\rho_0$, temperature scale $k_T$ (K); and melting point $T_m$ (K). The model is tuned to datasets, as described in Kojima (1967).

## A.4. GOU – Goujon and others (2003)

The densification scheme of Goujon and others (2003) builds upon metallurgy material science (Arnaud and others, 2000). The first stage of densification ($\rho/\rho_i < 0.6$) is based on (Alley, 1987),

$$ \frac{D(\rho/\rho_i)}{Dt} = \gamma \frac{\sigma}{(\rho/\rho_i)^2} \left(1 - \frac{5}{3}\frac{\rho}{\rho_i}\right), \qquad (A17) $$

where $\gamma$ is a function of grain-boundary viscosity and geometrical parameters, with units of inverse viscosity.

For $0.6 < \rho/\rho_i < 0.9$ (Zone 2), the densification equation is,

$$ \frac{D(\rho/\rho_i)}{Dt} = 5.3A \left(\frac{\rho/\rho_i}{D_0}\right)^{1/3} \left(\frac{a}{\pi}\right)^{1/2} \frac{\sigma^{-3}}{3}. \qquad (A18) $$
Fluidity A is given by

$$A = 7.89 \times 10^3 \exp\left[-\frac{E}{RT}\right] \quad \text{MPa}^{-3} \text{s}^{-1} \quad \text{(A19)}$$

The activation energy is $E = 60 \text{ kJ mol}^{-1}$. The densification includes an empirical coefficient $D_0$,

$$D_0 = 0.0026 \times T_s(\text{°C}) + 0.647, \quad \text{(A20)}$$

which is the relative density at the transition from Zone 1 to Zone 2. In the second zone, the overburden stress $\sigma$ due to the firn column above is all carried on limited contact areas between grains, and the effective pressure $\sigma^*$ on those grain contacts depends on the coordination number Z and on the average contact area a which is expressed as a dimensionless ratio relative to the contact area at the Zone 1–2 transition.

$$\sigma^* = \frac{4\pi\sigma}{aZ(\rho/\rho_i)}. \quad \text{(A21)}$$

For $0.9 < \rho/\rho_i > 0.95$ the densification equation is

$$\frac{D(\rho/\rho_i)}{Dt} = 2A \left[ \frac{\rho/\rho_i(1-\rho/\rho_i)}{\left(1-(1-\rho/\rho_i)^{1/3}\right)^3} \right] \left(\frac{2\sigma_{\text{eff}}}{3}\right)^3. \quad \text{(A22)}$$

The effective pressure is

$$\sigma_{\text{eff}} = \sigma + \sigma_{\text{atm}} - \sigma_b. \quad \text{(A23)}$$

For $\rho/\rho_i > 0.95$ the densification equation is:

$$\frac{D(\rho/\rho_i)}{Dt} = \frac{9}{4} A(1 - \rho/\rho_i)\sigma_{\text{eff}}. \quad \text{(A24)}$$

where $A = 1.2 \times 10^3 \exp(-E/RT) (\text{MPa}^{-1} \text{s}^{-1})$.

## A.5. LIG – Ligtenberg and others (2011)

The Ligtenberg and others (2011) model is based on the principles of Herron and Langway (1980) and is updated with the densification expressions of Arthern and others (2010). These densification equations are validated and tuned towards depth-density observations from Van den Broeke (2008). To make the model applicable on the complete ice sheet, snow-melt processes (retention, percolation, refreezing and run off) are included, although these are not used within FirnMICE. The main focus of the model is on ice-sheet mass balance and surface elevation changes.

Ligtenberg and others (2011) tuned the steady-accumulation derivation of (Arthern and others, 2010, Eqns (A14) and (A15) in this work) model, using the mean-surface temperature $\bar{T}_s$ and the same coefficients $c_0^{(ART)}$ and $c_1^{(ART)}$ for the low and high density firn:

$$c_0^{(\text{LIG})} = M_0 c_0^{(\text{ART})}, \quad \text{(A25)}$$

$$c_1^{(\text{LIG})} = M_1 c_1^{(\text{ART})}, \quad \text{(A26)}$$

The coefficients M tune the model for the low and high density firn to match depth-density profiles in Greenland and Antarctica,

$$M_0 = 1.435 - 0.151 \ln(\dot{b}) \quad \rho < 550 \text{ kg m}^{-3}, \quad \text{(A27)}$$

$$M_1 = 2.366 - 0.293 \ln(\dot{b}) \quad \rho > 550 \text{ kg m}^{-3}. \quad \text{(A28)}$$

Eqns (A25) and (A26) are substituted into Eqns (A7) and (A8) respectively to construct the firn-densification model.

## A.6. CMG – Cummings and others (2013)

The Cummings and others (2013) model was developed to simulate creation of ice lenses and the water transport in the ablation zones of ice sheets for ice-sheet mass-balance applications. The model for the internal energy $\theta = c_i T$ from the diffusion-advection equation is

$$\rho \left( \frac{\partial \theta}{\partial t} + w \frac{\partial \theta}{\partial z} \right) = \frac{\partial}{\partial z} \left( \frac{k}{c} \frac{\partial \theta}{\partial z} \right), \quad \text{(A29)}$$

where density $\rho$, thermal conductivity $k$, and heat capacity c vary in the presence of water and temperature. For the applications here, only cold firn is modeled and the thermal parameters reduce to those used by Arthern and others (2010),

$$k = 2.1 \left(\frac{\rho}{\rho_i}\right)^2, \quad c = 2009. \quad \text{(A30)}$$

The surface temperature evolves based on the equation

$$T_s = T_{\text{avg}} + 10(\cos(2\pi t) + 0.3 \cos(4\pi t)). \quad \text{(A31)}$$

The firn compaction velocity w is derived from Zwally and Li (2002, Eqn (3)),

$$w(z, t) = \int_{z_i}^{z} \frac{1}{\rho(z)} \frac{d\rho(z)}{dt} dz, \quad \text{(A32)}$$

where $z_i$ is the depth where the firn transitions into ice.

## A.7. BAR – Barnola and others (1991)

The Barnola and others (1991) model was designed to calculate the ice core $\Delta$age. For the first zone of densification, the Herron and Langway (1980) dynamic model is used. For $\rho > 550 \text{ kg m}^{-3}$, the densification model is

$$\frac{D\rho}{Dt} = \rho_i A_0 \exp\left[-\frac{Q}{RT}\right] f \sigma_{\text{eff}}^n, \quad \text{(A33)}$$

where $A_0$ is a constant, activation energy Q is 60 kJ mol$^{-1}$ and the empirical function f is written as $f_e$ for $\rho = 550 - 800 \text{ kg m}^{-3}$ and $f_s$ for $\rho > 800 \text{ kg m}^{-3}$

$$f_e = 10^{\alpha p^3 + \beta p^2 + \delta p + \gamma}, \quad \text{(A34)}$$

with $\alpha, \beta, \delta$ and $\gamma$ equal to -37.455, 99.743, -95.027 and 30.673. For $\rho > 800 \text{ kg m}^{-3}$,

$$f_s = \frac{3}{16} \frac{(1 - \rho/\rho_i)}{\left[1 - (1 - \rho/\rho_i)^{1/3}\right]^3}. \quad \text{(A35)}$$
### A.8. SIM – Simonsen and others (2013)

The Simonsen and others (2013) model was built on the empirical description of firn densification by Herron and Langway (1980) and the parameter updates of Arthern and others (2010). Further tuning of the model has been done to assimilate Greenland conditions, including the Flade Isblink and NGRIP ice cores and radar profiles of the EGIG-line (Simonsen and others, 2013). The model was developed for ice-sheet mass-balance studies and the model tuning has focused on the upper firn, but the model has been shown to be capable of giving a full description of the firn pack until the firn/ice transition. The model is aimed at giving an ice sheet wide description of changes in firn compaction over short time periods (<100 years at sub-annual resolution). The modular model framework allows for modules to be turned on/off or replaced, including a modular description of 1-D percolation of meltwater; however, this process is not used for

the FirnMICE experiments. The full model system (densification, heat transfer and percolation) is computationally efficient, which enables parameter studies by Monte Carlo inversion.

The Arthern and others (2010) coefficients $c_0^{(ART)}$ and $c_1^{(ART)}$ (Eqns (A14) and (A15) in this work) are additionally tuned by scalars $f_0$ and $f_1$

$$c_0^{(SIM)} = f_0 c_0^{(ART)}, \qquad (A36)$$

$$c_1^{(SIM)} = f_1 \frac{61.7}{b^{0.5}} \exp\left[-\frac{3800}{RT_a}\right] c_1^{(ART)}. \qquad (A37)$$

Equation (A37) includes a modification of the accumulation rate and Arrhenius temperature dependence. To formulate a densification model, $c_0^{(SIM)}$ and $c_1^{(SIM)}$ are substituted in Eqns (A7) and (A8).

MS received 23 February 2016 and accepted in revised form 12 September 2016; first published online 7 February 2017