Research article 14 Jan 2022
Research article  14 Jan 2022
On the contribution of grain boundary sliding type creep to firn densification – an assessment using an optimization approach
 ^{1}Institute of Applied Mechanics, Technische Universität Kaiserslautern, Kaiserslautern, Germany
 ^{2}Division of Solid Mechanics, Technische Universität Darmstadt, Darmstadt, Germany
 ^{3}AlfredWegenerInstitut HelmholtzZentrum für Polar und Meeresforschung, Bremerhaven, Germany
 ^{4}Faculty of Geosciences, University of Bremen, Bremen, Germany
 ^{1}Institute of Applied Mechanics, Technische Universität Kaiserslautern, Kaiserslautern, Germany
 ^{2}Division of Solid Mechanics, Technische Universität Darmstadt, Darmstadt, Germany
 ^{3}AlfredWegenerInstitut HelmholtzZentrum für Polar und Meeresforschung, Bremerhaven, Germany
 ^{4}Faculty of Geosciences, University of Bremen, Bremen, Germany
Correspondence: Timm Schultz (tschultz@rhrk.unikl.de)
Hide author detailsCorrespondence: Timm Schultz (tschultz@rhrk.unikl.de)
Simulation approaches to firn densification often rely on the assumption that grain boundary sliding is the leading process driving the first stage of densification. Alley (1987) first developed a processbased material model of firn that describes this process. However, often socalled semiempirical models are favored over the physical description of grain boundary sliding owing to their simplicity and the uncertainties regarding model parameters. In this study, we assessed the applicability of the grain boundary sliding model of Alley (1987) to firn using a numeric firn densification model and an optimization approach, for which we formulated variants of the constitutive relation of Alley (1987). An efficient model implementation based on an updated Lagrangian numerical scheme enabled us to perform a large number of simulations to test different model parameters and identify the simulation results that best reproduced 159 firn density profiles from Greenland and Antarctica. For most of the investigated locations, the simulated and measured firn density profiles were in good agreement. This result implies that the constitutive relation of Alley (1987) characterizes the first stage of firn densification well when suitable model parameters are used. An analysis of the parameters that result in the best agreement revealed a dependence on the mean surface mass balance. This finding may indicate that the load is insufficiently described, as the lateral components of the stress tensor are usually neglected in onedimensional models of the firn column.
 Article
(3477 KB) 
Supplement
(155 KB)  BibTeX
 EndNote
Firn densification models fall into two basic categories. Models in the first category, which includes most existing models, follow the socalled semiempirical approach of Herron and Langway (1980), which itself is based on Sorge's law (Bader, 1954) and the Robin hypothesis (Robin, 1958). Examples are the models of Arthern et al. (2010), Ligtenberg et al. (2011), and Simonsen et al. (2013). The empirical parameters of these models are typically adjusted on the basis of certain datasets of depth density profiles. In the second category of firn densification models, an attempt is made to quantify the physical processes related to firn densification. These processes include various types of creep and diffusion. Micromechanical models are used for smallscale investigations (Johnson and Hopkins, 2005; Theile et al., 2011; Fourtenau et al., 2020), whereas models based on continuum mechanics can be used for largescale simulations. Examples of the latter are the models of Arthern and Wingham (1998), Arnaud et al. (2000), and Goujon et al. (2003).
Alley (1987) first applied the theory of grain boundary sliding adopted from Raj and Ashby (1971) to firn densification at densities below the critical density of ρ_{c}=550 kg m^{−3}. The description of this process by Alley (1987) was subsequently used in other firn densification models (Arthern and Wingham, 1998; Arnaud et al., 2000; Goujon et al., 2003; Bréant et al., 2017). However, the assumption that grain boundary sliding is the dominant process in firn densification at densities below ρ_{c}=550 kg m^{−3} has been questioned numerous times (Ignat and Frost, 1987; Roscoat et al., 2010). For example, Theile et al. (2011) conducted experiments on a small number of snow samples and suggested that densification is more likely driven by processes within the grain rather than by the intergranular process of grain boundary sliding.
In this study, we aim to evaluate (i) whether the description of grain boundary sliding given by Alley (1987) is suitable for the simulation of firn densification at low density, (ii) how a modification of the constitutive relation introduced by Bréant et al. (2017) affects simulation results, (iii) whether hidden or additional dependencies on climatic or other conditions can be identified in the constitutive relation, and (iv) how the constitutive relation of Alley (1987) might be improved. Note that our study aims to assess the constitutive relation for grain boundary sliding proposed by Alley (1987). An evaluation that clarifies whether grain boundary sliding is the dominant process driving firn densification below the critical density of ρ_{c}=550 kg m^{−3} must be conducted using other methods. Experimental attempts to do so have been made by, for example, Kinosita (1967), Ignat and Frost (1987), and Theile et al. (2011). In contrast to these experimental investigations, a datadriven model approach is used in our study. Since the original study of Alley (1987) was published, the number of available data have increased greatly. The data include a large number of firn profiles and forcing data, and they allow us to simulate firn profiles at very high quality using additional modeling techniques.
To test the description of grain boundary sliding given by Alley (1987), we used a numerical model to simulate the evolution of a onedimensional firn column with time. The model uses variants of the constitutive relation of Alley (1987), all of which combine several model parameters in a single factor. We force the model with data provided by the regional climate model RACMO2.3 (Van Wessem et al., 2014; Noël et al., 2015), which represents the climate of recent decades at 159 locations where firn density measurements were made. These firn measurements are available through the Surface Mass Balance and Snow Depth on Sea Ice Working Group (SUMup) snow density subdataset (Koenig and Montgomery, 2019). By varying the factor used in the variants of the constitutive equation, we produced a large number of simulation results, which were compared to the corresponding density measurements. The quality of the factors used in the simulations was evaluated in terms of the deviation of the computed density profiles from the measured profiles. The factor values that yielded the best results reveal possible improvements in the description of grain boundary sliding for firn densification at low density. In the following sections, the constitutive equation for grain boundary sliding given by Alley (1987), the optimization scheme, and the density and forcing data used are described. A detailed description of the model is presented in the Appendix (Sect. A).
2.1 Constitutive relation for grain boundary sliding
We explain briefly the components and characteristics of the constitutive law of Alley (1987) describing the process of grain boundary sliding.
The factor of $\mathrm{2}/\mathrm{15}$ results from the geometric deviation pointed out by Alley (1987). Another geometric parameter, δ_{b}, describes the width of the grain boundary.
The following part of the equation describes the reciprocal bond or boundary viscosity (Raj and Ashby, 1971). The optimization approach of Alley (1987) was intended to identify the optimal values of the boundary viscosity. Alley (1987) compared the results of this optimization to the description of the boundary viscosity given by Raj and Ashby (1971), which includes the volume of the H_{2}O molecule Ω, the Boltzmann constant k_{b}, the temperature T, and the amplitude of grain boundary obstructions h. The latter is a measure of the roughness of the grain boundary. D_{BD} is an Arrhenius factor describing the rate of boundary diffusion. Values of the typical activation energy for this process, Q_{BD}, and the corresponding prefactor, A_{BD}, can be found in the literature (e.g., Maeno and Ebinuma, 1983) and are further discussed in Sect. 2.2. R is the universal gas constant.
The strain rate resulting from grain boundary sliding also depends on the grain radius r. The ratio of the grain radius to the neck radius μ was introduced by Arthern and Wingham (1998) and is assumed to be constant. Various methods can be used to determine the size of grains in crystalline materials (e.g., Gow, 1969). The model of Alley (1987) was developed assuming perfectly spherical grains. Although this assumption is not true of firn, it provides a reasonable basis for modeling. Therefore, throughout this paper, the grain radius r represents the radius of theoretical spherical grains.
The next factor in Eq. (1) describes the dependence on the inverse density relative to the ice density cubed. The factor of $\mathrm{1}(\mathrm{5}\mathit{\rho}/\mathrm{3}{\mathit{\rho}}_{\mathrm{ice}})$ causes the strain rate due to grain boundary sliding to decrease with increasing density until it vanishes at the critical density of ρ_{c}=550 kg m^{−3}. When the critical density is reached, close random packing is established (Anderson and Benson, 1963), and grains can no longer slide against each other; thus, the process of grain boundary sliding ends. Other deformation processes, in particular dislocation creep (Maeno and Ebinuma, 1983), result in further densification with increasing stress.
Alley (1987) suggested that additional processes contribute to densification below the critical density. It is feasible that the effect of grain boundary sliding on the strain rate decreases, whereas that of other processes increases. The studies of Arthern and Wingham (1998) and Bréant et al. (2017) use this description, in which only grain boundary sliding drives densification in the first stage of firn densification. In the study of Bréant et al. (2017), the constitutive relation of Alley (1987) is changed such that the transition into the second stage of densification is modified. We evaluate this modification in this work.
Finally, the stress in the vertical direction, σ_{zz}, resulting from the overburden firn drives grain boundary sliding. Whereas Alley (1987) used the product of the accumulation rate, acceleration due to gravity, and time since the deposition of a specific firn sample to describe the overburden stress, we use a more general form here (see Sect. A2, Eq. A5). The other physical properties affecting the process are density ρ, temperature T, and grain radius r.
2.2 Optimization
To test the material model developed by Alley (1987), we formulated variants of Eq. (1) and compared the model results to density measurements of various firn cores. These variants of the constitutive Eqs. (2) to (5) preserve its general form but group several material parameters into a single factor. As a result, the simulation result does not depend on those parameters but on the single factor. The factor is then varied to find an optimal simulation result that best reproduces the measured firn profile. The factor that yields the optimal simulation result depends on the measured firn density profile and the corresponding climate conditions. It is therefore sitespecific. This feature makes it possible to assess whether the description of grain boundary sliding given by Alley (1987) can be used to reproduce measured firn profiles, assuming an optimal set of parameters. It further allows us to analyze the sitespecific factors yielding the best simulation results for possible hidden dependencies.
Arnaud et al. (2000), Goujon et al. (2003), and Bréant et al. (2017) also incorporated the material parameters of the model of Alley (1987) into a single parameter. Bréant et al. (2017) also modified the factor of $\mathrm{5}/\mathrm{3}$ to change the density at which the deformation due to grain boundary sliding becomes zero. In the following, four variants, indicated by the subscripts $(\cdot {)}_{{\mathrm{v}}_{\mathrm{1}}}$ to $(\cdot {)}_{{\mathrm{v}}_{\mathrm{4}}}$, are presented:
Variant 1 (Eq. 2) and Variant 2 (Eq. 3) of the constitutive equation combine all material constants using the factors ${C}_{{\mathrm{v}}_{\mathrm{1}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{2}}}$, respectively. The Arrhenius factor for boundary diffusion, D_{BD} (see Eq. 1), is retained in these variants. Following Maeno and Ebinuma (1983), we use a value of Q_{BD}=44.1 kJ mol^{−1} for the boundary diffusion activation energy. This variable was defined by Maeno and Ebinuma (1983) as twothirds of the activation energy for lattice diffusion measured by Itagaki (1964). The corresponding prefactor is ${A}_{\mathrm{BD}}=\mathrm{3.0}\times {\mathrm{10}}^{\mathrm{2}}$ m^{2} s^{−1}. Alley (1987) assumed a similar value for the boundary diffusion activation energy.
In addition to the temperature T, the vertical strain rate ${\dot{\mathit{\epsilon}}}_{zz}$ depends only on the firn density ρ, grain radius r, and stress in the vertical direction σ_{zz}. Variant 2 differs from Variant 1 in that it uses the modification introduced by Bréant et al. (2017). This modification causes a theoretical end to the process of grain boundary sliding at a density of ${\mathit{\rho}}_{\mathrm{c}}^{*}=\mathrm{596}$ kg m^{−3}. It was introduced to obtain a better transition to the second stage of firn densification. The strain rate due to grain boundary sliding is therefore higher at the critical density when this modification of Bréant et al. (2017) is used.
To determine the effect of the Arrhenius factor, it is disregarded in variants 3 and 4, as shown in Eqs. (4) and (5):
In Variant 4, the modification of Bréant et al. (2017) is used, whereas Variant 3 uses the original formulation of Alley (1987). The goal of the optimization is to find optimal values of the factors C_{v} for every variant of the constitutive relation (Eqs. 2 to 5) to produce simulated density profiles that best reproduce the measured profiles. As an example, we explain the optimization process for one selected firn core in more detail. The upper part of ice core ngt03C93.2 (Wilhelms, 2000) is shown in Fig. 1a.
Every simulation begins with a spinup period in which constant values are used for forcing. The model is forced with prescribed values of temperature, accumulation rate, firn density, and grain radius at the surface. We check for steadystate conditions by comparing the change in density between time steps. If the maximum density change is smaller than Δρ_{max}<0.1 kg m^{−3}, we assume that the steady state is reached. In this case, a transient run using varying forcing data follows. We use a constant value of 48 time steps per year for spinup and transient simulation runs (see Sect. A7). The forcing at the location of ngt03C93.2 is shown in Fig. 1c. The resulting firn profile is then compared to the measured profile. We used the root mean square deviation (RMSD) between the measured and modeled density for comparison, as it is simple and easy to compute. To calculate the deviation, the simulated density values are interpolated linearly to the measurement locations along the profile. To ensure the quality of the results, we limited the calculation of the deviation to the domain defined by the location of the uppermost available measurement point and the oldest simulated horizon within the firn profile affected by forcing. For ngt03C93.2, this horizon is the 1958 surface at a depth of approximately 11 m below the surface, indicated by dashed horizontal lines in Fig. 1a. Only results located above the simulated 1958 surface are used to calculate the deviation. The reason is that the forcing data from RACMO2.3 (Van Wessem et al., 2014; Noël et al., 2015) begin in 1958 for Greenland. Consequently, the results are not affected by the spinup period. For firn profiles retrieved in Antarctica, climate forcing from RACMO2.3 begins in 1979. Thus, only those results located above the simulated horizon of 1979 are considered for comparison with the Antarctic firn profiles. The examination of other firn cores revealed that the surface of the oldest available forcing may be located at greater depth when the density of ρ_{c}=550 kg m^{−3} has already been reached. In those cases, the computation of the RMSD was limited to the domain showing density values below 540 kg m^{−3}. We decided to use a density threshold below the critical density because of the asymptotic nature of the resulting density profiles obtained using variants 1 and 3 of the constitutive equation (Eqs. 2 and 4). The value of 540 kg m^{−3} was chosen to ensure that the results obtained using the variants of the constitutive relation are comparable, whereas unique values of the factor C_{v} were quickly determined throughout the optimization.
As the implementation of our model is efficient and the approach is simple and reliable, we decided to determine the best factor C_{v} for the four variants of the constitutive equation by simply testing 250 values within certain ranges. These ranges are shown in Eqs. (6) and (7), which include and exclude the Arrhenius factor, respectively. Optimal factors can be found within these ranges for every analyzed firn profile. To ensure that this is the case, all simulations were performed multiple times using different ranges of the factors.
Figure 1b shows the RMSD plotted over the 250 tested values for the four different factors. The variants are colorcoded, and the best results are marked. The smallest values of the deviation are shown in the figure. The corresponding density profiles are shown in Fig. 1a.
The firn profile of ngt03C93.2 starts at a depth of approximately 1.3 m. Therefore, an appropriate surface density boundary condition must be found. As firn density profiles differ greatly, especially near the surface, the derivation of an appropriate surface density is difficult. Although Alley (1987) simulated the density starting at a depth of 2 m below the surface, we included this domain in our simulation so that we could apply transient surface forcing to our model. To find suitable values of the surface density, we included this parameter in our optimization. For each of the four variants and 250 factors C_{v}, we tested 21 values of the surface density between ρ_{0}=250 and ρ_{0}=450 kg m^{−3}, using steps of Δρ_{0}=10 kg m^{−3}. The best result was chosen. We used the method of testing 21 surface density values for all the analyzed firn profiles. We included profiles including measurements of the density at small depths. In this way, we established that the results are comparable. Profiles including nearsurface density values are, however, well represented. A total of $\mathrm{4}\times \mathrm{250}\times \mathrm{21}=\mathrm{21}\phantom{\rule{0.125em}{0ex}}\mathrm{000}$ simulations were performed for ice core ngt03C93.2 to find the optimal results shown in Fig. 1. The same procedure was applied to all 159 firn profiles analyzed in the study.
3.1 Firn profiles
To evaluate the description of grain boundary sliding given by Alley (1987), we used 159 firn profiles, 80 of which were retrieved in Greenland. The remaining 79 measurements were taken in Antarctica. The profiles are included in the SUMup snow density subdataset (Koenig and Montgomery, 2019). Individual references for all 159 firn profiles are listed in the Supplement. The dataset does not include the four profiles used in the study of Alley (1987), as the original data for these firn cores are unpublished. To obtain firn profiles suitable for this study from the dataset, we filtered it according to the following criteria.

Profiles must consist of at least 10 data points.

The overall length of each profile must exceed 3 m.

Profiles must start at a depth of less than 3 m below the surface.

The initial density of the profiles must not exceed ρ_{c}=550 kg m^{−3}.

The annual mean surface mass balance at the profile locations must be strictly positive.

Forcing data for at least 5 years must be available for the profile location.
Criteria 1 and 2 ensure the overall quality of the data, whereas criteria 3 and 4 ensure that the first stage of firn densification is included. As the model cannot handle melting and the study focuses on dry firn densification, the surface mass balance should be positive, as stated in the fifth criterion. However, a positive mean annual surface mass balance does not ensure that no melting occurs over the course of a year. It is not possible to distinguish between sites affected by melting and sites where no melting occurs using the data available for this study. The number of these sites, however, is expected to be small. Because of the method used, their influence on the overall result is therefore small. The forcing data come from the regional climate model RACMO2.3 (Van Wessem et al., 2014; Noël et al., 2015), which provides the surface mass balance. The model delivers data for the periods 1958 to 2016 and 1979 to 2016 for Greenland and Antarctica, respectively (see also Sect. 3.2). The density measurements used for model comparison should thus be retrieved during these periods. We used only datasets for which at least 5 years of forcing data are available. Furthermore, a number of density profiles were manually excluded from the filtered data. These profiles include those with very low spatial resolution, atypical profiles showing decreasing density with depth, and measurements with a surface density very close to the critical density of ρ_{c}=550 kg m^{−3}. As explained in Sect. 2.2 and illustrated in Fig. 1, we used only a certain domain for the comparison of the simulated and measured data. If this domain was found to be less than 2.5 m long for any of the tested variants of the constitutive equation, the firn profile in question was not analyzed further.
Figure 2 shows the locations from which the 159 density profiles were retrieved. The 80 measurements from Greenland are relatively uniformly distributed over the ice sheet. Coastal locations are not well covered owing to the requirement of a strictly positive surface mass balance. In the Antarctic, sites in East Antarctica are underrepresented. However, a wide variety of environments is covered, including the Filchner–Ronne Ice Shelf, the West Antarctic coast, and Dronning Maud Land.
3.2 Boundary conditions and forcing
To force the firn densification model, we need the surface values of density, temperature, accumulation rate, and grain radius at the locations of the 159 firn profiles. Although Alley (1987) used constant forcing, we followed the example of Arthern and Wingham (1998) and Goujon et al. (2003), who performed transient simulations.
As measured firn density profiles represent past climate conditions, the choice of forcing data in the proposed method is crucial. Uncertainties in the forcing will affect the simulation results and therefore the comparison with the measured firn profiles. Neither the model formulation nor the optimization scheme can compensate for these effects. We used data provided by the regional climate model RACMO2.3 (Van Wessem et al., 2014; Noël et al., 2015). Details on RACMO2.3, including the limitations of the model and the resulting data products, were presented by Van Wessem et al. (2014), Noël et al. (2015), and van Wessem et al. (2018). RACMO2.3 provides forcing data for the Greenland ice sheet covering 1958 to 2016. For Antarctica, the time period is shorter, running from 1979 to 2016. Data for the mean annual skin temperature and surface mass balance of the Greenland ice sheet are available at mean spatial resolutions of 11.3 and 1.0 km, respectively, for this study. The mean spatial resolutions of the mean annual skin temperature and surface mass balance for Antarctica are 8.0 and 28.5 km, respectively.
The time period for the transient simulation runs, as described in Sect. 2.2, is specified by the earliest data available from RACMO2.3 and the drilling date of the firn core under consideration. For ice core ngt03C93.2 (Wilhelms, 2000), which was retrieved in central Greenland in 1993, the simulation time covers the period from 1958 to 1993 (see Fig. 1c). Constant values of the surface temperature and surface mass balance for the preceding spinup period were calculated as mean values over this time range.
We present a second example to illustrate the effect of the temporal resolution of the forcing on the optimization results and why we used yearly averaged data provided by RACMO2.3. Figure 3 shows the depth density profile of the firn core retrieved at site 3 of the iSTAR Traverse in 2013 (Morris et al., 2017). The location of the site, at Pine Island Glacier in West Antarctica, is shown in Fig. 2. Instead of using forcing data from RACMO2.3, for this particular simulation we used ERA5Land monthly averaged data from 1981 to present (Muñoz Sabater, 2019; Hersbach et al., 2020), as they are freely available at monthly resolution. From these data, we computed the annual average data for a second simulation run. The forcing data at both resolutions are shown in Fig. 3c. Figure 3a shows the simulated firn profiles that best reproduce the measured density profile, which were identified using the optimization approach. The results obtained using the annually averaged forcing data are shown on the left, whereas those obtained using monthly averaged data from ERA5 are shown on the right. The data with higher temporal resolution reveal much more detail within the simulated firn density profiles. However, the aim of this study is not primarily to reproduce the analyzed measured firn profiles with the highest possible detail but to evaluate the constitutive relation of Alley (1987) using an optimization approach that identifies sitespecific optimal constitutive factors C_{v} (see Sect. 2.2). Figure 3b shows the RMSD of the simulated profiles from the measured profile over the range of tested optimization factors. Dashed lines represent simulations performed using the highresolution forcing data, whereas solid lines represent the annual averaged data. The difference between the optimization results is small. We therefore used the annual averaged data provided by RACMO2.3 suitable for this study, as the data cover a longer time period, especially for Greenland. As a result, we can analyze more firn profiles at greater detail. For ice core ngt03C93.2 (Wilhelms, 2000), the horizon of the year 1981, the earliest forcing available in ERA5, lies at a depth of approximately 5 m below the surface, as shown in Fig. 1a. The horizon of the earliest forcing available from RACMO2.3, the year 1958, is located at a depth of approximately 11 m below the surface. A much greater portion of the simulated firn profile is therefore affected by surface forcing. Furthermore, the use of yearly averaged data requires less overhead.
As noted in Sect. 2.2, we used 21 surface density values in the range 250 kg m^{−3} ≤ ρ_{0} ≤ 450 kg m^{−3} for every tested firn profile. The value that yielded the best result was then used for further analysis. Owing to a lack of relevant data, and for simplicity and ease of comparison, the grain radius at the surface was assumed to be the same at all locations and to be constant over time. We used a grain radius of r_{0}=0.5 mm on the basis of the measurements and empirical relation of Linow et al. (2017) and the assumption of Arthern and Wingham (1998). As climatic conditions, and therefore the surface grain size, vary among the investigated locations, this choice is a simplification. Owing to the use of the optimization approach, this parameter has a less significant effect, as it can be understood as a constant part of the grain radius, which is the same for every analyzed firn profile.
Figure 4 illustrates the range and distribution of surface boundary conditions at the investigated sites. The locations in Greenland show a higher mean surface temperature and surface mass balance than those in Antarctica. The surface density is higher at the Antarctic locations than at those in Greenland. Overall, a wide variety of typical climatic conditions for both ice sheets is covered.
3.3 Distribution and effects of input data
Ice core ngt03C93.2 (Fig. 1a) is an example of a highresolution density measurement showing extensive smallscale layering. Only a few of the 159 firn profiles are of such high quality and include this type of layering. Although our proposed model works at high temporal and spatial resolution, it does not cover layering, as shown in Fig. 1a. The density profile retrieved at site 3 of the iSTAR traverse (Morris et al., 2017) (Fig. 3) illustrates that the model, if it is forced with data of higher temporal resolution, still does not cover the measured density variability. Smallscale layering of firn appears to be driven by a number of different processes (Hörhold et al., 2011). An extension of the model to cover such processes may be introduced in the future. We would prefer the approach of Freitag et al. (2013), who introduced the concept of impuritycontrolled densification. Forcing data for this model are not globally available. However, the model in its current form does reproduce the mean density well, as demonstrated in Fig. 1a.
In this study, we focused on the initial stage of firn densification and the uppermost part of the firn column. We limited the model domain to a maximum depth of −25 m below the surface. Furthermore, we did not develop the density further after the critical density of ρ_{c}=550 kg m^{−3} was reached. This choice raises the question of whether the use of a Neumann boundary condition set to zero at the profile base to solve for the temperature is justified for this particular model setup (Sect. A4). The critical density of ρ_{c}=550 kg m^{−3} is usually reached within the upper 10 m of the firn column. The temperature within this domain is affected by surface conditions, which are covered by the forcing data. At greater depths, the temperature corresponds to the mean annual surface temperature and changes very little (e.g., Cuffey and Paterson, 2010, 399 pp.). There have been few highresolution temperature measurements of firn. Orsi et al. (2017) published a temperature profile of a 147 m length from borehole NEEM2009S1, which shows a temperature difference of little more than 1 K over the entire profile. Vandecrux et al. (2021) obtained the same result. This small temperature change at depths greater than approximately 10 m below the surface has little effect on our model approach. To test the temperature dependence of the optimization approach, we performed the simulations again using the sitespecific surface temperature forcing plus or minus 1 K. The effect on the optimal factors C_{v} depends on the variant of the constitutive relation, but it is generally small. For example, the greatest difference between the optimal factors obtained using the correct and adjusted forcing for ice core ngt03C93.2 (Wilhelms, 2000) were found using Variant 1. This variant yields a value of $\mathrm{\Delta}{C}_{{\mathrm{v}}_{\mathrm{1}},\mathrm{max}}=\mathrm{0.02}\times {\mathrm{10}}^{\mathrm{4}}$ K s^{2} kg^{−1}, which corresponds to twice the sampling space of the factors ${C}_{{\mathrm{v}}_{\mathrm{1}}}$ used in the optimization. Variant 4 yields the same results. Because of this small effect overall, the stronger effect of the surface temperature in the investigated domain, and the limitation of the comparison to the domain actually affected by surface forcing, a Neumann boundary condition at the profile base is justified despite the low depth.
Figure 5 shows the distribution of the RMSD calculated from the simulated firn profiles that best reproduce the measured profiles. The four plots in the figure represent the four variants of the constitutive equation for grain boundary sliding (Eqs. 2 to 5), as described in Sect. 2.2. Additionally, the median values of the data are shown. The distributions of the deviation for the four variants show little difference. The median values differ in a small range. Variants 2 and 3 show the smallest and largest values, respectively. The use of the modification introduced by Bréant et al. (2017) in the constitutive equation improves the agreement between the simulated and measured firn profiles by 6.2 %–6.8 %. To put the values in perspective, the deviation of the four bestfitting modeled firn profiles from ice core ngt03C93.2 (Wilhelms, 2000) displayed in Fig. 1 is approximately 28 kg m^{−3}. That is, more than half of the simulations show even better agreement with the corresponding firn density profiles than this one. An overview of the deviation among all 159 measured firn profiles and the corresponding best simulation results can be found in the Supplement.
The range of factors contributing to the bestfitting firn profile is smaller for variants 2 and 4 of the constitutive equation, which use the modification of Bréant et al. (2017), compared to variants 1 and 3, as shown in the box plots in Fig. 6. In contrast to factors ${C}_{{\mathrm{v}}_{\mathrm{1}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{2}}}$, factors ${C}_{{\mathrm{v}}_{\mathrm{3}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{4}}}$ incorporate the Arrhenius factor from the original description of grain boundary sliding given by Alley (1987). Therefore, direct comparison between the two groups of factors is not possible. The quartile coefficient of dispersion of the four variants, which is shown on the right side of Fig. 6, is a relative measure of the scatter of the values. The coefficient reveals that the factors obtained for variants 1 and 2 are defined in a narrower range than those of variants 3 and 4. All four sets of factors show slightly nonuniform distributions that tend toward smaller values.
To check for possible mean surface temperature dependence of the 159 factors found by optimization, these values are plotted against each other in Fig. 7. The values of the mean surface temperature were calculated from the forcing data for each firn profile site. The left plot illustrates the results for variants 1 and 2 in red and blue, whereas the right plot shows the factors for variants 3 and 4 using green and purple markers, respectively. The legend shows the Pearson correlation coefficient r_{Pearson}, a measure of the linear correlation of the two variables, and the distance correlation dCor (Székely et al., 2007). The distance correlation was designed by Székely et al. (2007) to overcome problems with the Pearson correlation coefficient. It describes the correlation of two vectors and is not limited to linear dependence. It is defined between zero and one, where zero indicates the independence of the variables.
The correlation between the resulting factors and the mean surface temperature is higher compared to the correlation with other properties. This statement is especially true for factors ${C}_{{\mathrm{v}}_{\mathrm{3}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{4}}}$ of variants 3 and 4, respectively. However, the Pearson correlation coefficient indicates only a linear correlation. The scatter of factors ${C}_{{\mathrm{v}}_{\mathrm{3}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{4}}}$ with respect to mean surface temperature might resemble a higherorder function, as indicated by the higher values of the distance correlation.
The values of the Pearson correlation coefficient and distance correlation shown in Fig. 8, which shows the factors with respect to mean surface mass balance, are higher than those in Fig. 7.
Like the mean surface temperature, the mean surface mass balance was calculated from the forcing data used during the simulations. The correlation coefficients for the results obtained using variants 3 and 4 exceed those for the results obtained using the variants of the constitutive equation incorporating the Arrhenius factor explicitly. Variant 3 shows the best indication of a relationship between surface mass balance and the factor yielding the density profile that best reproduces the corresponding field measurement.
In Fig. 8, a striking feature appears around a mean surface mass balance of 0.4 m weq. a^{−1}. The values of the factors yielding the simulated density profiles that best reproduce the measured profiles show a wide range for all four variants of the constitutive relation. These profiles are part of a study in western Greenland by Harper et al. (2012). The study region is relatively small, which explains the similar climatic conditions. Although the mean annual surface mass balance is positive, melting occurs throughout the year and affects the firn density.
Using the four variants of the constitutive relation of Alley (1987) (Eqs. 2 to 5), we generated density profiles that were in good agreement with most of the 159 density measurements. Uncertainties may result from the forcing, limited knowledge of the initial grain radius, and the poor constraint on the density at the surface. As measured firn densities represent past climate conditions, unrealistic forcing always results in a mismatch between the simulated and measured firn properties, independent of the optimization approach and physical model. However, when the proposed optimization scheme was used, the differences between the simulated and measured density profiles show distinct minima, indicating that the forcing represents the climatic history of the firn profiles, in principle.
The optimization scheme produces sitespecific values of the factors ${C}_{{\mathrm{v}}_{\mathrm{1}}}$ to ${C}_{{\mathrm{v}}_{\mathrm{4}}}$. As the optimization process is unique to each analyzed site and constitutive relation variant, the four simulated density profiles for each site are very similar, as illustrated in Fig. 5. This feature allows us to compare the factors obtained using the four variants. The differences between the factors do not arise primarily from differences in the simulated density profiles but reflect the differences in the variants; consequently, the results are almost the same.
However, owing to the nature of the optimization scheme, possible errors in the forcing data or other parameters such as the activation energy used to calculate the grain radius (see Eq. A11) are also included in the sitespecific optimal factors ${C}_{{\mathrm{v}}_{\mathrm{1}}}$ to ${C}_{{\mathrm{v}}_{\mathrm{4}}}$. If, for example, the surface temperature forcing at a site consistently deviates by 5 K during the simulation, this error can be compensated for by adjusting the specific optimization factor C_{v} (see Sect. 3.3). However, the large number of analyzed firn profiles compensates for such random errors. Systematic errors in the forcing data (van Wessem et al., 2018) cannot be identified. Therefore, improving the temporal resolution of the forcing and covering longer periods could result in better and more detailed simulation results, as shown in Fig. 3.
This study analyzed only dry firn densification. The current model cannot handle melting. We accommodate this feature by setting the annual mean surface mass balance at the investigated sites to be strictly positive (Sect. 3.1). However, this limitation means that we cannot ensure that no melting occurs over the course of a year. The results shown in Fig. 8 illustrate how this limitation affects the optimization results. The limitation is problematic, especially in recent years, when melting occurred over almost the entire Greenland ice sheet (e.g., Nghiem et al., 2012). The simulation of meltwater percolation through the firn and its interaction with firn densification is important, especially in the upper part of the firn body (e.g., Vandecrux et al., 2020). The proposed method could be improved by the application of this model approach in future investigations. However, we identified some correlations between the optimization results and the surface mass balance.
Note, again, that we have not investigated whether grain boundary sliding is indeed the dominant process during the first stage of firn densification. We assess whether a process with a functional dependence on density, firn overburden stress, temperature, and grain radius represents the observed density profiles well. Any other deformation process with the same functional dependence would be equally well suited. Nevertheless, by maintaining the general structure of the constitutive relation of Alley (1987), we conclude that this description of grain boundary sliding is a good basis for a physically based model describing firn densification up to the critical density, ρ_{c}=550 kg m^{−3}.
Comparing the results for variants 1 and 2, as well as those for variants 3 and 4, we found that an adjustment of $\left(\mathrm{1}\frac{\mathrm{5}}{\mathrm{3}}\mathit{\rho}/{\mathit{\rho}}_{\mathrm{ice}}\right)$ to $\left(\mathrm{1}+\frac{\mathrm{0.5}}{\mathrm{6}}\frac{\mathrm{5}}{\mathrm{3}}\mathit{\rho}/{\mathit{\rho}}_{\mathrm{ice}}\right)$ results in better reproduction of the measured density profiles. This result must be reviewed carefully in terms of the study design and the background of a physicsbased model describing firn densification. As Alley (1987) pointed out, grain boundary sliding might be the dominant process driving firn densification at low densities, but it is presumably not the only one. The constitutive law of Alley (1987) is designed such that the densification due to grain boundary sliding becomes zero at a density of ρ_{c}=550 kg m^{−3}, owing to the densest packing of spheres obtained at that density and increasing accommodation incompatibilities. The modification of Bréant et al. (2017) changes this behavior such that grain boundary sliding vanishes at a density of approximately ${\mathit{\rho}}_{\mathrm{c}}^{*}=\mathrm{596}$ kg m^{−3}, which could have advantages for the transition into the next stage of firn densification. We suggest that a simultaneous decrease in grain boundary sliding and increase in one or more other processes would provide a good characterization. Specifically, dislocation creep drives densification at higher density (Maeno and Ebinuma, 1983) because of increasing stresses. The onset of dislocation creep at densities below ρ_{c}=550 kg m^{−3} would not necessarily affect the entire bulk firn matrix but increasing volume fractions of the porous matrix.
Variants 1 and 2 of the constitutive relation of Alley (1987) incorporate the Arrhenius factor for boundary diffusion, D_{BD}, from the description of the bond viscosity given by Raj and Ashby (1971). In the formulation of variants 3 and 4, we neglected the Arrhenius factor. As shown in Fig. 5, the difference in the resulting RMSD is small regardless of whether D_{BD} is considered. This result is reasonable, as we determined individual factors C_{v} for each site. The similarity of the results allows us to compare the factors ${C}_{{\mathrm{v}}_{\mathrm{1}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{2}}}$, which were determined using the Arrhenius factor D_{BD}, to factors ${C}_{{\mathrm{v}}_{\mathrm{3}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{4}}}$, which were obtained using the variants of the constitutive relation without the Arrhenius factor. As D_{BD} is a function of temperature, it is reasonable to examine the dependence of the factors on the mean surface temperature, which is shown in Fig. 7.
The factors ${C}_{{\mathrm{v}}_{\mathrm{3}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{4}}}$, which were determined using the variants without D_{BD}, show a stronger dependence on the mean surface temperature than the factors ${C}_{{\mathrm{v}}_{\mathrm{1}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{2}}}$. In addition, the factors ${C}_{{\mathrm{v}}_{\mathrm{1}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{2}}}$ show less dispersion than ${C}_{{\mathrm{v}}_{\mathrm{3}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{4}}}$, as shown in Fig. 6. The inclusion of D_{BD} in the constitutive relation improves the determination of these factors. It is therefore a meaningful description within the constitutive relation. Although the inclusion of the Arrhenius factor improves the determination of ${C}_{{\mathrm{v}}_{\mathrm{1}}}$ and ${C}_{{\mathrm{v}}_{\mathrm{2}}}$, we still see some dependence on the mean surface temperature. A better determination of the parameters of the Arrhenius factor may resolve this dependence. If this is not the case, another dependence on the temperature may be introduced to improve the constitutive relation for grain boundary sliding.
We interpret the dependence on surface mass balance as indicating that the load is currently not represented well. The stress is represented by a secondorder tensor. A firn column represented by a onedimensional modeling approach would be surrounded by neighboring firn columns, resulting in lateral confinement that limits deformation in the horizontal direction. The horizontal components of the stress tensor are not zero. As firn is a compressible material, the determination of these horizontal stress components is not trivial. The frequently used term “overburden pressure” is misleading, as the mechanical pressure is defined as the spherical part of the Cauchy stress tensor (e.g., Haupt, 2002, p. 301) and is not, in general, identical to the normal stress in the vertical direction. With increasing depth, the magnitude of the horizontal stress components and their effects would increase. Modeling approaches that consider the full stress tensor were reported by Greve and Blatter (2009), Salamantin et al. (2009), and Meyer and Hewitt (2017). It might be worth using the constitutive relation for grain boundary sliding of Alley (1987) in such a modeling context. This approach will not necessarily result in a full threedimensional model, as the problem can be formulated assuming axial symmetry, which would require adjustment of the constitutive relation. A more extensive interpretation of the factors that yielded the best agreement is therefore challenging. The determination of a single factor C_{v} for one or all of the variants is not useful. It would not result in better simulation results compared to other published firn densification models. The sitespecific values of the factors determined using the presented optimization approach simply show the differences in the variants of the constitutive relation.
Using variants of the constitutive relation for grain boundary sliding of Alley (1987) and an efficient optimization scheme, we reproduced 159 firn density profiles reasonably well. Thus, we conclude that the description of grain boundary sliding introduced by Alley (1987) is appropriate for the simulation of firn densification at low density.
The modification of the constitutive relation of Bréant et al. (2017) results in slightly better simulation results when only the first stage of firn densification is considered. Further factors, including the transition from the first to the second stage, must be considered to answer questions regarding in which domain and to what extent different processes drive firn densification.
In our optimization approach, we use a single factor representing various model parameters and search for the factor value resulting in the best agreement between the simulated and measured firn profiles. Consequently, the sitespecific simulation results are independent of the possibly deficient model parameters, which are considered collectively. It is not possible to derive a distinct value for the factor representing the climatic conditions at all locations of the investigated firn profiles. We found a linear dependence of the factors on the sitespecific surface mass balance (Fig. 8). By contrast, we did not find a clear linear dependence on temperature (Fig. 7), but the results show that a sitespecific parameter is required.
As the amount of surface accumulation affects the load conditions, we assume it is not represented well in the model. Unlike other firn densification models, the physical simulation of grain boundary sliding does not depend directly on the surface mass balance but depends on the actual stress. Further interpretation of the obtained factors is difficult using the presented simulation setup. The description of grain boundary sliding given by Alley (1987) could be improved by using a higherdimensional approach including the horizontal components of the stress tensor. Modeling approaches of this kind include those of Greve and Blatter (2009), Salamantin et al. (2009), and Meyer and Hewitt (2017).
We would like to emphasize that optimization of any type is possible only because of the enormous efforts of the SUMup team (Koenig and Montgomery, 2019), which has made a vast number of firn core data available. These data are strategically crucial for advances in firn densification modeling, reinforcing the recommendations of FirnMICE (Lundin et al., 2017) for enhanced efforts toward physically based models.
A1 Numerical treatment of densification
All model equations are solved on an adapting onedimensional grid that is updated at every time step. The approach is based on an updated Lagrangian description, where the update velocity of the grid is the material flow velocity. This results in materialfixed coordinates. The Lagrangianlike description affords very high spatial and temporal resolution in the simulations. It can be shown that by integrating the local Eulerian form of the mass balance in one dimension over a material control volume with moving boundaries z_{1}(t) and z_{2}(t), we obtain (Ferziger and Perić, 2002, p. 374)
Here ρ represents the density, z is the vertical coordinate, t is the time, and v is the material flow velocity, whereas v_{b} represents the grid velocity or the velocity of the integration boundary. When the grid velocity equals the material flow velocity (v_{b}=v), the second term of Eq. (A1), which describes the advection, vanishes. The resulting equation is equal to the Lagrangian form of the mass balance (Ferziger and Perić, 2002, p. 374). On a onedimensional grid consisting of a number of grid points, as illustrated in Fig. A1, we denote the grid point velocity as v_{b}, which is equal to the flow velocity v_{b}≡v. The locations of all grid points are updated at each time step by integrating the grid point velocity v_{b} using a forward Euler scheme. Thus, advection is represented entirely by the adapting grid.
The grid point velocity v_{b} is calculated using the constitutive equation for grain boundary sliding, as described in Sect. 2.1, and the definition of the strain rate in one dimension. The description of grain boundary sliding provides the strain rate in the vertical direction along the grid, ${\dot{\mathit{\epsilon}}}_{zz}$, as a function of vertical stress σ_{zz}, density ρ, temperature T, and grain radius r:
The strain rate of a material line element can be defined as the spatial derivative of the velocity, as shown in Eq. (A2) (Haupt, 2002, 32–38). On a onedimensional grid, which is defined as a number of grid points, the space between neighboring grid points can be considered a material line element (see Fig. A1). Therefore, the grid point velocity v_{b} is easily computed by integrating the strain rate ${\dot{\mathit{\epsilon}}}_{zz}$ in the vertical direction along the length of the grid cell:
To implement Eq. (A3), an integration constant determined by a suitable boundary condition is needed. It is reasonable to apply a Dirichlet boundary condition that forces the grid point velocity v_{b} to be zero at either the top or the base of the computational domain to represent a fixed reference point at the top or the base of the modeled firn profile, respectively. All other points defining the adapting grid are moving with respect to this anchor point. In this study, we placed the anchor point at the base of the simulated firn profiles (z_{0} in Fig. A1). The depth coordinates of the profiles shown in the figures were adjusted for better readability.
To represent accumulation at the top of a simulated firn profile, an inflow boundary condition must be implemented. To this end, at each time step an additional grid point is added at the top of the grid. Its coordinate within the grid ${z}_{n+\mathrm{1}}(t+\mathrm{\Delta}t)$ is calculated as
where a_{0} is the accumulation rate given as meters per second (m s^{−1}) water equivalent, Δt is the length of the time step, ρ_{water} is the density of water, and ρ_{0} is the density of the deposited snow. The expression ${z}_{n+\mathrm{1}}(t+\mathrm{\Delta}t)$ is the position of the firn surface at the last time step z_{n} plus the thickness of the firn layer deposited during the last time step. This thickness is defined by the time step Δt and the timedependent accumulation rate a_{0}(t). We use an accumulation rate given in the unit of meter water equivalent per second, which must be converted to the unit of meter firn equivalent taking the sitespecific surface firn density ρ_{0} and water density ρ_{water} into account.
In this process, the number of grid points increases. Therefore, grid points are removed from the model domain at its base when a maximum number is reached.
A2 Stress
To evaluate the stress in the vertical direction, σ_{zz}, we use the local form of the static linear momentum balance in its Eulerian description. We neglect the acceleration component, as the changes in velocity are assumed to be small; thus, we have
The stress σ_{zz} can be easily calculated by integrating the product of density ρ and acceleration due to gravity g along the simulated profile. We assumed that the surface of the profile is tractionfree.
A3 Density
As noted in Sect. A1 and illustrated by Eq. (A1), the change in density integrated over a control volume with respect to time must be zero. That is, the mass in a control volume cannot change. As the positions of the grid points, and thus the material control volume, do change, the density changes accordingly. The evaluation of the density integral over a control volume at two time steps yields
As the space between two neighboring grid points can be understood as a material line element (Haupt, 2002, 32–38), Eq. (A6) can be rewritten in the form of Eq. (A7), where $\mathrm{}\mathrm{d}Z\mathrm{}=\mathrm{}{z}_{\mathrm{2}}\left(t\right){z}_{\mathrm{1}}\left(t\right)\mathrm{}$ and $\mathrm{}\mathrm{d}z\mathrm{}=\mathrm{}{z}_{\mathrm{2}}(t+\mathrm{\Delta}t){z}_{\mathrm{1}}(t+\mathrm{\Delta}t)\mathrm{}$ are the lengths of a material line element in the reference configuration and its image in the current configuration, respectively:
By rearranging Eq. (A7), we obtain the formulation of the density change with time depending on the definition of the strain ε_{zz} in one dimension (Haupt, 2002, p. 34):
The evolution of the density can therefore be computed by integrating the strain rate ${\dot{\mathit{\epsilon}}}_{zz}$ (Sect. 2.1) over time.
A4 Temperature
As noted in Sect. A1, all advection in the model domain is represented by the moving grid. Therefore, the description of the temperature evolution reduces to simple heat diffusion:
Following Paterson (1994), we assume a constant heat capacity of c_{p}=2009 J kg^{−1} K^{−1}, and following the example of Zwinger et al. (2007), we assume the densitydependent thermal conductivity described by Sturm et al. (1997) as
The temperature profile is initialized using a constant mean value computed from the sitespecific surface forcing. To solve for temperature, a Neumann boundary condition is used at the lower boundary at the profile base. The vertical gradient of the temperature is set to be zero.
A5 Grain radius
Alley (1987) used measured grain size data to fit his simulation results to four firn profiles. As information about grain size is sparse, we use a modeling approach to describe the grain radius. The evolution of the grain radius r is simulated using the wellknown description of Stephenson (1967) and Gow (1969), as given by Arthern et al. (2010). Stephenson (1967) and Gow (1969) described the grain size in terms of the mean crosssectional area. By contrast, Arthern et al. (2010) assumes the mean crosssectional area to be $A=(\mathrm{2}/\mathrm{3})\mathit{\pi}{r}^{\mathrm{2}}$ and formulates the grain growth rate as
This formulation enables simple calculation of the grain radius r. The values of the activation energy, E_{g}=42.4 kJ mol^{−1}, and prefactor, ${k}_{\mathrm{0}}=\mathrm{1.3}\times {\mathrm{10}}^{\mathrm{7}}$ m^{2} s^{−1}, of the Arrhenius factor are based on data reported by Paterson (1994) and were also adapted from Arthern et al. (2010). Unlike Arthern et al. (2010), we do not use the mean annual temperature but the actual temperature T(z,t) along the simulated profile. To solve Eq. (A11), suitable boundary conditions must be provided. We used a constant surface grain radius. See Sect. 3.2 for further information on the boundary conditions.
A6 Age
For comparison (Sect. 3.2) and general interest, the firn age χ was also simulated. Because advection is represented by an adapting grid, the description is very simple. It is calculated as
Newly deposited snow has an age of zero, which is represented by a Dirichlet boundary condition. The age discretized at the grid points then increases according to the time step.
A7 Time discretization
Time is discretized using constant time steps. For this study, a value of 48 time steps per year was found to be a good compromise between economical simulation and desirable resolution and was used in all simulations. The grid resolution depends on the time step, as a new grid point is generated at each time step to represent accumulation, as described in Sect. A1 and shown in Eq. (A4). Timedependent properties such as density, temperature, grain radius, and age were developed using an explicit Euler scheme.
The code used to simulate firn profiles and to perform the optimization is available via the openaccess repository zenodo at https://doi.org/10.5281/zenodo.5771956 (Schultz, 2021).
The supplement related to this article is available online at: https://doi.org/10.5194/tc161432022supplement.
TS developed the numerical approach and code and conducted and analyzed all simulations. All authors jointly developed the concept of the modeling approach, discussed the results, and wrote the manuscript.
The contact author has declared that neither they nor their coauthors have any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We acknowledge the fantastic community effort SUMup, which provides a database of firn core data. We thank Sepp Kipfstuhl (AWI) and Johannes Freitag (AWI) for discussions on firn densification.
This research has been supported by Deutsche Forschungsgemeinschaft in the priority program SPP1158 (grant no. 403642112).
This paper was edited by Kaitlin Keegan and reviewed by C. Max Stevens and one anonymous referee.
Alley, R. B.: Firn Densification by GrainBoundarySliding: A First Model, Journal de Physique, 48, C1249–C1256, https://doi.org/10.1051/jphyscol:1987135, 1987. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac, ad, ae, af, ag, ah, ai, aj, ak, al, am, an, ao, ap, aq
Amante, C. and Eakins, B. W.: ETOPO1 1 ArcMinute Global Relief Model: Procedures, Data Sources and Analysis, NOAA Technical Memorandum NESDIS NGDC24, National Geophysical Data Center, NOAA, https://doi.org/10.7289/V5C8276M, 2009. a
Anderson, D. L. and Benson, C. S.: Ice and Snow: Properties, Processes and Applications, in: chap. The densification and diagenesis of snow, MIT Press, Cambridge, MA, 391–411, ISBN 13 9780262110099, 1963. a
Arnaud, L., Barnola, J. M., and Duval, P.: Physical modeling of the densification of snow/firn and ice in the upper part of polar ice sheets, Physics of Ice Core Records, 285–305, available at: http://hdl.handle.net/2115/32472 (last access: 13 January 2022), 2000. a, b, c
Arndt, J. E., Schenke, H. W., Jakobsson, M., Nitsche, F. O., Buys, G., Goleby, B., Rebesco, M., Bohoyo, F., Hong, J., Black, J., Greku, R., Udintsev, G., Barrios, F., ReynosoPeralta, W., Taisei, M., and Wigley, R.: The International Bathymetric Chart of the Southern Ocean (IBCSO) Version 1.0 – A new bathymetric compilation covering circumAntarctic waters, Geophys. Res. Lett., 40, 3111–3117, https://doi.org/10.1002/grl.50413, 2013. a
Arthern, R. J. and Wingham, D. J.: The Natural Fluctuations of Firn Densification and their Effect on the Geodetic Determination of Ice Sheet Mass Balance, Climatic Change, 40, 605–624, https://doi.org/10.1023/A:1005320713306, 1998. a, b, c, d, e, f
Arthern, R. J., Vaughan, D. G., Rankin, A. M., Mulvaney, R., and Thomas, E. R.: In situ measurements of Antarctic snow compaction compared with predictions of models, J. Geophys. Res., 115, F03011, https://doi.org/10.1029/2009JF001306, 2010. a, b, c, d, e
Bader, H.: Sorge's Law of Densification of Snow on High Polar Glaciers, J. Glaciol., 2, 319–323, https://doi.org/10.3189/S0022143000025144, 1954. a
Bréant, C., Marinterie, P., Orsi, A., Arnaud, L., and Landais, A.: Modelling firn thickness evolution during the last deglaciation: constraints on sensitivity to temperature and impurities, Clim. Past., 13, 833–853, https://doi.org/10.5194/cp138332017, 2017. a, b, c, d, e, f, g, h, i, j, k, l, m
Cuffey, K. M. and Paterson, W. S. B.: The Physics of Glaciers, 4th Edn., ButterworthHeinemann, ISBN 13 9780123694614, 2010. a
Ferziger, J. H. and Perić, M.: Computational Methods for Fluid Dynamics, 3. rev. Edn., Springer, Berlin, Heidelberg, New York, Barcelona, Hong Kong, London, Milan, Paris, Tokyo, ISBN 3540420746, 2002. a, b
Fourtenau, K., GilletChaulet, F., Martinerie, P., and Faïn, X.: A MicroMechanical Model for the Transformation of Dry Polar Firn Into Ice Using the LevelSet Method, Front. Earth. Sci., 8, 101, https://doi.org/10.3389/feart.2020.00101, 2020. a
Freitag, J., Kipfstuhl, S., Laepple, T., and Wilhelms, F.: Impuritycontrolled densification: a new model for stratified polar firn, J. Glaciol., 59, 1163–1169, https://doi.org/10.3189/2013JoG13J042, 2013. a
Goujon, C., Barnola, J.M., and Ritz, C.: Modeling the densification of polar firn including heat diffusion: Application to closeoff characteristics and gas isotopic fractionation for Antarctica and Greenland sites, J. Geophys. Res., 108, 4792, https://doi.org/10.1029/2002JD003319, 2003. a, b, c, d
Gow, A. J.: On the Rates of Growth of Grains and Crystals in South Polar Firn, J. Glaciol., 8, 241–252, https://doi.org/10.3189/S0022143000031233, 1969. a, b, c
Greve, R. and Blatter, H.: Dynamics of Ice Sheets and Glaciers, in: Advances in Geophyiscal and Environmental Mechanics and Mathematics, SpringerVerlag, Berlin, Heidelberg, https://doi.org/10.1007/9783642034152, 2009. a, b
Harper, J., Humphrey, N., Pfeffer, W. T., Brown, J., and Fettweis, X.: Greenland icesheet contribution to sealevel rise buffered by meltwater storage in firn, Nature, 491, 240–243, https://doi.org/10.1038/nature11566, 2012. a
Haupt, P.: Continuum Mechanics and Theory of Materials, SpringerVerlag, Berlin, Heidelberg, New York, ISBN 13 9783540431114, 2002. a, b, c, d
Herron, M. M. and Langway, C. C.: Firn Densification: An Empirical Model, J. Glaciol., 25, 373–385, https://doi.org/10.3189/S0022143000015239, 1980. a
Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.N.: The ERA5 gloabal reanalysis, Q. J. Roy. Meterol. Soc., 146, 1999–2049, https://doi.org/10.1002/qj.3803, 2020. a, b
Hörhold, M. W., Kipfstuhl, S., Wilhelms, F., Freitag, J., and Frenzel, A.: The densification of layered polar firn, J. Geophys. Res.Earth, 116, F01001, https://doi.org/10.1029/2009JF001630, 2011. a
Ignat, M. and Frost, H. J.: Grain Boundary Sliding in Ice, J. Phys. Colloques, 48, C1189–C1195, https://doi.org/10.1051/jphyscol:1987127, 1987. a, b
Itagaki, K.: SelfDiffusion in Single Crystals of Ice, J. Phys. Soc. Jpn., 19, 1081, https://doi.org/10.1143/JPSJ.19.1081, 1964. a
Johnson, J. B. and Hopkins, M. A.: Identifying microstructural deformation mechanisms in snow using discreteelement modeling, J. Glaciol., 51, 432–442, https://doi.org/10.3189/172756505781829188, 2005. a
Kinosita, S.: Compression of Snow at Constant Speed, in: Physics of Snow and Ice: proceedings, 1, 911–927, International Conference on Low Temperature Science, I. Conference on Physics of Snow and Ice, II. Conference on Cryobiology, 14–19 August 1966, Sapporo, Japan, 1967. a
Koenig, L. and Montgomery, L.: Surface Mass Balance and Snow Depth on Sea Ice Working Group (SUMup) snow density subdataset, Greenland and Antarctica, 1950–2018, Tech. rep., Arctic Data Center [data set], https://doi.org/10.18739/A26D5PB2S, 2019. a, b, c
Ligtenberg, S. R. M., Helsen, M. M., and van den Broeke, M. R.: An improved semiempirical model for the densification of Antarctic firn, The Cryosphere, 5, 809–819, https://doi.org/10.5194/tc58092011, 2011. a
Linow, S., Hörhold, M. W., and Freitag, J.: Grainsize evolution of polar firn: a new empirical grain growth parameterization based on Xray microcomputer tomography measurements, J. Glaciol., 58, 1245–1252, https://doi.org/10.3189/2012JoG11J256, 2017. a
Lundin, J. M. D., Stevens, C. M., Arthern, R., Buizert, C., Orsi, A., Ligtenberg, S. R. M., Simonsen, S. B., Cummings, E., Essery, R., Leahy, W., Harris, P., Helsen, M. M., and Waddington, E. D.: Firn Model Intercomparison Experiment (FirnMICE), J. Glaciol., 63, 401–422, https://doi.org/10.1017/jog.2016.114, 2017. a
Maeno, N. and Ebinuma, T.: Pressure Sintering of Ice and Its Implications to the Densification of Snow at Polar Glaciers and Ice Sheets, J. Phys. Chem., 87, 349–365, https://doi.org/10.1021/j100244a023, 1983. a, b, c, d, e
Meyer, C. R. and Hewitt, I. J.: A continuum model for meltwater flow through compacting snow, The Cryosphere, 11, 2799–2813, https://doi.org/10.5194/tc1127992017, 2017. a, b
Miller, H. and Schwager, M.: Accumulation rate and stable oxygen isotpoe ratios of ice core ngt03C93.2 form the North Greenland Traverse, PANGAEA [data set], https://doi.org/10.1594/PANGAEA.218274, 2004. a
Morris, E. M., Muvaney, R., Arthern, R. J., Davies, D., Gurney, R. J., Lambert, P., De Rydt, J., Smith, A. M., Tuckwell, R. J., and Winstrup, M.: Snow Densification and Recent Accumulation Along the iSTAR Traverse, Pine Island Glacier, Antarctica, J. Geophys. Reas.Earth, 122, 2284–2301, https://doi.org/10.1002/2017JF004357, 2017. a, b, c, d
Muñoz Sabater, J.: ERA5Land hourly data from 1981 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CFS) [data set], https://doi.org/10.24381/cds.e2161bac, 2019. a, b
Nghiem, S. V., Hall, D. K., Mote, T. L., Tedesco, M., Albert, M. R., Keegan, K., Shuman, C. A., DiGirolamo, N. E., and Neumann, G.: The extreme melt across the Greenland ice sheet in 2012, Geophyis. Res. Lett., 39, L20502, https://doi.org/10.1029/2012GL053611, 2012. a
Noël, B., van de Berg, W. J., van Meijgaard, E., Kuipers Munneke, P., van de Wal, R. S. W., and van den Broeke, M. R.: Evaluation of the updated regional climate model RACMO2.3: summer snowfall impact on the Greenland Ice Sheet, The Cryosphere, 9, 1831–1844, https://doi.org/10.5194/tc918312015, 2015. a, b, c, d, e, f, g
Orsi, A. J., Kawamura, K., MassonDelmotte, V., Fettweis, X., Box, J. E., DahlJensen, D., Clow, G. D., Landais, A., and Severinghaus, J. P.: The recent warming trend in North Greenland, Geophys. Res. Lett., 44, 6235–6243, https://doi.org/10.1002/2016GL072212, 2017. a
Paterson, W. S. B.: The Physics of Glaciers, 3rd Edn., ButterworthHeinemann, Oxford, Birlington, ISBN 9780750647427, 1994. a, b
Raj, R. and Ashby, M. F.: On Grain Boundary Sliding and Diffusional Creep, Metallurg. Trans., 2, 1113–1127, https://doi.org/10.1007/BF02664244, 1971. a, b, c, d
Robin, G. D. Q.: Glaciology III: Seismic Shooting and Related Investigations, NorwegianBritishSwedish Antarctic Expedition, 1949–52, Scientific Results Vol. 5, Norsk Polarinstitutt, Oslo, 1958. a
Roscoat, S. R. D., King, A., Phillip, A., Reischig, P., Ludwig, W., Flin, F., and Meyssonnier, J.: Analysis of Snow Microstructure by Means of XRay Diffraction Contrast Tomography, Adv. Eng. Mater., 13, 128–135, https://doi.org/10.1002/adem.201000221, 2010. a
Salamantin, A. N., Lipenkov, V. Y., Barnola, J. M., Hori, A., Duval, P., and Hondoh, T.: Snow/Firn Densification in Polar Ice Sheets, in: Physics of Ice Core Records II: Papers collected after the 2nd International Workshop on Physics of Ice Core Records, 2–6 February 2007, Sapporo, Japan, 195–222, 2009. a, b
Schultz, T.: Implementation of a firn densification model including an optimization approach as used in Schultz et al., The Cryosphere, Zenodo [code], https://doi.org/10.5281/zenodo.5771956, 2021. a
Simonsen, S. B., Stenseng, L., Aðalsgeirsdóttir, G., Faustio, R. S., Hvidberg, S., and LucasPicher, P.: Assessing a multilayered dynamic firncompaction model for Greenland with ASIRAS radar measurement, J. Glaciol., 59, 545–558, https://doi.org/10.3189/2013JoG12J158, 2013. a
Stephenson, P. J.: Some Considerations of Snow Metamorphism in the Antarctic Ice Sheet in the Light of Ice Cycle Studies, in: Physics of Snow and Ice, vol. 1, Part 2, edited by: Ōura, H., Hokkaido University, Institute of Low Temperature Science, Sapporo, 725–740, 1967. a, b
Sturm, M., Holmgren, J., König, M., and Morris, K.: The thermal conductivity of seasonal snow, J. Glaciol., 43, 26–41, https://doi.org/10.3189/S0022143000002781, 1997. a
Székely, G. J., Rizzo, M. L., and Bakirov, N. K.: Measuring and testing dependence by correlation of distances, Ann. Stat., 35, 2769–2794, https://doi.org/10.1214/009053607000000505, 2007. a, b
Theile, T., Löwe, H., Theile, T. C., and Schneebeli, M.: Simulating creep of snow based on microstructure and the anisotropic deformation of ice, Acta Material., 59, 7104–7113, https://doi.org/10.1016/j.actamat.2011.07.065, 2011. a, b, c
Van Wessem, J. M., Reijmer, C. H., Morlighem, M., Mouginot, J., Rignot, E., Medley, B., Joughin, I., Wouter, B., Depoorter, M. A., Bambler, J. L., Lenaerts, J. T. M., Van De Berg, W. J., Van Den Broeke, M. R., and Van Meujgaard, E.: Improved representation of East Antarctic surface mass balance in a regional atmospheric climate model, J. Glaciol., 60, 761–770, https://doi.org/10.3189/2014JoG14J051, 2014. a, b, c, d, e, f, g
van Wessem, J. M., van de Berg, W. J., Noël, B. P. Y., van Meijgaard, E., Amory, C., Birnbaum, G., Jakobs, C. L., Krüger, K., Laenarts, J. T. M., Lhermitte, S., Ligtenberg, S. R. M., Medley, B., Reijmer, C. H., van Tricht, K., Trusel, L. D., van Ulft, L. H., Wouters, B., Wuite, J., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 2: Anarctica (1979–2016), The Cryosphere, 12, 1479–1498, https://doi.org/10.5194/tc1214792018, 2018. a, b
Vandecrux, B., Mottran, R., Langen, P. L., Fausto, R. S., Olesen, M., Stevens, C. M., Verjans, V., Leeson, A., Ligtenberg, S., Kuipers Munneke, P., Marchenko, s., van Pelt, W., Meyer, C. R., Simonsen, S. B., Heilig, A., Samimi, S., Marshall, S., Machguth, H., MacFerrin, M., Niwano, M., Miller, O., Voss, C. I., and Box, J. E.: The firn meltwater Retention Model Intercomparison Project (RetMIP): evaluation of nine firn models at four weather station sites on the Greenland ice sheet, The Cryosphere, 14, 3785–3810, https://doi.org/10.5194/tc1437852020, 2020. a
Vandecrux, B., Colgan, W., Solgaard, A. M., Steffensen, J. P., and Karlsson, N. B.: Firn Evoluation at Camp Century, Greenland: 1966–2100, Front. Earth Sci., 9, 1–16, https://doi.org/10.3389/feart.2021.578978, 2021. a
Wilhelms, F.: Density of ice core ngt03C93.2 from the North Greenland Traverse, PANGAEA [data set], https://doi.org/10.1594/PANGAEA.56560, 2000. a, b, c, d, e, f, g
Zwinger, T., Greve, R., Gagliardini, O., Shiraiwa, T., and Lyly, M.: A full Stokesflow thermomechanical model for firn and ice applied to the Gorshkov crater glacier, Kamchatka, Ann. Glaciol., 45, 29–37, https://doi.org/10.1029/2006JF000576, 2007. a