Snowpack and firn densification in the Energy Exascale Earth System Model (E3SM) (version 1.2)

Earth’s largest island, Greenland, and the Antarctic continent are both covered by massive ice sheets. A large fraction of their surfaces consist of multi-year snow, known as firn, which has undergone a process of densification since falling from the atmosphere. Until now this firn densification has not been fully accounted for in the U.S. Department of Energy’s Energy Exascale Earth System Model (E3SM). Here, we expand the E3SM Land Model (ELM) snowpack from 1 m to up to 60 m to enable more accurate simulation of snowpack evolution. We test four densification models in a series of century-scale land 5 surface simulations forced by atmospheric re-analyses, and evaluate these parameterizations against empirical density-versusdepth data. To tailor candidate densification models for use across the ice sheets’ dry-snow zones, we optimize parameters using a regularized least squares algorithm applied to two distinct stages of densification. We find that a dynamic implementation of a semi-empirical compaction model, originally calibrated to measurements from the Antarctic peninsula, gives results more consistent with ice core measurements from the cold, dry snow zones of Greenland and Antarctica, compared to when using the 10 original ELM snow compaction physics. In its latest release, the Community Land Model (CLM) (version 5) provides updated snow compaction physics that we test in ELM, resulting in top 10 m firn densities that are in better agreement with observations than densities simulated with the semi-empirical model. Below 10 m, however, the semi-empirical model gives results that more closely match observations, while the current CLM(v5) compaction physics predict firn densities that increase too slowly with depth and are thus unable to simulate pore close off (a phenomenon of particular interest to paleoclimate studies). Because snow 15 and firn density play roles in snowpack albedo, liquid water storage, and ice sheet surface mass balance, these improvements will contribute to broader E3SM efforts to simulate the response of land ice to atmospheric forcing and the resulting impacts on global sea level.


Snowpack in Earth system and land surface models
Snow compaction in Earth system and land surface models is routinely simulated using one-dimensional parameterizations of the general forṁ 60 which relates vertical strain rates˙ to overburden pressures σ by a viscosity function η. As in the land component CLM(v4) of the Community Earth System Model (CESM) (Lawrence et al., 2011), the original ELM (version 1) further represents snow strain (i.e., fractional compaction) rates as the sum of three terms representing destructive metamorphism, overburden pressure, and melt (Anderson, 1976;Oleson et al., 2013). Densification due to destructive metamorphism (denoted by the subscript "dm"), which includes the settling of snow grains as they age, is calculated as a temperature (T ) dependent, piecewise-defined function of density (ρ), expressed by its engineering strain equivalent 1 ∆z for each vertical layer with thickness ∆z, constants c 3 = 2.78 × 10 −6 s −1 , c 4 = 0.04 K −1 , c 5 = 46 × 10 −3 m 3 kg −1 , T f = 273.15 K, and a density threshold ρ dm (100 kg m −3 ) above which the densification rate tapers off. At any given depth z, the fractional compaction rate due to overburden (load) pressure (denoted by the subscript "op") is calculated from the integrated 70 column density above z, denoted by σ, by expanding eq. (1) as 1 ∆z for each vertical layer with thickness ∆z, constants c 1 = 0.08 K −1 , c 2 = 0.023 m 3 kg −1 , and a viscosity coefficient η 0 = 9 ×10 5 kg s −1 m −2 . The snowpack model includes a one dimensional temperature model, where temperatures T (z) are calculated using an energy balance scheme in conjunction with an implicit finite difference (Crank-Nicolson) method (Jordan, 1991). The 75 radiative transfer of energy is simulated using the Snow, Ice, and Aerosol Radiative Model that includes the evolution of the ice effective grain size (Flanner and Zender, 2006;Flanner et al., 2007).
In CLM(v5), the updated fractional compaction rate due to overburden pressure is adapted from Vionnet et al. (2012), expressed as 1 ∆z threshold ρ dm (from 100 to 175 kg m −3 ), from eq. (2), above which densification due to destructive metamorphism tapers off. By updating the CLM accordingly, van Kampenhout et al. (2017) eliminated a bias in the depth where the bulk density is equal to 550 kg m −3 , which was previously much too shallow. This characteristic depth is commonly used in firn studies as a single valued proxy for a given site's full density profile and also represents the transition from the first to the second stages of densification (Herron and Langway, 1980).

Firn densification models
Empirical firn densification models typically employ analytic functions that assume a steady-state density profile. They commonly define a critical density (usually 550 kg m −3 ) above and below which represent two distinct stages of densification. Herron and Langway (1980), for example, demonstrate how their empirical firn densification model can predict observed density-depth relationships for the first two stages of densification given the mean annual temperature, annual accumulation 95 rate, and surface density. When the climate is not changing, mean annual temperatures and accumulation rates are constant, eventually leading to firn well-approximated by the steady-state condition. Assuming a steady-state, a small parcel of firn has a vertical velocity w relative to the surface, such that where A is the accumulation rate (in terms of SWE per year) and ρ is the bulk density of firn at a given depth z (Bader, 1954).

100
Neglecting wind shear, the one-dimensional (kinematic) densification rate can then be derived from the material derivative with ∂ρ/∂t = 0 in a steady-state. Substituting eq. (5) into eq. (6) then gives the advective densification rate, which is closely related to the volumetric strain rate˙ via 105 This steady-state approach is useful for deriving realistic density profiles and vertical strain rates in dry snow zones, but does not provide a dynamic representation of physical processes simulated in modern ESMs. ingly, the densification rate ∂ρ/∂t can be expressed in terms of temperature T , bulk density ρ, overburden pressure P , and an effective snow grain radius r e , such that with activation energy E c (60 kJ mol −1 ), universal gas constant R (8.31 J K −1 mol −1 ), ice density ρ i (917 kg m −3 ), and k c = 9.2 × 10 −9 kg −1 m 3 s for ρ ≤ 550 kg m −3 or k c = 3.7 × 10 −9 kg −1 m 3 s for ρ > 550 kg m −3 . Adjustment of the rate 115 4 https://doi.org/10.5194/gmd-2020-247 Preprint. Discussion started: 13 August 2020 c Author(s) 2020. CC BY 4.0 License. coefficient k c for ρ ≤ 550 kg m −3 is necessary to capture greater densification rates during the first stage of densification due to grain-boundary sliding (Alley, 1987). Because their study was confined to sites on or near the relatively warm Antarctic Peninsula, however, the model should not be applied in ESMs to represent the cold interior regions of Greenland and Antarctica without further calibration.

Firn model development
In this section, we describe model development resulting in a set of ELM configurations used to simulate the accumulation of snow and firn densification on land ice. Our modifications in E3SM began with the snow model in ELM(v1), which was inherited from the CLM(v4.5) used in CESM(v1). After adopting the 12 layer grid introduced into CLM(v5) by van Kampenhout et al. (2017), we expanded and modified their improved vertical grid to improve the spatial resolution at snow depths of 10 to 125 60 m. This improved spatial resolution is necessary to better simulate relevant processes at depths (and overburden pressures) more typical of firn.

130
upper-most nine layers vary in their minimum and maximum allotted thicknesses and conform to the top nine layers used in the 12 layer grid described by van Kampenhout et al. (2017). Expanding their 12 layer grid to a new 16 layer grid, we modified layers 10 and 11 (counting from the top) and appended onto the bottom an additional 4 so that layers 10-15 have an equal (7.67 m) spacing. When the snowpack becomes deep enough to fill these 15 layers, a semi-infinite (bottom-most) layer is created.
This new layering scheme yields a vertical resolution of roughly 8 m for firn depths down to 60 m while maintaining its semi-135 infinite capability. It also maintains the variable spacing near the snowpack surface that is needed to resolve high temperature gradients and simulate radiative transfer.
Equipped with a 16 layer vertical snowpack grid, the simulation of firn in the ELM is enabled by removing the upper bound on the mass of snow (1 m SWE). This allows for the growth of semi-infinite firn columns as snow accumulates. With these changes in place, a baseline ELM configuration (hereafter referred to as "A'76," short for Anderson (1976)) is used to evaluate 140 the original snowpack and firn densification scheme. For dry snow-zones (columns without any melting), vertically dependent, fractional compaction rates in this configuration are calculated as the sum of eq. (2) and eq. (3).

Improved ice sheet surface densities
In the A'76 configuration (Sect. 3.1.1), snowfall that accumulates during a time step is assigned a density independent of the wind speed. This temperature-only dependence is one of the model deficiencies attributed to ice sheet surface densities that 145 are too low. To address this and other missing physics, two key improvements from van Kampenhout et al. (2017) were implemented into ELM. These include a wind speed dependence in the fresh snow density parameterization and the addition of a densification term that incorporates compaction (up to 350 kg m −3 ) due to drifting snow. For the purpose of comparing experimental results to the current CLM(v5) configuration, we updated the ELM accordingly, including the overburden compaction parameterization from Vionnet et al. (2012), shown in eq. (4). These updates, combined with CLM's current 12 layer 150 grid as well as a destructive metamorphism (dm) density threshold ρ dm = 175 kg m −3 in eq. (2), are fully described by van Kampenhout et al. (2017), which we refer to here as "vK'17." 3.1.3 Compaction due to snowpack load pressure Next, compaction due to the overburden pressure from snow loading was updated according to Arthern et al. (2010) by replacing eq. (3) with eq. (8). Here, the overburden pressure (or "load") is represented in the column as a function of the overlying, 155 integrated column density σ, so that where g is the acceleration due to gravity (9.81 m s −2 ) and the factor ρ i /ρ approximates the increased grain-load stress due to pore space within the snowpack (Cuffey and Paterson, 2010). From the "vK'17" configuration (Sect. 3.1.2), another configuration, hereafter referred to as "A'10" (short for Arthern et al. (2010)), is defined using this updated (semi-empirical) 160 overburden pressure model, which also includes the major improvements from van Kampenhout et al. (2017), described in Sect. 3.1.2. In this configuration, we keep the density threshold ρ dm in eq. (2) at 175 kg m −3 and apply our 16 layer grid ( Fig.   1).
Later (but not included in A'10), we added a fresh snow compaction term, calibrated for dendritic snow, from , where Because ELM does not specify snow grain shape or type, we added this term only for snow having a low-enough snow grain size (i.e., where its layer-dependent optical sphere equivalent radius is less than roughly 80 µm) to be considered dendritic.

Statistical modeling
Implementation of candidate firn densification models into ELM provides a direct method of evaluation against observations.

170
Although such comparisons provide complete tests, they require century scale simulations that, because of their computational expense, hinder the outcome of sensitivity studies needed for development. To circumvent this problem, statistical modeling was used to estimate empirical compaction rates, which are then used to calibrate the more physically-based models discussed above (A'76, vK'17, and A'10).
To approximate a statistical model of the dry snow zones across Antartica and Greenland, Monte Carlo experiments were 175 conducted to generate thousands of plausible firn density-versus-depth profiles. First, random mean-annual temperatures (T ) below -25 • C were drawn from the left tail of a Gaussian distribution (with an estimated global mean µ = 14.9 • C and standard deviation = 16 • C) selected to give a distribution of temperatures crudely representative of Earth's land surface. Next, accumulation rates (A) were drawn at random from a lognormal distribution selected to give values representative of warm (T > -51 • C) or cold (T ≤ -51 • C) dry snow zones, with 0.07 < A < 0.4 or A < 0.07 m SWE yr −1 , respectively. Valid mean 180 annual temperature and accumulation rate pairs were then combined with a surface density (ranging from 300 to 380 kg m −3 ) drawn at random from a uniform distribution. These values were inserted in the empirical model of Herron and Langway (1980) from which we derive empirical steady-state strain rates (with a vertical resolution of 10 cm) from eq. (7).
From our estimated empirical strain rate-versus-depth data, we optimized the previously described densification model coefficients (from A'76, vK'17, and A'10) by applying a regularized least squares algorithm for two stages of densification 185 (above and below ρ = 550 kg m −3 ). Normalized covariances (i.e., correlation R) between empirical strain rates and those from optimized firn densification parameters were inferred to evaluate the skill of each set of parameters.

Firn density simulations
Although statistical modeling enables a basic evaluation of experimental densification models, it does so using empirically derived, steady-state density profiles, rather than those derived from the densification models themselves. To calculate density 190 profiles that are truly stable for a given densification model, a numerical scheme is needed to integrate accumulating firn to a steady state. The Common Infrastructure for Modeling the Earth (CIME) provides the software necessary to fully simulate the accumulation of snow and densification of firn on land ice with ELM. In ELM stand-alone mode, the surface boundary condition is provided by atmospheric re-analysis data that include diurnally (3 hour) varying precipitation, solar radiation, temperature, and wind speed from the Climate Research Unit and the National Center for Environmental Prediction  CEP) (Viovy, 2018). The following historical climate simulations were conducted using coarse-resolution, stand alone ELM (an "I-compset" at ne11 resolution).

Equilibrium climate
To spin-up the model, we initialized (starting on January 1st) snowpack conditions using a thin (50 mm SWE) cover on ELM's glacier columns (including those representing the Greenland and Antarctic Ice sheets) as well as all land columns located north This initial, cold start, condition removes any prior assumption about the natural firn density profiles, resulting in firn densities representative of densification processes simulated entirely within the model.

Twentieth century climate
Twentieth century simulations are initialized using conditions at the end of 260-year equilibrium climate simulations (i.e., as 205 "restart runs"). One-hundred years of firn densification over land ice were then simulated using atmospheric forcing data from

Results and Discussion
To evaluate firn density simulations, we test a series of compaction parameterizations in ELM with a semi-infinite vertical snowpack grid ( simulation results (Fig 2.). In these initial simulations, densities are highly variable in the top 3 m. Densities increase from less than 300 kg m −3 at the surface to as much as 600 kg m −3 within the first 5 m for relatively warm dry snow zones where the mean annual temperature is within a couple degrees of -25 • C. Compared to the widely used steady-state empirical model from Herron and Langway (1980), ELM-simulated densities become too large below the top two meters. High densities persist 230 down to a depth of 20 m for all temperature regimes, a model deficiency also identified and described in van Kampenhout et al. (2017).
To improve firn density simulations, we implemented key changes from van Kampenhout et al. (2017)  To better understand what drives rapid densification near the surface, we approximate vertical strain rates (negative for deformation) and firn age from Sorge's Law (Bader, 1954). Assuming a steady-state, we compare densities and fractional 245 compaction rates as a function of firn age (Fig. 4). These data show rapid densification within the first ten years. Near the surface, fractional compaction rates, which are equivalent to negative strain rates, can become as large as 10 −6 s −1 (9 % per day), after which they decrease by three orders of magnitude (to below 3 % per year) after just 20 years. Densification tapersoff at lower densities (around 450 kg m −3 ) for colder climates, a temperature-dependent effect enhanced with the model from Arthern et al. (2010).

250
To test our statistical calibration method, we repeated ELM experiments, but replaced the pressure compaction equation with an optimized version described by Vionnet et al. (2012). We also set the density threshold ρ dm in eq. (2) to 150 kg m −3 , reverting back to its original value given by Anderson (1976). This optimized configuration produces compaction rates at least moderately correlated (R 2 = 0.4) with empirical modeling, though only for the second stage of densification (ρ >550 kg m −3 ).
For reference, we compare steady-state density profiles from this additional ELM experiment, referred to as "vK'17+," to the 255 A'10 ELM experiment ( Fig. 3 and Fig. 4). Note also (from Table 1) that this vK'17+ configuration includes the dendritic snow compaction term from eq. (10) , though this addition did not have a noticeable effect on the density profile below 1 m. Like in the original CLM (i.e., our vK'17 run), results from vK'17+ show slightly improved top 10 m near-surface densities for warm regions, but under-estimate densification as depths increase.
Although the semi-empirical model improves the density profile, its implementation into ELM, does not go far enough to 260 represent the main drivers of snow compaction that occur at relatively low densities (< 500 kg m −3 ). With simple represen-tations of compaction due to destructive metamorphism and other relevant processes (e.g., grain-boundary sliding), modeling challenges near the surface arise where densities range from 300 to 500 kg m −3 and vary due to sub-grid scale snow microstructure properties not yet accounted for in firn models (Lundin et al., 2017). Despite widespread over-densification of near-surface firn, colder regions develop steady-state density profiles that vary too weakly with depth. Their effective strain 265 rates, however, are similar to those predicted with empirical modeling for low accumulation rates (i.e., those less than 0.15 m SWE yr −1 ). These results indicate that the 1901-1920 atmospheric forcing data contain accumulation rates (from snowfall) that are too high in colder regions of the ice sheet. Unrealistically high accumulation rates, and the resulting faster downward advection of firn relative to the surface, explains why kinematic densities (density-versus-depth relationships) appear too low while undergoing realistic compaction rates. 3. adding a constant compaction term equal to (-)1.18 × 10 −10 s −1 .
As expected, these changes, representing an a priori calibration, better capture densification below 10 m, slightly improving the steady-state comparison between our simulations and Herron and Langway (1980). Because our a priori calibration uses a mean annual temperature domain that includes the coldest regions of the Antarctic ice sheet, we might expect further improved 280 density profiles for warm regions by restricting our calibration to mean annual temperatures greater than -34 • C. Also worth mentioning is that similar results can be produced using the original pressure compaction term from Anderson (1976) when increasing the viscosity coefficient η 0 by factor of 50. This is because the viscosity η has essentially the same functional form as the one from Vionnet et al. (2012), but with slightly different coefficients (for which we adjust via a least squares regression).

Regression modeling 285
To improve model accuracy, we use Monte Carlo experiments to calibrate densification model coefficients. Covariance matrices, indicating model specific a priori variances and inter-model covariances (and correlations), can be generated for multiple experimental parameterizations without conducting full ELM simulations (Table 2). Diagonal entries of the multi-model covariance matrix give model-specific variances (SD 2 ), while off-diagonal entries give covariances that are normalized to calculate empirical correlations (R). Surprisingly, a slightly modified destructive metamorphism expression from eq. (2) (Anderson, 290 1976), by itself, results in compaction rates that are more highly correlated with the empirical model (Herron and Langway, 1980) than that from unmodified, combined parameterizations that include both destructive metamorphism and compaction due to overburden pressure. We found this higher correlation by changing the density tapering constant (from 46 to 12.5 cm 3 g −1 ) in eq. (2). This modified destructive metamorphism expression is most highly correlated with empirical compaction rates for densities greater than 600 kg m −3 , but does not co-vary for densities less than 450 kg m −3 . Given that its original purpose 295 was to model low density (less than 250 kg m −3 ) compaction rates undergoing temperature dependent metamorphism, its moderately high correlation with empirical strain rates at relatively high densities does not justify a sufficiently quantitative densification mechanism.
Using the semi-empirical, pressure-driven model developed by Arthern et al. (2010), we are able to capture more variability (R 2 = 0.67) in strain rates in the colder regions across the interior of the Greenland and Antarctic Ice Sheets (Fig 5). This result 300 is encouraging, as this semi-empirical model is constructed with data from the Antarctic Peninsula and Berkner Island.
After applying the method of least squares for two stages of densification, the optimized models have smaller variances than the empirical model. A lower model variance occurs when it does not covary with the empirical model. This effectively reduces a model's prediction risk if it does not also result in an increased bias. Although further statistical inference can assess this bias-variance tradeoff, we can directly evaluate model performance with ELM simulations.

305
Across the ice sheets' dry snow zone temperature domain, roughly -60 to -25 • C, these results indicate negative correlations between overburden pressure and empirical strain rates. Consequently, regression models, which we use to estimate model coefficients, minimize pressure compaction terms to calibrate for decreasing strain rates with increasing overburden pressure.
Without including either a destructive metamorphism term or grain-size dependence in the model, the regression coefficient for a pressure term can even be negative. While firn densification is usually considered analogous to pressure sintering, where 310 porous materials deform under stress, Meussen et al. (1999) determined that the density itself is the most crucial variable in approximating the densification rate. These findings, as well as other firn model inter-comparison studies (Lundin et al., 2017;Verjans et al., 2019), emphasize the need for further development of constitutive, density-dependent stress versus strain relationships which can also account for microphysical, grain-scale features.
While the semi-empirical firn densification model is tailored for Antarctic Peninsula sites, we hypothesize its implemen-315 tation in ELM will improve density profiles on a site-by-site basis after optimizing its coefficients from temperatures and accumulation rates that are representative of its target domain, including the ice sheet interior. Even without optimization, we find improvements in the simulated density profile in equilibrium simulations, though our analysis with ELM thus far is limited to a generalized comparison with a broad (climate) perspective rather than to a more site-specific comparison against direct observations. Although these developments are ongoing, we might expect improved near-surface densities, based on our 320 statistical computing, with the following changes to the A'10 experiment: 1. To further improve surface densities, add the expression for dendritic snow compaction from eq. (10) , as we did for vK'17+; 2. Set the destructive metamorphism coefficient c 3 = 5.8 × 10 −7 s −1 , from eq. (2); 3. Set the creep coefficient k c = 1.4 × 10 −9 kg −1 m 3 s when ρ ≤ 550 kg m −3 and k c = 1.2 × 10 −9 kg −1 m 3 s when ρ > 550 325 kg m −3 , from eq. (8); 4. Add a constant compaction term, as we did for vK'17+, but equal to (-)2.02 × 10 −10 s −1 when ρ ≤ 550 kg m −3 and equal to (-)2.7 × 10 −11 s −1 when ρ > 550 kg m −3 .

End of 20th century Greenland and Antarctic Ice Sheet case studies
To evaluate the ELM snow and firn densification mechanics against recent observations, we continued the steady-state experi- to recent firn core measurements (Fig. 6). We average data in this manner (i.e., over the whole ice-sheet for Greenland versus only 4 grid cells nearest Siple Dome) to most closely represent model results with respect to observations. In general, the Our ELM simulations are currently limited by their spatial resolution. Because calculations across spatially varying grid cells are independent, i.e., they have no ability to influence one another, we seek a balance between computational expense versus 340 adequate spatial sampling. The coarse ne11 global resolution enables short wall-clock time, multi-century duration simulations while still providing multiple grid cells on both ice sheets that develop deep (> 60 m) firn columns. These coarse-resolution simulations do, however, present challenges in evaluating regions near coasts as atmospheric forcing data are averaged over large grid-cells that include both land and ocean surface types. Near Siple Dome, Antarctica, for example, this large gridcell remapping lead to a cold bias, resulting in too-slow densification. Therefore, we adjusted our Siple Dome comparisons 345 to include grid-cells away from the coast that better represent atmospheric conditions and result in a more realistic density simulation.
Encouragingly, our simulation results compare well with firn density measurements and indicate an improved capability in the ELM. Moving forward, we have yet to address whether the new ELM firn model satisfactorily represents ice sheet SMB.
To address this question, we should focus on the near surface layers, as they contain the primary SMB components. To capture 350 melt water percolation, retention, and refreezing, it could be necessary to model the upper most 20 m, but accurate modeling of firn densities all the way down to pore close-off might be unnecessary. As a result, the optimized version (vK'17+) is likely the better choice for implementation into the next major release of the E3SM. Although it shows a low density bias at depth, this bias might actually be preferable to a high bias if we aim to partition run-off regimes -which are characterized by thick, dense, near-surface ice lenses -from percolation regimes.

355
On the other hand, we have yet to test in ELM an optimized version of the semi-empirical model, which could be progressively tailored to accurately reproduce near surface densities. Furthermore, instead of using Monte Carlo integration to determine a theoretical mean annual temperature-accumulation probability distribution, an improved method would use data from reanalyses or regional climate modeling as inputs to the empirical model from Herron and Langway (1980). Such an approach could be tailored to dry-snow zones of the GrIS, for example, by restricting data to regions where the surface tem-360 perature never exceeds 0 • C.
The next major step, however, is a full evaluation of the model's firn densification mechanics that include the effects of liquid water. Such an evaluation remains challenging, as observations in these regions are limited. Such limitations call for more observational studies where the ice sheet is beginning to melt. In the near future, this could be the entire GrIS, so accurate modeling of these regimes will be crucial.

Conclusions
We considered different compaction models to improve firn simulation in ELM, and compared simulated density profiles with density measurements of firn cores. Expanding this evaluation to include regions across the dry snow zones of the Greenland and Antarctic Ice Sheets, we examined the variation of density with depth in equilibrium climate simulations and compared results with steady-state empirical models (Herron and Langway, 1980). Our analysis, similar to that by van Kampenhout Ultimately, this study seeks to enable better predictions of SLR as a direct result of surface melt and mass loss from the GrIS. Accurate partitioning of the SMB terms, however, remains a challenge in ESMs. Although regional climate models 380 currently enable a higher spatial resolution, limitations in our understanding of supra-glacial hydrology still remain a source of uncertainty in determining Greenland's future contribution to SLR. Future developments will seek to improve the capability of simulating supra-glacial hydrology in the E3SM, including better quantifying surface melt, percolation, refreezing, and the build-up of perennial aquifers.
Code and data availability. The E3SM (version 1) source code is maintained and available at https://github.com/E3SM-Project. Source 385 code modifications and model development that will reproduce these results are available from Edwards et al. (2020). ELM simulation data, associated python analysis scripts, and the offline statistical firn model referenced in this manuscript are archived and publicly accessible from Schneider (2020). The CLM (version 5) source code is maintained and available at https://github.com/ESCOMP/CTSM.
Author contributions. AS collected published data, modified E3SM source code, conducted ELM simulations and formal analyses, developed the statistical model, curated data, and prepared the figures, tables, and paper. CZ and SP supervised and co-led the project conceptualization, funding acquisition, and administration and helped write the paper.
Competing interests. The authors are not aware of any conflicts of interest regarding the publication of this manuscript.     . End of 20th century firn density profiles simulated with ELM compared to firn core measurements from Greenland (Mosley- Thompson et al., 2001) and Siple Dome (Lamorey, 2003). Line graphs represent ELM simulations with two different densification models: a semi-empirical formulation (Arthern et al., 2010), in black; and a calibrated, 16-layer version of the current CESM Land Model configuration, in gray. ELM simulation results from Greenland represent regional means calculated from 8 model grid cells with mean annual temperatures less than or equal to -25 • C. ELM simulation results from Antartica represent regional means calculated from 4 grid-cells located near Siple Dome.