How Well Does Grain Boundary Sliding Represent Densification Rates in the Upper Firn?
Welcome!
To show you our work on grain boundary sliding I made this site. It is based on WebSlides. As though WebSlides provides mobile support, I recommend to switch to a bigger screen, if you are using a mobile device at the moment. Some elements, especially interactive ones, are not optimized for use on mobile devices. For the best user experience I recommend entering full screen mode. If the text still does not fit the screen, try to zoom out a little.
If you encouter any troubles using this site, please send me a mail. Briefly describe your problem and the system you are using. I'll try to fix it. Nevertheless, keep in mind that I'm simply an interested enthusiast and not a web developer. And now, just take a look around, have fun, and contact me via the chat or mail if you have any questions. If you'd like to, leave a comment.
Timm

↑ ↓ use up and down arrow keys or scroll to navigate through the slides

🖱 move your mouse to the upper border of the screen to see some contact information

🖱 move your mouse to the lower middle of the screen to see navigation symbols

2/12 click on the slide number to open an overview of all slides

links are shown in blue, try to click them
The Idea
The idea that grain boundary sliding is driving firn densification at low densities was questioned numerous times. For example by Theile et al. (2011). Therefore we wanted test if the description of grain boundary sliding is suitable to simulate firn densification at low depths by adopting the idea of Alley (1987) and reevaluate the constitutive relation using measured firn profiles.
Not only the amount of firn density profiles available for comparision has grown. Forcing data covering several decades can also help to create more detailed simulations. We use a complete model setup solving for various physical properties, including for example the temperature. The Surface Mass Balance and Snow Depth on Sea Ice Working Group (SUMup) snow density subdataset by Koenig & Montgomery (2019) provides firn profiles from Greenland and Antarctica for comparision. Forcing data for the simulation comes from the regional climate model RACMO (van Wessem et al. (2014), Noël et al. (2015)).
We use a more global form of Alley's constitutive relation described by Arthern & Wingham (1998) and formulate variants of it incorporating a single factor. This factor is then changed until an optimal match between the simulated firn profile and its corresponding measured density profile is found.
The Concept of Grain Boundary Sliding
The principle of grain boundary sliding can be understood by taking a look at a bumpy grain bundary as shown in the figure on the right adapted from Raj & Ashby (1971).
As grains are loosley packed in the first stage of firn densification they can slide against each other. Resulting in shear stress at the grain boundary, indicated with the symbol \( \tau_0 \) on the right hand side.
This shear stress leads to pressure differences along the bumpy grain boundary and finally to mass transport due to diffusion.
Alley's Approach to Grain Boundary Sliding
In an 1987 paper Alley proposed that the process of grain boundary sliding, possibly among other processes, drives the densification of firn in the upper meters of the firn column, below the critical density of 550 kg m^{3}. The paper presents a constitutive relation describing grain boundary sliding in an continuum mechanical approach.
This constitutive relation was then tested by fitting it to four different firn profiles available at the time and evaluating resulting model parameters. On the right hand side an image extract from Alley's work is shown. The figure illustrates the four firn profiles and the fitted density curves from the original paper, as well as simulation results from two other firn densification models for comparision. Recreation of the fits computed by Alley did not succeeded entirely, as can be seen in case of Site A and South Pole.

Site
Temperature
Accumulation

Ridge BD
246.65 K
0.083 m weq. a^{1}

Site A
243.65 K
0.284 m weq. a^{1}

South Pole
222.15 K
0.070 m weq. a^{1}

Dome C
218.85 K
0.034 m weq. a^{1}
Constitutive Relation
The constitutive relation that describes the process of grain boundary sliding in an approach following continuum mechanics and used in this study is shown on the right. Alley (1987) used it first for the description of firn. Arthern & Wingham (1998) developed it further to the form shown. The Arrhenius law describes the amout of boundary diffusion.
$$ \dot{\varepsilon} _{zz} = \frac{2}{15} \quad \delta_b \quad \frac{8 \, D_\mathrm{BD} \, \Omega}{k_b \, T \, h^2} \quad \frac{1}{r \, \mu^2} \quad \left( \frac{\rho _\mathrm{ice}}{\rho} \right)^3 \quad \left( 1  \frac{5}{3} \frac{\rho}{\rho_\mathrm{ice}} \right) \quad t_{zz} $$
\( D_\mathrm{BD} = A_\mathrm{DB} \exp \left( \frac{Q_\mathrm{BD}}{R \, T} \right) \)

\( \dot{\varepsilon} _{zz} \)
vertical strain rate

\( \delta_b \)
width of grain boundary

\( D_\mathrm{BD} \)
boundary diffusion

\( \Omega \)
volume of \( H_2 O \) molecule

\( A_\mathrm{BD} \)
boundary diffusion coefficient

\( k_b \)
Boltzmann's constant

\( T \)
Temperature

\( h \)
amplitude grain boundary obstructions

\( r \)
grain radius

\( Q_\mathrm{BD} \)
activation energy boundary diffusion

\( \mu \)
ratio grain radius / neck radius

\( \rho_\mathrm{ice} \)
ice density

\( \rho \)
firn density

\( t_{zz} \)
stress in vertical direction

\( R \)
gas constant
Properties Solved

Density
The evolution of the density is determined by the constitutive equation for grain boundary sliding as defined by Alley (1987) in the form of Arthern & Wingham (1998) and variants of this description. $$ \dot{\varepsilon} _{zz \, \mathrm{v_1}} = C _\mathrm{v_1} D_\mathrm{BD} \frac{1}{T} \frac{1}{3} \left( \frac{\rho_\mathrm{ice}}{\rho} \right) ^3 \left( 1  \frac{5}{3} \frac{\rho}{\rho_\mathrm{ice}} \right) t_{zz} $$ 
Grain Size
The grain size is caluclated using a simple temperature dependend evolution equation using parameters by Arthern et al. 2010.
$$ \frac{\partial r^2}{\partial t} = k_0 \exp \left( \frac{E_g}{R \, T} \right) $$ 
Age
With an Lagrangianlike description, the evolution of the age can be expressed in the following way.
$$ \frac{\partial \chi}{\partial t} = 1 $$ 
Temperature
Due to the model implementation the temperature can be solved by simple heat diffusion.
$$ \rho c_p \left( \frac{\partial T}{\partial t} \right) + \frac{\partial}{\partial z} \left( k \left( \rho \right) \frac{\partial T}{\partial z} \right) = 0 $$ 
Heat Capacity
We assume the heat capacity to be constant following Paterson (1994).
$$c_p = 2009 \, \mathrm{J \, kg}^{1} \mathrm{K}^{1}$$ 
Thermal Conductivity
The description of the thermal conductivity also follows Paterson (1994) with an additional dependency on the density. $$k (\rho) = 2.1 \left( \frac{\rho}{\rho_\mathrm{ice}} \right) \, \mathrm{W \, m}^{1} \, \mathrm{K}^{1}$$
Model Implementation
We follow the approach of a simple one dimensional model, representing a firn column. In principle the set of equations is solved in an Eulerian reference frame. Nevertheless, the grid representing the firn profile is updated at every time step. The grid point velocity \(v_b\) is set to equal the advection velocity \(v\). This scheme results in a Lagrangian like model description. This can easily be shown by integrating the local form of the mass balance in an Eulerian reference frame over a test "volume" in the zdirection (e.g Ferzinger & Perić, 2002).
$$ \frac{d}{dt} \int _{z_1 (t)} ^{z_2 (t)} \rho \, dz + \int _{z_1 (t)} ^{z_2 (t)} \frac{\partial}{\partial z} \left( \rho \left( v  v_b \right) \right) \, dz = 0 $$
When the grid point velocity equals the advection velocity \(v \equiv v_b\) the advection term drops and the equation results in the Lagrangian form of the mass balance.
$$ \frac{d}{dt} \int _{z_1 (t)} ^{z_2 (t)} \rho \, dz = 0 $$
Accumulation can easily represented by adding grid points at the top of the simulated firn profile.
Forcing
Data for the surface mass balance and the surface temperature, driving the simulation, is provided by the regional climate model RACMO (van Wessem et al., 2014; Noël et al., 2015). Data from RACMO covers the period from 1958 to 2016 in case of the Greenland Ice Sheet and from 1979 to 2016 in case of Antarctica. The figure on the right hand side shows the available forcing for the location of ice core ngt03C93.2 (Wilhelms, 2000), drilled in central Greenland in 1993.
Every simulation starts with a spin up. During this spin up we use constant average values calculated from the site specific forcing data to drive the model until a steady state is reached. The end of the spin up marks the date of the earliest available forcing, which is 1958 in case of Greenland and 1979 in case of Antarctica. Then follows the transient part of the simulation in which we use the site specific forcing from RACMO.
The surface density is determined within the optimisation scheme. Values of the surface density are assumed to be site specific and constant over time. Due to the lack of data we use a constant value of 0.5 mm for the grain radius at the firn profile's surface. Test have shown rather small influence of this parameter to the overall simulation result.
Firn Profiles
To assess the description of grain boundary sliding, we used 159 firn profiles from Greenland and Antarctica incorporated in the
Surface Mass Balance and Snow Depth on Sea Ice Working Group (SUMup) snow density subdataset
Browse and explore the data using the slider. Note that not the entire length of the profiles is shown.
map data: Amante & Eakins (2009), Arndt et al. 2013, SCAR Antarctic Digital Database.
Variants of the Constitutive Relation
The following four variants of the constitutive equation are tested.
The factors denoted by C_{v1} to C_{v4} incorporate various material constants, that might be flawed.
Nevertheless the general structure of the equation is preserved.
Variant 1
\[ \color{#bf616a} \dot{\varepsilon} _{zz \, \mathrm{v_1}} = C _\mathrm{v_1} D_\mathrm{BD} \frac{1}{T} \frac{1}{r} \left( \frac{\rho_\mathrm{ice}}{\rho} \right) ^3 \left( 1  \frac{5}{3} \frac{\rho}{\rho_\mathrm{ice}} \right) t_{zz} \]Variant 2
\[ \color{#81a1c1} \dot{\varepsilon} _{zz \, \mathrm{v_2}} = C _\mathrm{v_2} D_\mathrm{BD} \frac{1}{T} \frac{1}{r} \left( \frac{\rho_\mathrm{ice}}{\rho} \right) ^3 \left( 1 + \frac{0.5}{6}  \frac{5}{3} \frac{\rho}{\rho_\mathrm{ice}} \right) t_{zz} \]Variant 3
\[ \color{#a3be8c} \dot{\varepsilon} _{zz \, \mathrm{v_3}} = C _\mathrm{v_3} \frac{1}{T} \frac{1}{r} \left( \frac{\rho_\mathrm{ice}}{\rho} \right) ^3 \left( 1  \frac{5}{3} \frac{\rho}{\rho_\mathrm{ice}} \right) t_{zz} \]Variant 4
\[ \color{#ebcb8b} \dot{\varepsilon} _{zz \, \mathrm{v_4}} = C _\mathrm{v_4} \frac{1}{T} \frac{1}{r} \left( \frac{\rho_\mathrm{ice}}{\rho} \right) ^3 \left( 1 + \frac{0.5}{6}  \frac{5}{3} \frac{\rho}{\rho_\mathrm{ice}} \right) t_{zz} \]Optimisation Approach
For every measured firn profile we test different values for the surface density and the factor incorporated in the different variants of the constitutive relation. We then force the model using data for the specific site of the measurement. The resulting firn profile is compared to the measured one by calculating the root mean square deviation (RMSD).
Only few density profiles start right at the surface.
In case of the example shown on the right, the upper part of ice core ngt03C93.2 (Wilhelms, 2000), data is available at a depth of about 1.3 metres.
This raises the question of the correct surface density.
As our model implementation is fast, we decided to follow the simple approach of testing surface densities between
On the right hand side you can see an example for this process using
firn profile: Wilhelms (2000)
Comparision
The figure shows the density profile of ice core ngt03C93.2 (Wilhelms, 2000) as well as optimisation results for the four variants of the constitutive equation in different colours. Additionally the surfaces of 1958, 1970, 1980 and 1993 are shown, indicated by horizontal dashed lines. 1958 marks the earliest year for which forcing data is available. In 1993 the core was drilled. Whereas the horizontal lines, left of the vertical dashed line, originate from the model results, the black dashed horizontal lines show the same horizons as reconstructed from the core itself by Miller & Schwager (2004). The past surfaces show good agreement, indicating a good simulation result.
A crucial part of the optimisation scheme is the comparision between measured and modelled density profiles. To guarantee a high quality in the computation of the deviation, we restrict the simulated profiles to the domain defined by the profile surface and the horizon representing the surface of the earliest available forcing. In case of the shown example this is the surface of 1958. Everything below this surface results from the spin up and is not affected by the forcing.
In other cases the critical density of
Best Fits
The interactive plot on the right shows the best fits for all investigated firn profiles using
Horizontal lines illustrate the horizons representing the firn surface at the year of the earliest available forcing.
1958 in case of Greenland and 1979 in case of Antarctica.
If no horizontal lines appear in the plot, they are located outside the plotted area.
For example in case of
In most cases a rather good match between simulated and measured firn profiles can be found.
All profiles are clipped at a depth of
Correlation
Try which values show the most correlation yourself using the control elements. You can choose from various properties and results for the different variants of Alley's constitutive law.
We found an interesting correlation between the factors resulting in the best match between measured and modelled densities (factor) and the mean surface mass balance (mean SMB).
The value "r" in the upper left corner of the plot shows the Pearson correlation coefficient for the plotted values.
Dependency on the SMB
The linear correlation of the factors resulting in the best matches between simulated and measured firn profiles with the surface mass balance is relatively high.
This applies for all four variants of the constitutive relation.
In case of
We interpret this dependency of the optimisation result on the SMB as an incomplete description of the stress regime within the constitutive relation. Another look at the constitutive relation considering theory of continuum mechanics reveals the problem. In future models also horizontal components of the stress have to be included. A further interpretation of the resulting factors without the consideration of these components is difficult. Future work therefore has to describe the process of grain boundary sliding using an axially symmetrical three dimensional approach
Another Look at the Constitutive Relation
$$ \dot{\varepsilon} _{zz} = \frac{2}{15} \delta _b \frac{8 D_\mathrm{DB} \Omega}{k_b T h^2} \frac{1}{r \mu^2} \left( \frac{\rho_\mathrm{ice}}{\rho} \right) ^3 \left( 1  \frac{5}{3} \frac{\rho}{\rho_\mathrm{ice}} \right) t _{zz} $$A closer look at the constitutive relation by Alley (1987) in the form of Arthern & Wingham (1998) reveals, that it resembles the constitutive relation of an linearviscous but incompressible fluid. The thermodynamic pressure can be neglected, leading to the form in which strain rate and stress are connected only by the viscosity. However, firn is compressible! The density of an incompressible fluid would be constant. The computation of densification is only possible, due to the one dimensional simulation approach and neglects the additional term in the description of a linearviscous compressible fluid. This problem has been identified before and was described in detail by Greve & Blatter (2009). It may explain the observed dependency on the surface mass balance.

\( \boldsymbol{T} \)
Cauchy stress tensor

\( \boldsymbol{D} \)
spatial strain rate tensor

\( p \)
thermodynamic pressure

\( \boldsymbol{1} \)
identiy tensor

\( \rho \)
mass density

\( \lambda \)
Lamé's first parameter

\( \eta \)
shear viscosity

\( \mathrm{tr} \)
trace
Compression of an Incompressible Fluid?
Using a one dimensional simulation approach and taking a look at a one dimensional "control volume", it is rather easy to calculate the stress due to overburden firn at any point. The stress in vertical direction can then be used to compute the strain rate in the same direction. By intergrating the strain rate over time, one can determine a change in length of the "control volume". Using the mass conversation in the form dervied before, the lenght change results in a change of the density.
Taking a look at a three dimensional case, the process works basically in the same manner. But a change of the control volume's height, would result in elongation in the other directions. The density would not change! To fix this behaviour, the constitutive relation for grain boundary sliding has to be expanded so it describes a linearviscous compressible fluid. This will be the scope of further work.
linearviscous incompressible fluid
linearviscous compressible fluid
Any Questions?
Join me in the chat, write me a mail, leave a comment, or contact me in any other way.
Timm Schultz M.Sc.
Technische Universität Kaiserslautern
Chair of Applied Mechanics
PO 3049
D67653 Kaiserslautern
✉ tschultz@rhrk.unikl.de
☎ +49 631 205 3042
Further Reading
We described this work in a paper submitted to The Cryosphere, which is open for disscussion until 07 May 2021. If you like to, you can read some more about our approach to grain boundary sliding and firn densification there. We would appreciate if you take part in the discussion.
References

Alley, R. (1987).Firn Densification by GrainBoundary Sliding: A First Model. J. Phys. Colloques, 48, C1249  C1256.

Amante, C. and Eakins, B. W. (2009). ETOPO1 1 ArcMinute Global Reflief Model: Procedures, Data Sources and Analysis. NOAA Technical memorandum NESDIS NGDC24. National Geophysical Data Center, NOAA.

Arndt, J. E., Schenke, H. W., Jakobsson, M., Nitsche, F. O., Buys, O. Goleby, B., Rebesco, M., Bohoyo, F., Hong, J., Black, J., Greku, R., Udintsev, F. B., ReynosoPeralta, W., Taisei, M. and Wigley, R. (2013). The International Bathymetric Chart of the Southern Ocean (IBSCO) Version 1.0  A new bathymetric compilation covering circumAntarctic waters. Geophysical Research Letters, 40 (12), 31113117.

Arthern, R. J., Vaughan, D. G., Rankin, A. M., Mulvaney, R. and Thomas, E. R. (2010). In situ measurements of Antarctic snow compaction compared with predictions of models. J. Geophys. Res., 115 (F03011).

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. Climate Change, 40, 605624.

Bréant, C., Martinerie, P., Orsi, A., Arnaud, L. and Landais, A. (2017). Modelling firn thickness evolution during the last deglaciation: constraints on sensitivity to temperature and impurities. Clim. Past, 13, 833853.

Ferziger, J. H. and Perić, M. (2002). Computational Methods for Fluid Mechanics. Springer, Berlin, Heidelberg, New York, Barcelona, Hong Kong, Milan, Paris, Tokyo, 3., rev. edn.

Greve, R. and Blatter, H. (2009). Dynamics of Ice Sheets and Glaciers. Springer, Berlin, Heidelberg.

Haupt, P. (2002). Continuumn Mechanics and Theory of Matericals. Springer, Berlin, Heidelberg, New York. Second Edition.

Koenig, L. and Montgomery, L. (2020). Surface Mass Balance and Snow Depth on Sea Ice Working Group (SUMup) snow density subdataset, Greenland and Antarctica, 19502018. Arcitc Data Center.

Miller, H. and Schwager, M. (2004). Accumlation rate and stable oxygen isotope raios of ice core ngt03C93.2 from the North Greenlan Traverse. PANGAEA.

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. (2015). Evaluation of the updated regional climate model RACMO2.3: summer snowfall impact on the Greenland Ice Sheet. The Cryosphere, 9 (5), 18311844.

Paterson, W. S. B. (1994). The Phyiscs of Glaciers. Third Edition, Pergamon.

Raj, R. and Asby, M. F. On grain boundary sliding and diffusional creep. Metall Mater Trans, B2, 11131127.

Theile, T., Löwe, T. C., Theile, T. C. and Scheebeli, M. (2011). Simulating creep of snow based on microstructure and the anisotropic deformation of ice. Acta Materialia, 59 (18), 71047113.

van Wessem, J. M., Reijmer, C. H., Morlighem, M., Mouginot, J., Rignot, E., Medley, B., Joughin, I., Wouters, B., Depoorter, M. A., Bamber, J. L., Lenaerts, J. T. M., can de Berg, W. J., van den Broeke, M. R. and van Meijgaard, E. (2014). Improved representation of East Antarctic surface mass balance in a regional atmospheric climate model. Journal of Glaciology, 60 (222), 761770.

Wilhelms, F. (2000). Density of ice core ngt03C93.2 from the North Greenland Traverse. Pangaea.
Acknowledgements
This work has been funded by the German Reasearch Foundation (project number: 403642112), it is part of the priority program 1158 "Antarctic Research with Comparative Investigations in Arcitc Ice Areas".