photo credit: Sepp Kipfstuhl

powered by WebSlides


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.


  • ↑ ↓ 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



  • 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 \)


  • \( 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

Alley (1987), Arthern & Wingham (1998)

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 Lagrangian-like 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 z-direction (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.



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

Koenig & Montgomery (2019).

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 Cv1 to Cv4 incorporate various material constants, that might be flawed. Nevertheless the general structure of the equation is preserved. Variant 1 and variant 2 include the Arrhenius equation DBD describing the amount of boundary diffusion, while in case of variant 3 and variant 4 this term is described implicitly within the factor. Variant 2 and variant 4 differ from variant 1 and variant 3 by the use of an addition term 0.5/6 introduced by Bréant et al. (2017). This leads to a different behaviour of the constitutive relation at densities near the ciritical density of 550 kg m-3.

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 250 kg m-3 and 450 kg m-3. We save the combination of surface density and factor leading to the best match between simulation and measurement for further analysis.

On the right hand side you can see an example for this process using variant 1 of the constitutive relation. You can try to recreate the best match, shown in red, by using the control elements. Note that this example does not picture the entire range of values tested in the study.

firn profile: Wilhelms (2000)


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 550 kg m-3 is reached above the horizon representing the earliest available forcing. In these cases the domain for the computation of the root mean square deviation between simulated and measured density is restricted to the domain showing density values below 540 kg m-3. We decided to use a smaller density threshold because of the asymptotic character of the density profiles, resulting from the used constitutive relation.

Best Fits

The interactive plot on the right shows the best fits for all investigated firn profiles using variant 1 and variant 3 of the constitutive relation. Results for variant 2 and variant 4 look very similar, the resulting factors however are different due to the different formulations.

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 profile 63, the former surfaces are located at a depth of about -30 m. You can find them by dragging the plot.

In most cases a rather good match between simulated and measured firn profiles can be found. All profiles are clipped at a depth of -25 m and/or a density of 550 kg m-3, due to issues regarding the size of the plot.


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 variant 3 and variant 4 the scattering of values might resemble a higher order function. This is not surprising as the Arrhenius law describing the amount of boundary diffusion is included within the factor for these variants. Nevertheless, the linear correleation with the SMB is higher than that with the mean surface temperature.

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 linear-viscous 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 linear-viscous 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} = -p \left( \boldsymbol{x}, t \right) \boldsymbol{1} +2 \eta \boldsymbol{D} $$
constitutive relation of a linear-viscous incompressible fluid
$$ \boldsymbol{T} = - p \left( \boldsymbol{x}, t \right) \boldsymbol{1} + 2 \eta \boldsymbol{D} \color{#bf616a} + \lambda \left( \mathrm{tr} \boldsymbol{D} \right) \boldsymbol{1} $$
constitutive relation of a linear-viscous compressible fluid
  • \( \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} \)


Haupt (2002)

Compression of an Incompressible Fluid?

mass conversation: $$ \frac{d}{dt} \int _{z_1 (t)} ^{z_2 (t)} \rho \, dz = 0 $$

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 linear-viscous compressible fluid. This will be the scope of further work.


linear-viscous incompressible fluid


linear-viscous 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
D-67653 Kaiserslautern
☎ +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.



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".

photo credit: Sepp Kipfstuhl

powered by WebSlides