Glacier Flow

This example shows how you can set up a glacier flow model in principle. It exemplifies several modeling steps which are important for glacier modeling:

•

The creation of the geometry

•

The modeling of a non-Newtonian fluid

•

A way to investigate and compare different aspects of the model within the same model file

The cryosphere is the part of the climate system that contains frozen water and makes up 80% of our fresh water. Using the COMSOL Multiphysics® software, we can simulate classical ice flow to analyze cryosphere dynamics and assess climate change effects such as shrinking glaciers and rising sea levels.

In 1773, André Bordier, a Swiss naturalist, used the term “fluid” to describe the movement of mountain glaciers for the first time. However, it took more than a century for scientists to agree on a unified description for the dynamics of glaciers.

One of the most confusing aspects of glaciers is the observation that ice exhibits both viscous and plastic behavior, depending on the glacier. British physicist John Glen observed and described this intermediate behavior using a nonlinear relationship between stress and strain. Known as shear thinning, this classical behavior applies to many different fluids (for example, ketchup and blood).

The life of any mountain glacier can be schematically described as follows:

•

In the accumulation zone, snow piles up at a high altitude where the temperature is low, and compresses into ice

•

The ice starts deforming and flowing down the slope under its own weight

•

In the ablation zone at lower altitude, where the temperature is higher, the ice melts away

Figure 1 shows a sketch of a typical mountain glacier with the accumulation zone where the snow piles up and the ablation zone where the ice melts to water.

Glaciologists devide glaciers in different categories, depending on their thermal structure:

•

Cold glaciers where the temperature of the ice is below the pressure melting temperature throughout the glacier, except for maybe a thin surface layer.

•

Temperate glaciers where the whole glacier is at the pressure melting temperature, except for some seasonal variations and freezing at the surface layer.

•

Polythermal glaciers, where some parts of the glacier are cold and some are temperate. The cold parts are usually the highest accumulation area as well as the upper part of an ice column, whereas the surface and the base are temperate. In polythermal glaciers, cold and temperate ice is separated by the cold-temperate transition surface (CTS).

Figure 1: Sketch of a typical mountain glacier.

Most of the alpine glaciers are temperate, but at higher altitudes, also polythermal glaciers can be found. Cold glaciers can typically be found in polar regions like Greenland and Antarctica. Polythermal glaciers are common at high latitudes and high altitudes (like the Scandinavian Mountains, European Alps, Himalaya).

In this example a realistic geometry - a 2D cross section inspired by the Arolla glacier in the Pennine alps in Switzerland - is used to simulate the nonisothermal flow of the ice mass downslope, under its own weight and subject to basal sliding. Two versions of the glacier model are investigated: A cold glacier version and a temperate glacier version.

Model Definition

In this example the geometry is inspired by the Haut glacier d’Arolla in the Swiss Alps, which is 5 kilometers long and up to 200 meters thick with an average slope of 15%. Figure 2 shows a photograph of the Arolla glacier.

Figure 2: Aerial photo of the Arolla glacier in the Pennine Alps in Switzerland.

The model is built in 2D as a cross section through the flow line of the glacier. It is constructed from two datasets, one for the ice surface and one for the bedrock surface. Figure 3 shows the model geometry, colored by the ice thickness.

Figure 3: A 2D cross section of the Arolla glacier is constructed as model geometry. Here it is colored by the ice thickness and the y-axis is stretched compared to the x-axis to have a better view of the vertical differences.

The ice flow is simulated using the Stokes equations describing the so-called creeping flow:

(1)

(2)

where ρ is the density, u is the velocity vector, p is the pressure, g is the acceleration of gravity, and μ is the dynamic viscosity.

Using Glen’s flow law, the viscosity μ can be described as

(3)

where

is the shear rate and is classically defined as the norm of the strain rate tensor

(4)

To calculate the flow rate factor A in Equation 3 an Arrhenius law can be used:

(5)

with A0 being the flow rate constant, Q the activation energy, and R the universal gas constant. According to the literature, A0 and Q are a matter of debate. In this case we use,

(6)

and

(7)

To calculate the temperature distribution in the model, the energy balance equation is solved:

(8)

where cp is the heat capacity and k is the thermal conductivity of ice. Both are functions of temperature and are described using the following numerical value equations (see, for example, Ref. 4):

(9)

and

(10)

In terms of fluid, the inflow and outflow boundary conditions are the normal constraints, corresponding to the applied pressure of the ice, which is not included in the domain. It simply corresponds to an assigned hydrostatic (or cryostatic) pressure. The upstream boundary weighs on the domain, thus contributing to a streamwise velocity, while the downstream boundary resists to the flow. The surface of the glacier is a free surface.

If the temperature at the ice-bed contact is at the pressure melting temperature Tm — which is the case if the glacier is temperate — the glacier can slide over the base. If the temperature is below Tm a no slip condition would be appropriate. In case of a temperate glacier, a slip length of 50 m is defined at the lower boundary.

In terms of heat transfer, the surface is influenced by the ambient temperature and is subject of convective heat exchange and radiative heat exchange with the environment. For a cold glacier, the boundary in contact with the bedrock is subject to a geothermal heat flux, which could be modeled as a boundary condition. Typical geothermal heat fluxes are of the order of 40–120 mW/m2. In this case a geothermal heat flux of 120 mW/m2 is chosen.

In case the glacier is temperate, the temperature at the bedrock can be assumed to be at pressure melting point temperature, which is

(11)

with Ttp and ptp being the triple point temperature and pressure of water (see Ref. 2). Heat is allowed to leave and enter the domain at the inflow and outflow boundaries.

A mapped mesh is used that is consistent with the aspect ratio of the geometry.

The external weather conditions are an important input data for geophysical simulations. Accessing the ASHRAE 2017 database, we can import the average external temperature and wind velocities at a given time of the year for more than 6000 weather stations all over the world. Here, we use the data from the Grand Saint Bernard station in the Swiss Alps, located at about 25 km of the Arolla glacier on the first of February at noon. The ambient temperature is imposed at the glacier surface and the wind velocities are used to simulate a convective heat flux at the surface.

Results and Discussion

The results are first evaluated after 10 years, when the flow has reached a stationary state. Figure 4 shows the velocity for the “cold glacier” conditions. At the bedrock the velocity is zero and there is hardly any in- and outflow at the upper and lower boundaries. However, there is some uncertainty regarding basal sliding which is not present in this case: Depending on the ice thickness and the geothermal heat flux, even for cold glacier conditions it is possible that a water layer builds up at the bedrock as the heat can only be transported upwards very slowly.

Figure 4: Velocity magnitude after 10 years of simulated time for a given heat flux at the bottom (cold glacier model).

The temperature after 10 years of simulation time (Figure 5) shows that the main heat source is the external temperature being in exchange with the glacier surface. Here the highest temperatures occur during summer and the lowest during winter season.The ablation zone of the glacier is colder on average than the accumulation zone. The temperature at the bedrock is higher than in the middle of the glacier due to the geothermal heat flux, however, it is still well below pressure melting temperature.

Figure 5: Temperature distribution after 10 years of simulated time for a given heat flux at the bottom (cold glacier model).

The next figures show the results for the temperate glacier conditions. Here the temperature at the bedrock is set to pressure melting point temperature and therefore it can be assumed that the ice can slide over the bedrock. The velocity field which is displayed in Figure 6 shows this as the velocity is non-zero at the bottom of the glacier. In general the velocities are higher in this case.

Figure 6: Velocity magnitude after 10 years of simulated time for the bottom temperature being at pressure melting point (temperate glacier model).

The temperature distribution after 10 years of simulation time is shown in Figure 7. Again the ablation zone is colder on average than the accumulation zone and the temperature field shows a line structure which is due to a combination of the seasonal heating and cooling from the surface and the influence of the mesh and time step size. With refined mesh the structures still appear on a smaller scale.

Figure 7: Temperature distribution after 10 years of simulated time for the bottom temperature being at pressure melting point temperature (temperate glacier model).

Figure 8 shows the tangential velocity component along the glacier surface for both the cold and the temperate glacier scenario. As already mentioned the velocities are much bigger for the temperate scenario.

Figure 8: Tangential velocity component along the glacier surface after 10 years of simulated time.

References

1. A. Aschwanden and H. Blatter “Mathematical modeling and numerical simulation of polythermal glaciers,” J. Geophys. Res., vol. 114 (F1), p. F01027, 2008. (10.1029/2008JF001028.)

2. A. Aschwanden, “Thermodynamics of Glaciers”, McCarthy Summer School, 2010. (https://glaciers.gi.alaska.edu/sites/default/files/mccarthy/Notes_thermodyn_Aschwanden.pdf)

3. R. Greve, and H. Blatter, Dynamics of Ice Sheets and Glaciers, Springer, 2009.

4. C. Ritz, “Time dependent boundary conditions for calculation of temperature fields in ice sheets,” The Physical Basis of Ice Sheet Modelling, vol. 170, pp. 207–216, 1987.

Subsurface_Flow_Module/Heat_Transfer/glacier_flow_2d

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Global Definitions

Parameters 1

In the window, under click .

Start by setting the parameters defining the material properties of ice and the constants needed to calculate the ice flow. You can either import them from the file "glacier_flow_2d_parameters.txt" from the model’s Application Librabries folder or enter them as follows:

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
rho_ice	910[kg/m^3]	910 kg/m³	Density of ice
mu_ice	5e12[Pa*s]	5E12 Pa·s	Dynamic viscosity of ice
LSlip	50[m]	50 m	Slip length
T_init	-10[degC]	263.15 K	Initial temperature
betaCC	9.8e-8[K/Pa]	9.8E-8 m·s²·K/kg	Clapeyron constant
n_ice	3	3	Rheological exponent for ice
T_tp	0.01[degC]	273.16 K	Temperature at triple point of water
p_tp	611.657[Pa]	611.66 Pa	Pressure at triple point of water
q_geo	120[mW/m^2]	0.12 W/m²	Geothermal heat flux

Glacier Base

Now import the data to specify the glacier geometry and define the functions to calculate the flow viscosity and the pressure melting temperature.

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click

Browse to the model’s Application Libraries folder and double-click the file glacier_flow_2d_arolla01.txt.

Click

In the text field, type Glacier Base.

Locate the section. In the table, enter the following settings:


Function	Unit
int1	m

In the table, enter the following settings:


Argument	Unit
t	m

Glacier Surface

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click

Browse to the model’s Application Libraries folder and double-click the file glacier_flow_2d_arolla02.txt.

Click

In the text field, type Glacier Surface.

Locate the section. In the table, enter the following settings:


Function	Unit
int2	m

In the table, enter the following settings:


Argument	Unit
t	m

Now define the flow rate constant A0 and the activation energy Q to calculate the flow rate factor (Equation 5).

In the toolbar, click

and choose .

In the window for , type A0 in the text field.

In the text field, type A0.

Locate the section. In the text field, type 263.15.

In the text field, type 3.985e-13.

In the text field, type 1.916e3.

Click to expand the section. Clear the check box.

Right-click and choose .

In the window for , type Q in the text field.

In the text field, type Q.

Locate the section. In the text field, type 60e3.

In the text field, type 139e3.

Definitions

Define the thermal conductivity, the heat capacity, and the viscosity of ice as variables depending on the temperature T. As T is calculated in , the variables have to be defined locally, that is, in , too.

Variables 1

In the window, expand the node.

Right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
k_ice	9.828[W/(mK)]exp(-0.0057[1/K]*T)	W/(m·K)	Conductivity of ice
cp_ice	146.3[J/(kgK)]+7.253[J/(kgK^2)]*T	J/(kg·K)	Heat capacity of ice

m_ice

In the toolbar, click

and choose .

In the window for , type m_ice in the text field.

In the text field, type m_ice.

Locate the section. In the text field, type abs(A0(T)*exp(-Q(T)/(R_const*T)))^(-1/3)*0.5 according to Equation 5.

In the text field, type T.

Locate the section. In the text field, type Pa*s.

In the table, enter the following settings:


Argument	Unit
T	K

Locate the section. In the table, enter the following settings:


Argument	Lower limit	Upper limit	Unit
T	260	280	K

Click

to plot the function and compare with the figure below.

Specify the pressure melting point temperature Tm according to Equation 11.

T_m

In the toolbar, click

and choose .

In the window for , type T_m in the text field.

In the text field, type T_m.

Locate the section. In the text field, type T_tp-betaCC*(p-p_tp).

In the text field, type p.

Locate the section. In the text field, type K.

In the table, enter the following settings:


Argument	Unit
p	Pa

Locate the section. In the table, enter the following settings:


Argument	Lower limit	Upper limit	Unit
p	0	120000	Pa

Click

Ambient Properties 1 (ampr1)

Now define the ambient properties. As there is no data for the Arolla glacier directly, select a station where similar weather conditions can be expected.

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Click .

In the dialog box, select in the tree.

Click .

In the window for , locate the section.

Find the subsection. In the table, enter the following settings:


Day	Month
01	01

Find the subsection. In the table, enter the following settings:


Hour	Minute	Second
00	00	00

Locate the section. From the list, choose .

Use automatic scaling of the axes to get a better view of the vertical changes within the glacier.

Axis

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Click

Geometry 1

Build the geometry using the interpolation functions defined above.

Parametric Curve 1 (pc1)

In the window, expand the node.

Right-click and choose .

In the window for , locate the section.

In the text field, type 5000.

Locate the section. In the text field, type s.

In the text field, type int1(s).

Click

Parametric Curve 2 (pc2)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 5000.

Locate the section. In the text field, type s.

In the text field, type int2(s).

Click

Line Segment 1 (ls1)

In the toolbar, click

and choose .

On the object , select Point 1 only.

In the window for , locate the section.

Find the subsection. Click to select the

toggle button.

On the object , select Point 1 only.

Click

Line Segment 2 (ls2)

In the toolbar, click

and choose .

On the object , select Point 2 only.

In the window for , locate the section.

Find the subsection. Click to select the

toggle button.

On the object , select Point 2 only.

Click

Convert to Solid 1 (csol1)

In the toolbar, click

and choose .

Select the object only.

Click the

button in the toolbar.

In the window for , click

Form Union (fin)

In the toolbar, click

Materials

Define the material ice in the next step. Introduce it as empty material node, first. As soon as the physics has been defined, the material node menu will show you which properties are needed for the simulation.

Ice

In the window, under right-click and choose .

In the window for , type Ice in the text field.

Laminar Flow (spf)

Now define the flow properties.

In the window, under click .

In the window for , locate the section.

Select the check box.

In the pref text field, type ampr1.p_amb.

Specify the rref vector as.
This guarantees that the pressure at the ice surface is equal to zero or the ambient pressure, respectively.


x	x
int2(x)	y

Fluid Properties 1

In the window, under click .

In the window for , locate the section.

Find the subsection. From the list, choose . For the enter 1e-15 [1/s].

Wall 2

In the toolbar, click

and choose .

Select Boundary 3 only.

In the window for , locate the section.

From the list, choose .

Select the check box.

In the Ls text field, type LSlip. This boundary condition is needed for the temperate glacier model. For the cold glacier flow it has to be disabled.

Open Boundary 1

In the toolbar, click

and choose .

Select Boundary 4 only.

In the window for , locate the section.

Clear the check box.

Open Boundary 2

In the toolbar, click

and choose .

Select Boundaries 1 and 2 only.

Pressure Point Constraint 1

In the toolbar, click

and choose .

Select Point 2 only.

In the window for , locate the section.

Clear the check box.

Heat Transfer in Fluids (ht)

In the window, under click .

In the window for , locate the section.

In the dz text field, type 500[m].

Initial Values 1

In the window, under click .

In the window for , locate the section.

In the T text field, type T_init.

Heat Flux 1

In the toolbar, click

and choose .

Select Boundary 3 only.

In the window for , locate the section.

In the q0 text field, type q_geo.

Open Boundary 1

In the toolbar, click

and choose .

Select Boundaries 1, 2, and 4 only.

In the window for , locate the section.

From the Tustr list, choose .

Heat Flux 2

In the toolbar, click

and choose .

Select Boundary 4 only.

In the window for , locate the section.

From the list, choose .

In the L text field, type 5000.

From the U list, choose .

From the pA list, choose .

From the Text list, choose .

Surface-to-Ambient Radiation 1

In the toolbar, click

and choose .

Select Boundary 4 only.

In the window for , locate the section.

From the ε list, choose . In the associated text field, type 0.97.

From the Tamb list, choose .

Temperature 1

In the toolbar, click

and choose .

Select Boundary 3 only.

In the window for , locate the section.

In the T0 text field, type T_m(p). This boundary condition is needed for the temperate glacier model. For the cold glacier flow it has to be disabled.

Materials

Having defined the physics, now fill in the empty expressions in the node.

Ice (mat1)

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Fluid consistency coefficient	m_pow	m_ice(T)	Pa·s	Power law
Flow behavior index	n_pow	1/n_ice	1	Power law
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	k_ice	W/(m·K)	Basic
Density	rho	rho_ice	kg/m³	Basic
Heat capacity at constant pressure	Cp	cp_ice	J/(kg·K)	Basic