On the Contribution of Grain Boundary Sliding to Firn Densification-an Assessment using an Optimisation Approach

Physics based simulation approaches to firn densification often rely on the assumption that grain boundary sliding, first introduced by Alley (1987) to firn, is the leading process driving the first stage of densification. However, often so called semi empirical models are favored against the description of grain boundary sliding due to simplicity and uncertainties regarding model parameters. In this study, we are assessing the applicability of grain boundary sliding to firn using a numeric firn densification model and an optimisation approach, for which we formulate variants of the constitutive relation by Alley (1987). 5 The efficient model implementation based on an updated Lagrangian numerical scheme enables us to perform a large number of simulations testing different model parameters, to find simulation results suiting 159 firn density profiles from Greenland and Antarctica best. For most of the investigated locations a good agreement of simulated and measured firn density profiles was found. This implies that the constitutive relation by Alley (1987) characterises the fist stage of firn densification well, if suitable model parameter are used. An analysis of the parameters that lead to best matches reveals a dependency on the mean 10 surface mass balance. This may indicate an insufficient description of the load situation, as horizontal components of the stress tensor are usually neglected in one dimensional models of the firn column.

The second category of firn densification models tries to quantify the physical processes related to the densification of firn. These processes incorporate different types of creep and diffusion. Micro mechanical models are used for small scale investigations (Johnson and Hopkins, 2005;Theile et al., 2011;Fourtenau et al., 2020) while continuum mechanics based models can be used for large scale simulations. Examples of the latter are models by Wingham (1998), Arnaud et al. (2000) and Goujon et al. (2003). These models again are based on the concepts first introduced by Maeno and Ebinuma 30 (1983) who adopted modelling approaches from sintering theory, developed in the field of metallurgy, to firn.
The theory of hot isostatic pressing, which applies for firn (Maeno and Ebinuma, 1983), is formulated for well compacted powders as used for many purposes in industry (eg. Kang, 2005). However, the densification of firn starts with far less regular grain geometries and conglomerates of grains. Therefore the first stage of firn densification, namely up to the critical density of ρ c = 550 kg m −3 (Herron and Langway, 1980), has to be described by other processes. 35 Alley (1987) first applied grain boundary sliding, adopted from Raj and Ashby (1971), to firn densification at low density.
Since then the description of this process was used in various other firn densification models (Arthern and Wingham, 1998;Arnaud et al., 2000;Goujon et al., 2003;Bréant et al., 2017). Nevertheless, the assumption that grain boundary sliding is the dominant process in firn densification at densities below ρ c = 550 kg m −3 was questioned numerous times. For example Theile et al. (2011), by conducting experiments on a small number of snow samples, tried to point out that the densification is more 40 likely driven by processes within the grain than by the inter granular process of grain boundary sliding.
When Alley (1987) first published the description of grain boundary sliding for firn, he tested the material model by fitting the model results to four firn profiles available at this time, evaluating the resulting model parameters. As indicated by the title of his paper and pointed out in its discussion, grain boundary sliding might not be the only process driving the densification of firn at low density and model parameters might differ, but by using the given constitutive law it is possible to recreate measured 45 depth density profiles to a satisfying degree.
In this study we aim at evaluating whether grain boundary sliding and its description is suitable for the simulation of firn densification at low density. In the sense of Alley (1987) we compare model results to measured data. In contrast to Alley (1987) the amount of available data became much larger. Alongside to a large number of firn profiles this includes forcing data which, together with additional modelling techniques, allows us to simulate firn profiles at a very high quality. In order to test the description of grain boundary sliding by Alley (1987) we use a numeric model, simulating the evolution of a one dimensional firn column with respect to time. The model incorporates variants of the constitutive relation of Alley (1987), all of which combine several model parameters in a single factor. We then force the model with data provided by the concept of an updated Lagrangian description with the update velocity being the flow velocity. This results in material fixed coordinates. The Lagrangian like description allows for a very high spatial and temporal resolution in the simulations. It can be shown that 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) leads to (Ferziger and Perić, 2002):

70
Here ρ describes the density, z the vertical coordinate, t the time and v is the flow velocity, whereas v b represents the velocity of the grid or integration boundary. When the boundary velocity equals the material velocity v b = v, the second part of Equation (1), describing advection, vanishes. The resulting equation is equal to the Lagrangian form of the mass balance (Ferziger and Perić, 2002). On a one dimensional grid, build up by a number of grid points, as illustrated in Fig. 1, we define the grid point velocity to be v b and to equal the flow velocity v b ≡ v. The location of all grid points is updated in every time step by 75 integrating the grid point velocity using a forward Euler scheme. In this way advection is entirely represented by the adapting grid.
The grid point velocity is calculated using the constitutive equation of grain boundary sliding as described in Section 3.2 and the definition of the strain rate in one dimension. The description of grain boundary sliding provides the strain rate in vertical direction along the gridε zz as a function of the vertical stress t 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 Equation (2) (Haupt, 2002). On a one dimensional grid, defined by a number of grid points, the space between neighbouring grid points can be considered as a material line element (see Fig. 1). Therefore the grid point velocity v b can easily be computed by integrating the strain rateε zz in vertical direction along the length of the grid cell: For the implementation of Equation (3) an integration constant determined by a suitable boundary condition is needed. It is reasonable to apply a Dirichlet boundary condition, forcing the grid point velocity to be zero, at either the top or the base of the computational domain, representing a fixed reference point at the top or the base of the modelled firn profile, respectively. All other points defining the adapting grid are moving with respect to this anchor point. In case of the present study we decided to 90 place the anchor point at the base of the simulated firn profiles (z 0 in Fig. 1). Depth coordinates of profiles shown in following figures were adjusted for better readability.
For the representation of accumulation at the top of a simulated firn profile an inflow boundary condition has to be implemented. To achieve this, in every time step an additional grid point is added at the top of the grid. Its coordinate within the grid z n+1 (t + ∆t) is calculated by z n+1 (t + ∆t) = z n (t + ∆t) + ∆t a 0 (t) ρ water ρ 0 , with a 0 the accumulation rate given in m s −1 water equivalent, ∆t the length of the time step, ρ water the density of water and ρ 0 the density of the deposited snow. The position of a new grid point (z n+1 (t + ∆t) in Fig. 1) 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 time step ∆t and the time dependent accumulation rate a 0 (t). As we use an accumulation rate given in the unit of 100 meter water equivalent per second, it has to be converted to the unit of meter firn equivalent taking the site specific surface firn density ρ 0 and water density ρ water into account.

Grain Boundary Sliding
Equation (5) shows the constitutive law by Alley (1987), describing the process of grain boundary sliding. In the following the different components and characteristics of the equation will be explained briefly.
The factor of 2/15, seeming arbitrary at first, results from the geometric deviation pointed out by Alley (1987). Another geometric parameter δ b describes the width of the grain boundary. The following term describes the reciprocal bond viscosity (Raj and Ashby, 1971) where Ω stands for the volume of the H 2 O molecule, k b is the Boltzmann constant, T resembles the temperature and h describes the amplitude of the grain boundary obstructions, a measure for the roughness of the grain 110 boundaries. The coefficient D BD is given by an Arrhenius law describing the rate of boundary diffusion depending on the temperature. Corresponding values for the activation energy and a prefactor have to be provided.
The strain-rate resulting from grain boundary sliding also depends on the grain radius r. The ratio of grain radius to neck radius µ was introduced by Arthern and Wingham (1998) and is assumed to be constant. The next factor of Equation (5) describes the dependency on the inverse relative density to the power of three. The characteristics of the depth density profiles 115 resulting from the model description are greatly influenced by this and the following term of the constitutive law. The factor of 5/3 corresponds to the inverse relative density of ρ c = 550 kg m −3 , resulting in a kind of fade out behaviour ending with a strain rate of zero when reaching this critical density.
Finally the stress in vertical direction t zz resulting from the overburden firn is driving grain boundary sliding. The other physical properties influencing the process are density ρ, temperature T , and grain size r.

Stress
For the evaluation of the stress in vertical direction t zz we use the local form of the static linear momentum balance in its Eulerian description. We neglect the acceleration term, as changes in velocity are assumed to be small, leading to ∂t zz ∂z + ρ g = 0 .
Computation of the stress t zz can easily be achieved by integrating the product of density ρ and acceleration due to gravity g 125 along the simulated profile. We assume the surface of the profile to be traction free.

Density
As pointed out in Section 3.1 and illustrated by Equation (1) the change of the density integrated over a test volume with respect to time has to be zero. Or in other words, the mass incorporated in a material control volume cannot change. As the position of the grid points, and therefore the material control volume, does change, the density is changing accordingly. The evaluation of 130 the density integral over a control volume at two time steps leads to As the space between two neighbouring grid points can be understood as a material line element (Haupt, 2002), Equation (7) can be rewritten in the form of Equation (8). Where |dZ| = |z 2 (t) − z 1 (t)| and |dz| = |z 2 (t + ∆t) − z 1 (t + ∆t)| are the lengths of a material line element in the reference configuration and its image in the current configuration respectively: Sorting Equation (8) leads to the formulation of the density change with respect to time depending on the definition of the strain ε zz in one dimension (Haupt, 2002) The evolution of density can therefore be computed by integration of the strain rateε zz (Section 3.2) in time.

Temperature
As pointed out in Section 3.1 all advection in the model domain is represented by the moving grid the model is discretized by.
Therefore the description of 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 a density dependent thermal 145 conductivity, described by

Grain Size
To describe the evolution of grain size we use the well known description of Stephenson (1967) and Gow (1969) as given in Arthern et al. (2010) for mean grain radius r rather than for mean grain cross-sectional area. This is suitable as the grain radius 150 is used in the constitutive equation by Alley (1987). It is written in differential form by The square of the grain radius increases with a rate in the form of an Arrhenius law, where R is the gas constant and T represents the absolute temperature. Values for the activation energy E g = 42.4 kJ mol −1 of the process and the corresponding pre-factor k 0 = 1.3 × 10 −7 m 2 s −1 were adapted from Arthern et al. (2010) following Paterson (1994). In contrast to Arthern 155 et al. (2010) we do not use the mean annual temperature but the actual temperature T (z, t) along the simulated profile.

Age
For reasons of comparison (Section 5.2) and general interest additionally the firn age χ is simulated. Again due to the fact that advection is represented by the adapting grid, the description is very simple. It is calculated from

160
New deposited snow has an age of zero, prescribed in the form of a Dirichlet boundary condition. The age discretized at the grid points then increases according to the time step.

Time Discretisation
Time is discretised using constant time steps. For this study 48 time steps per year have shown to be a good compromise between simulation costs and resolution and was used throughout all simulations. The grid resolution depends on the time 165 step as a new grid point is generated in every time step representing accumulation as described in Section 3.1 and shown in Equation (4). Time dependent properties such as the density, temperature, grain size and age were developed using an explicit Euler scheme.

Optimisation
To test the concept of the material model developed by Alley (1987) we vary Equation (5) (14) to (17)) preserve its general form, but group several material parameters into a single factor. Arnaud et al. (2000), Goujon et al. (2003) as well as Bréant et al. (2017) also summed the material parameters of the model by Alley (1987) up into a single parameter . In the study by Bréant et al. (2017) additionally the factor of 5/3 was modified to change the fade out behaviour (see Section 3.2) of the constitutive equation. In the following the four variants, indicated by the subscripts (·) v1 175 to (·) v4 , are shown: Variant 1 (Equation (14)) and Variant 2 (Equation (15)) of the constitutive equation combine all material constants using the factors C v1 and C v2 , respectively. Excluded from the factor is the Arrhenius law for D BD with its activation energy Q and the corresponding pre-factor A. To preserve the additional dependency on the temperature T constant values of A = 3.0 × 10 −2 m 2 s −1 and Q = 44.1 kJ mol −1 are used (Alley, 1987;Paterson, 1994;Arthern and Wingham, 1998). Except for the temperature T , the vertical strain rateε zz only depends on the firn density ρ, the grain radius r and the stress in vertical direction t zz . Variant 2 differs from Variant 1 by the use of the modified fade out factor introduced by Bréant et al.

185
(2017), leading to a theoretical fade out of grain boundary sliding at a density of ρ * c = 596 kg m −3 . To test its influence, Variants 3 and 4, given in Equations (16) and (17), disregard the Arrhenius law: Aim of the optimisation is to find optimal values of the factors C vi for every variant of the constitutive relation (Equations (14) to (17)) resulting in a simulated density profile that resembles the measured profiles in the best possible way. As an example, we explain the optimisation process for one selected firn core in more detail. The upper part of ice core ngt03C93.2
Every simulation starts with a spin up in which constant values are used for the forcing. When steady state is reached, a transient run using evolving forcing data follows. The forcing at the location of ngt03C93.2 is shown in Fig. 2 (c). The resulting firn profile is then compared to the measured profile. Unlike other recent studies (Lundin et al., 2017;Verjans et al., 2020) we use the objective measure of the root mean square deviation between measured and modelled density, which allows  (14) and (16)).

210
As the implementation of our model is efficient and the approach is as well simple and reliable, we decided to determine the best factor C vi for the four variants of the constitutive equation, by simply testing 250 values within certain ranges. These  (18) and (19). They differ because of the inclusion or disregard of Arrhenius law, respectively.
The range boundaries were determined by performing several tests. The corresponding density profiles are shown in Fig. 2 (a).

220
As can be seen in Fig. 2  is available online and merges more than three thousand density measurements covering the period from 1950 to 2018 from about 180 different sources. To obtain firn profiles relevant for this study from the dataset, we filter it based on the following conditions: 1. Profiles have to consist of at least ten data points.
2. The overall length of the profiles has to exceed three meters. 3. Profiles have to start at a depth of less than three meters below surface.
4. The starting density of the profiles must not to exceed ρ c = 550 kg m −3 .
5. The surface mass balance at the profile locations has to be positive.
6. Forcing data of at least five years have to be available for the profile location.

Boundary Conditions and Forcing
For the forcing of the firn densification model values for surface density, surface temperature, accumulation rate and surface grain size at the locations of the measured firn profiles (Section 5.1) are needed. We use the regional climate model RACMO2.3 (Van Wessem et al., 2014;Noël et al., 2015) for this purpose. In case of Greenland RACMO2.3 provides annual data from 1958 to 2016. At every location of a measured depth density profile we interpolate the mean annual skin temperature and surface better comparison options the grain size at the surface is assumed to be the same at every location and to be constant over time.
We chose to use a grain radius of r 0 = 0.5 mm. with characteristic values of the surface density obtained and used in other studies (Sugiyama et al., 2012;Lundin et al., 2017;275 Tian et al., 2018;Weinhart et al., 2020). In comparison locations in Greenland show a higher mean surface temperature and surface mass balance than locations in Antarctica. This results in lower values of the surface density in Greenland, which is plausible. In general a wide variety of typical climatic conditions for both ice sheets is covered. In order to solve for the temperature a Neumann boundary condition is used at the profile base. The first derivative of the temperature is forced to be zero.  firn profiles from ice core ngt03C93.2 (Wilhelms, 2000) displayed in Fig. 2 is about 28 kg m −3 . That means, more than half of the resulting simulations show even better agreement with the corresponding firn density profiles than this one.

290
Furthermore, the example of ice core ngt03C93.2 shows not only a good match of simulated and measured density, but also quite a good agreement to the measured age structure as illustrated by four different horizons in Fig. 2 (a). Similar results were obtained from a subsample of 40 firn profiles of which the age structure is available.
The range of factors resulting in the best fitting firn profile is smaller for Variants 2 and 4 of the constitutive equation, which use the modification by Bréant et al. (2017), compared to Variants 1 and 3, as shown in form of box plots in Fig. 6. In contrast to 295 factors C v1 and C v2 , factors C v3 and C v4 incorporate the Arrhenius law from the original description of grain boundary sliding by Alley (1987). Therefore a direct comparison between the two groups of factors is not possible. The quartile coefficient of dispersion, shown for the four variants on the right side of Fig. 6, is a relative measure for the scatter of the values. The coefficient reveals that the resulting factors of Variants 1 and 2 are defined in a range narrower versus Variants 3 and 4. All four resulting sets of factors show a slightly non uniform distribution, tending towards smaller values.    (14) to (17)), derived using the optimisation scheme described in Section 4. On the right side the quartile coefficient of dispersion vr is shown for the variants in corresponding colors. The quartile coefficient of dispersion is calculated using the first Q1 and third Q3 quartile values of the data sets as shown in black colour and represents a robust relative measure of dispersion.  to Cv 4 and the mean surface temperature, is given.
To find a possible dependency of the resulting factors, leading to the best match between simulated and measured density profiles, Fig. 7 shows the 159 factors found by the optimisation, plotted against the mean surface temperature calculated from the forcing data specific for each firn profile's site. The left plot illustrates the results for Variants 1 and 2 in red and blue colour, while the right part of the figure shows the resulting factors for Variants 3 and 4 using green and yellow markers respectively.
Lines represent a linear fit to the four scatter plots, the legend features the Pearson correlation coefficient, a measure for the 305 linear correlation of the two variables. The linear correlation of the resulting factors with the mean surface temperature is quite high compared to other properties. This is especially true for factors C v3 and C v4 of Variants 3 and 4. However the Pearson correlation coefficient only indicate the linear correlation. The scatter of factors C v3 and C v4 with respect to the mean surface temperature might resemble a higher order function.
Even higher values of the Pearson correlation coefficient can be found in Fig. 8, which shows the resulting factors in the 310 same manner as Fig. 7 but with respect to the mean surface mass balance. Just as the mean surface temperature, the mean surface mass balance is calculated from the forcing data used during the simulations. The correlation coefficients for the results of Variants 3 and 4 exceed the ones for the variants of the constitutive equation incorporating the Arrhenius law explicitly. The best indication of a linear relationship between the factor resulting in the density profile best matching the corresponding field measurement and the surface mass balance can be seen for Variant 3.

Discussion
Given the uncertainties in forcing, the limited knowledge on initial grain size and the poor constraint on density at the surface, we find a reasonable good match between grain boundary sliding based simulated density profiles and observed ones. If the goal of a modelling approach is to simulate the upper firn layer down to a density of ρ c = 550 kg m −3 well, the description of grain boundary sliding first introduced to firn by Alley (1987) is a good basis for a physically based model. We were able to find 320 distinct minima of the deviation between simulated and measured density profiles using the described optimisation scheme.
Furthermore, the simulated firn age is in good agreement with the measurement. This indicates that not only the final state of the firn profile at the time of the measurement is well represented, but also its history. A factor often neglected in firn model development before. Unfortunately the measured age structure is only available for a rather small subset of the firn density profiles used in this study.

325
However, it has to be mentioned that we have not investigated whether grain boundary sliding is indeed the dominating process during the first stage of firn densification. We only assess whether a process with a functional dependence on density, firn overburden stress, temperature and grain size, is representing observed density profiles well. Any other deformation process with the same functional dependence would be equally well suited.
Comparing the results for Variants 1 and 2, as well as Variants 3 and 4, we find that the adjustment of 1 − 5 3 ρ/ρ ice 330 to 1 + 0.5 6 − 5 3 ρ/ρ ice is leading to better matches with measured density profiles. With regard to the study design and the background of a physics based model describing firn densification, this result has to be reviewed carefully. As Alley (1987) points out, grain boundary sliding might be the dominant process driving the firn densification at low densities, but presumably not the only one. The constitutive law by Alley (1987) is designed in a way that the densification due to grain boundary sliding becomes zero at the density of ρ c = 550 kg m −3 , motivated by densest packing of spheres and increasing accommodation 335 incompatibilities. The modification by Bréant et al. (2017) changes this behaviour in the way that the process grain boundary sliding vanishes at a density of about ρ * c = 596 kg m −3 , which could have advantages for the transition into the next stage of firn densification. We suggest a simultaneous decline of grain boundary sliding and increase of one or more other processes would provide a good characterisation. Namely dislocation creep drives the densification at higher density (Maeno and Ebinuma, 1983) due to increasing stresses. An onset of dislocation creep at densities lower than ρ c = 550 kg m −3 affecting not necessarily 340 the entire bulk firn matrix, but increasing volume fractions of the porous matrix should be investigated further in future.
As pointed out in Section 6 and shown in Figure 7, a considerable correlation of the resulting factors C v3 and C v4 with the mean surface temperature exists. The same correlation is not as high for Variants 1 and 2, in which the Arrhenius law from the description of the bond viscosity is incorporated. Also the quartile coefficient of dispersion describing the scattering of factors and shown in Fig. 6 is lower for Variants 1 and 2 compared to results of Variants 3 and 4. The additional dependency on 345 the temperature in form of the Arrhenius law is captured within the factors C v3 and C v4 , but results in a worse determination of these parameters. A good determination of the activation energy for boundary diffusion and the corresponding prefactor is therefore crucial for the description of grain boundary sliding.
As the mean surface temperature correlates well with the mean surface mass balance, a dependency of the resulting factors on the mean surface mass balance, as shown in Fig. 8, seems obvious. But the quite strong correlation of results from Variants 1 350 and 2, incorporating the additional temperature dependency in form of the Arrhenius law, and higher values of the Pearson correlation coefficient indicate a clear dependency of the resulting factors on the mean surface mass balance. We interpret the dependency of the mismatch on the surface mass balance such, that the load situation is currently not represented well.
Due to the lateral confinement we expect horizontal stresses to play an increasing role with depth. The description of the load situation can only be improved solving an axially symmetrical three dimensional problem and considering further processes 355 contributing to densification. Further interpretation of the factors resulting in the best match between simulated and measured density profiles is difficult because there exists no direct dependency on the surface mass balance in the described model.

Conclusions
Using variants of the constitutive relation for grain boundary sliding by Alley (1987) and a efficient optimisation scheme, we were able to reproduce a sample of 159 firn density profiles reasonably well. Even though available only for a rather 360 small subsample of the analysed firn profiles, the comparison of the measured age structure to the simulation also shows good agreement. Thus we conclude, that the description of grain boundary sliding as introduced by Alley (1987) is suitable for the simulation of firn densification at low density.
In our optimisation approach we use a single factor representing various model parameters and search for the factor value leading to the best match between simulated and measured firn profiles. In this way the simulation result is independent of 365 the possibly deficient, now substituted, model parameters. We find a linear dependency of the factors on the surface mass