Modelling subglacial water using Navier-Stokes ﬂow

(1)

Modelling subglacial water using Navier-Stokes flow

Masterarbeit im Studiengang Marine Geosciences Fachbereich Geowissenschaften

Universität Bremen von

Hans Leonard Rossmann March 14, 2015

First supervisor: Prof. Dr. Angelika Humbert

Second supervisor: Dr. Torsten Bickert

(2)

(3)

__________________________

Name

Declaration acc. to § 10 Paragraph 11 Common Part of the Master Examination Regulations

I hereby declare that I wrote my Master Thesis independently and that I did not use other sources and auxiliary means than the ones indicated.

This Master Thesis is not submitted in another examining procedure.

Place/Date: _____________________________________

Signature: _____________________________________

(4)

ii

(5)

Abstract

This Master thesis investigates the formation and evolution of subglacial drainage systems underneath ice sheets. The presence of water underneath ice masses has enormous effects on ice dynamics as it acts as a lubricant. Ice sheets play a big role in ocean circulation and climate dynamics. So far little is known about the formation and evolution of subglacial water. However, various elements of hydrological systems are suggested to exist ranging from thin water sheets to channels and lakes. This Master thesis attempts to fill this gap and creates a conceptual thermo-mechanical model of ice and water dynamics, which is solved using the finite element method.

While the Stokes-equation is used for ice, the Navier-Stokes equation is introduced for the simulation of the water flow. Thus the acceleration terms are included. A coupling of water and ice flow has to overcome the problem that water flows on short time scales while ice on long once. This problem is solved with an adaptive timestepping mode. In order to evolve a water layer the melt rates are calculated by a state of the art enthalpy formulation, using the advantage of solving for temperature and water content at the once.

The results demonstrate that it is possible to model both ice and water flow within one approach using the Navier-Stokes flow and the enthalpy formulation. The experiments investigate the effects of geothermal heat flux and inflow velocities on the melt rate behaviour. The results show that thin water layers develop into a channelized flow form when a certain threshold in the water velocity is overcome. Moreover, cavities form in regions of increased geothermal heat flux. With increasing inflow velocity the cavity propagates further in the direction of flow. The results demonstrate a strong relationship between water flow velocities and changes in melt dynamics. Further, it illustrates that the subglacial water system is under constant change. Hence, this study gives an important insight into melt dynamics under ice masses and their evolution.

iii

(6)

iv

(7)

1. Introduction

Glacier ice covers nearly 10% of the today’s Earth surface (Cuffey and Paterson, 2010).

These glaciated areas play a major role in the Earth system. Changes in the cryosphere have effects the Earth’s energy budget via changes in albedo or changes in sea level due to accumulation or release of fresh water. The effect of fresh water into the ocean has two major effects. First, when fresh water enters into the regions of deep water formations (e.g. the North Atlantic) it changes the whole Meridional Overturning Circulations (MOC) due to density decreases of the water masses (Rahmstorf, 2006; Delworth, 2008). This triggers changes in climate conditions and heat transport via ocean circulations (Rahmstorf, 2006). Second the change in sea level due to mass loss or gain from the land ice sheets.

The Greenland and the Antarctic ice sheets sum up to a potential seal level rise of 65.66 m (Solomon et al., 2007), when completely molten. The ice sheets will not melt completely in the near future but an increased volume change has been recorded in the last years. The volume loss of Greenland and Antarctica combined is -503 ± 107 km³yr⁻¹ between 2003 and 2009 (Helm et al., 2014). In Greenland the mass loss divides into ∼50% by melting and ∼50% by calving. Whereas, in Antarctica the major mass loss is due to calving (ca.

98%) and only 2% due to melting (Shepherd et al., 2012). Therefore, it is on of interest to analyse the contribution of these mass losses.

On an ice sheet mass is accumulated via precipitation in form of snow. It can lose mass due to sublimation, melting including run off and calving and transport over the grounding line into an ice shelf. Since the last two are the dominant process in Antarctica, it is an aim of this study to investigate effects, which contribute to the mass losses. In an ice sheet mass is transported via outlet glaciers or ice streams form the interior to the edges.

Here, glaciers either terminated on land, on sea or in ice shelves. For the later one, the transport of ice over the grounding line is important. The grounding line is the region where ice stream or glaciers lose contact to the bedrock and the ice gets afloat. At that point the ice already contributes to sea level change due to the hydrostatic equilibration.

The transport of ice in these outlet glaciers and ice streams is a result of their ice dynamics. Ice can flow via creep flow, a deformation on a grain scale basis (Cuffey and Paterson, 2010), which has flow rates ranging form1m y r⁻¹to several hundred meters per year(Rignot et al., 2011) depending on the applied stress. Whereas, the high end velocities in ice streams, which can reach up to kilometres per year (see Fig. 1.2), are not only due to creep deformation of the ice itself but due to a combination of creep flow and sliding over the bedrock. The sliding occurs if the ices masses are not frozen to the ground (so called temperate base) and liquid water (subglacial water) is present at the base. Liquid water

(10)

2 Chapter 1. Introduction

Grounding line

Ocean Ice divide

Geothermal heat source

Figure 1.1.: Schematic drawing of mass transported in ice sheets and velocity fields.

has an enormous lubrication effect on the ice. Ice mass dynamics are strongly influenced by the conditions at the bedrock. The subglacial water has different forms of appearance underneath an ice mass, which is described in the next section. As the ice dynamics is influenced by the presence of water, so is the heat transfer and temperature regime of an ice mass. If subglacial water flows between the bedrock and the ice mass, it can transport heat faster than ice. This has effects on the melt dynamics at the base and therefore on the evolution of the subglacial drainage systems. Hence, it is benefit to study the effects of melt dynamics linked to water flow to give assumption about the formation and effect of subglacial water systems.

Up to now only few field surveys have been done, as it is difficult to localize the water bodies underneath ice masses. There have been dye trace measurements, analysis of glacier relict areas, bore hole measurements and camera surveys (Benn and Evans, 2010). These studies suffer a long-term evolution record of these systems, which is hardly achievable.

Nevertheless, since the last decade the existence of subglacial water under ice sheet has been proven (Kleiner and Humbert, 2014; Benn and Evans, 2010; Cuffey and Paterson, 2010; Siegert, 2000) and its important effects on ice dynamics and melt dynamics is not negligible (Fricker et al., 2007). In order to examine the changes in subglacial hydrological systems over a longer period numerical models are necessary and can achieve this goal. They can predict changes in these subglacial systems, according to the implemented equations.

The aim of this study is now to create a numerical model, which links ice flow, water flow and heat transfer in one. This enables to study the evolution of subglacial systems. The

(11)

1.1. Subglacial water 3

Figure 1.2.: Surface velocities of Antarctica with clearly visible high velocities in outlet glaciers and ice streams (Rignot et al., 2011)

Navier-Stokes equation is common in Computational Fluid Dynamics (CFD). It can predict the flow in laminar and turbulent regimes. Further, it will enable to model the high viscous flow of ice and the low viscous flow of water. Using the same equation to model different flow brings the advantage that these flow can be modelled at the same time. Up to now in numerical models only parameterizations of the subglacial water flow are implemented like in Kleiner and Humbert (2014). This master thesis therefore fills the gap of modelling the evolution of subglacial systems with the help of a numerical model. The following quote by Benn and Evans underlines the need of further investigation in glaciers hydrology: "In the last decade or so, considerable advances have been made in our understanding of glacial hydrological systems, and our ability to represent them in glacier models has improved immensely. However, the statement made by Arnold et al. (1998), that much work remains to be done in terms of field observations, modelling and development of theory of glacier hydrological systems, is no less true today than when it was written."(Benn and Evans, 2010).

1.1. Subglacial water

Subglacial water describes all possible water flows beneath ice masses or glaciers. The water origins either from above the glacier, created by melting processes or from geothermal heating or frictional heating by the ice movement itself or strain heating. This raises the temperature to the pressure melting point. All surplus energy is than used for melting of

(12)

4 Chapter 1. Introduction

the ice. Once the water is present it can either gather in hydro potential sinks and form lakes underneath the ice sheet or it can be transported. There are five basic drainage systems, which can transport water. Bulk movement of water with deforming tilt, Darcian pore water flow, Dendritic channel network, Braided canal network, Linked cavity and Thin water film (see Fig.1.3(b) to 1.3(d)) (Benn and Evans, 2010; Glasser, 2013; Fountain and Walder, 1998).

(a) Braided canal network (b) Dendritic channel network

(c) Linked cavity (d) Thin water film

Figure 1.3.: Schematic drawings of subglacial drainage systems (modified after Benn and Evans (2010))

From Alpine glacier it is know that changes in these drainages systems control the ice dynamics. This linkage is well studied for surging glaciers. This type of glacier has a periodical fast and slow flow phase. One of the intensively studied glacier is the Variegated Glacier, in southern Alaska, (Murray, 2003). The surge behaviour of this glacier is linked to the change form a linked cavity system during the fast flow to a channelised drainage system which terminates the fast flow period (Murray, 2003; Kamb, 1987). These drainage systems are present underneath ice sheets as well. So the following thought experiments are potential effects of the named drainage systems. The effect of a thin film of water could lead to an increased ice velocities, since the enhance lubrication effect Fountain and

(13)

1.2. Aims and Objectives 5

Walder (1998). Linked cavities could have the same effect as the thin film as they appear over a broad region. They are supposed to be more stable than the thin film. In these both systems the water flow is slow. They are regarded as inefficient drainage systems.

Therefore, the heat in these both systems is mainly transported via a diffusion process.

In contrast to them the channelized systems are fast flowing systems, where water is drained rapidly. Therefore, heat is transported mainly via advection. Further, turbulent flow regimes are able to establish, which influences the heat transport, too. Their localized appearance should have a lower effect than the linked cavities or the thin film has on the ice dynamics. Overall the subglacial drainage systems can shape the thermal and the dynamic regime of glaciers and ice masses.

1.2. Aims and Objectives

The aim of this master thesis is to generate a numerical model, which describes the ice and water flow at the same time. Further a heat transfer module shall be implemented in order to model melt and refreezing process in subglacial systems.

The coupling of heat transfer and flow modelling with the Navier-Stokes equation has not jet been implemented in ice models. There is the difficulty of coupling the fast flow of water with the slow of ice. Therefore, one part of the thesis investigates the possible linkage between these two flow types. Further the state of the art Enthalpy formulation from (Aschwanden et al., 2012) is used to model the thermodynamics. The aim of this thesis is to create a conceptual two-dimensional numerical model, which reflects the conditions of an ice sheet and incorporates the named physical principles. Additionally different experiments shall give an insight into melt dynamics underneath ice masses.

The key scientific research questions are here: How can subglacial channels form underneath ice masses? What are the effects of flowing water on the heat transfer in subglacial systems?

How do cavities form and evolve underneath the ice. Since these in numerical model and field observations these question have not be answered, this Master thesis attempts to show that the flow velocities in subglacial drainage systems have enormous impact on their evolution and creation. The presents of water at the base and its velocities are in constantly changing system.

(14)

(15)

2. Theory

In this Chapter the theoretically foundation is given for the numerical model. On long time scales glacier ice behaves like a fluid. Therefore, it can be described using continuum mechanics. With the help of fluid dynamics, which is a branch of continuum mechanics, the governing equations describe the flow and thermal state of ice and water. The governing equations are the mass balance equation (section 2.1.1), momentum balance equation (section 2.1.2) and the energy balance equation (section 2.1.4). Constitutive equations (Section 2.2) are needed to complete the system. Here only the most important equation are presented as a full derivation of the equation would exceed the scope of this thesis.

Furthermore, they are described in detail by (eg. Greve and Blatter (2009)). The theory of the enthalpy formulations (heat transfer) was first described by Aschwanden et al. (2012) and only the most important equations are displayed here. The numerical method of finite elements is applied to solve the partial differential equations.

2.1. Balance Equations

2.1.1. Mass balance

Pure glacier ice and water are incompressible fluids. The mass balance of an incompressible continuum can be written as:

d i v(v) = 0, (2.1)

which states that the velocity field v is source and sink free.

2.1.2. Momentum balance

The momentum balance engages from Newton’s second law. The change of momentum over time has to be in equilibrium with all forces applied on the body. These forces can either be external volume forces f acting on each element of the body such gravity or Coriolis force, or internal stresses σ which act on the body’s boundary surface. The momentum balance for momentum density (ρv) reads:

∂ρv

∂t +d i v(ρv v) = d i v(σ) +f. (2.2)

(16)

8 Chapter 2. Theory

The equation 2.2 is as well called the Navier-Stokes equation. It can be rewritten as followed in its incompressible form:

ρ·(∂_tv+v∇v) =−∇p+ρg+η∇²v. (2.3) The Navier-Stokes equation links the velocity term (∂_tv) and the turbulent term (v∇v) of a fluid parcel on the left hand side to the pressure derivatives (−∇p), the gravity (ρg) acting on the parcel and the fluids viscosity term (η∇²v) on the right hand side. The ratio between inertial forces (−∇p +ρg) to viscous forces (η∇²v) is expressed by the dimensionless Reynold’s number (Re) for a given flow condition. When the Re number exceeds a critical value the flow type changes from laminar flow to turbulent flow.

Water: The Re number of water is high since the viscosity is low (η_{w ater} ≈ 1.78 · 10⁻³P a s). Water can flow in a laminar or turbulent regime depending on the flow velocity.

In order to calculate the flow of water the whole Navier-Stokes equation has to be solved, thus

Ice: On the other side, ice has a high viscosity, 17 magnitudes larger than water, and therefore a low Reynold’s number (Re ≈10⁻¹⁰) (Lliboutry, 1987). Thus the inertial forces, on the left hand side of the equ. 2.2, can be neglected. Equation 2.2 than simplify to:

d i v(σ) + f = 0, (2.4)

with the Cauchy stress tensor σ, which defines the state of the stress at a point inside a deformable material. In fluid dynamics it is common to split the tensor into a deviatoric partσ^D and a pressure part pI:

σ = σ^D−pI, (2.5)

where I is the identity tensor and p = −¹₃tr (σ) denotes the pressure. The volume force f contains gravitation force, from the rotating Earth centrifugal force and Coriolis force.

The Coriolis force can be neglected, since the glacier ice is flowing with slow velocities.

The gravitational and centrifugal force are combined into the effective force of gravity ρg.

The vector of the gravitational acceleration has the form ofg = (0,0,−g), withg = 9.81 m s⁻². This leads to the Stokes equation:

d i v(σ−pI) = −ρg (2.6)

2.1.3. Enthalpy formulation

For the description of the energy balance the approach by Aschwanden et al. (2012) is followed in order to determine cold, temperate ice and water with the same equation. As

(17)

2.1. Balance Equations 9

well as above the most important equation are shown from the papers by Aschwanden et al. (2012) and Kleiner et al. (2015). In thermodynamics literature (Moran and Shapiro (2006)) the specific enthalpy E is defined as E =U+p/ρ, where U is the specific internal energy and p the pressure. Since the work associated with changing the volume of the material is not included, it can be set E = U with the SI unit J k g⁻¹. In the following the enthalpy is described as functions of temperature T depending on the phase of the material. For cold ice, below the pressure melting point (T_pmp), the specific enthalpy of ice E_i is defined as:

E_i = Z T

T₀

C_i(T)d T , (2.7)

where C_i(T) is the heat capacity of ice and T₀ is the reference temperature. If this temperature is lower than all modelled temperatures the enthalpy will be always positive.

The specific enthalpy of liquid water E_w reads as:

Ew =

Z Tpmp(p)

T₀

Ci(T)d T +L+ Z T

Tpmp(p)

Cw(T)d T , (2.8) where C_w(T) is the heat capacity of water and L the latent heat of fusion. A further advantage of the enthalpy formulation is that temperate ice (T = Tpmp) can be described with a resulting water content (ω). This leads in case of temperature below the pressure melting point (T < T_pmp) to a water content of zero (ω = 0), or for temperature at the pressure melting point (T = T_pmp) to water content between zero and 100%

(0 ≤ ω ≤ 1).The specific enthalpy of mixture, including cold ice, temperate ice and liquid water reads as:

E = E(T , ω, p) = (1−ω)·E_i(T) +ωE_w(T , p). (2.9) These leads to two cases:

E =

( E_i(T), T < T_pmp(p)

E_s(p) +ωL, T = T_pmp(p)and0≤ω≤1 . (2.10) There are the following transfer rules for temperature and water content to enthalpy:

E(T , ω, p) =

( ci(T −Tr ef), i f E < Epmp

E_pmp+ωL, i f E ≥ E_pmp (2.11)

2.1.4. Enthalpy balance

The enthalpy balance for ice reads (Aschwanden et al., 2012):

∂(ρ_iE_i)

∂t = −∇ ·(ρ_iE_iv +q_i) +Q_i−X

w

, (2.12)

(18)

10 Chapter 2. Theory

and similar for liquid water:

∂(ρ_wE_w)

∂t = −∇ ·(ρ_wE_wv +q_w) +Q_w +X

w

, (2.13)

where P

w is the enthalpy exchange rate between components,Q_i and Q_w are the dissipation heating rates. The advective and non-advective enthalpy fluxes of the ice component are ρ_iE_iv and q_i, respectively and same for the liquid component. From equ. 2.12 and 2.13 leads to the balance for the total enthalpy flux:

ρd(E)

d t = −∇ ·q+Q, (2.14)

where q = q_i +q_w. The rate, Q , at which dissipation of strain releases heat into the mixture is:

Q = tr(σ·σ^D) (2.15)

2.2. Constitutive equations and rheology of ice

The balance equations, derived above, are principles. They need determining equations in order to represent specific materials. These constitutive equations close the above balance equations since there are 3 sets of equations from mass balance and momentum balance but 6 unknowns (Pressure, velocity and stress).

2.2.1. Incompressible Newtonian flow

Any isotropic fluid can be described by the relation between the deviatoric viscous stresses (σ^D) and the strain rates tensor ().˙

σ^D = 2η,˙ (2.16)

where η is the viscosity and the strain rate tensor , written in component form:˙

˙ _{i j} = 1

2 ∂u_i

x_j + ∂u_j x_i

.

If the viscosity is constant, which is the case for liquid water at given temperature, than the stress is linear proportional to the strain rate. Fluids, which have a linear relationship between strain rate and stress, are called Newtonian fluids.

2.2.2. Incompressible Non-Newtonian fluids and the rheology of ice

In glaciers ice does not occur as on single crystal. It is rather build up a great assemblage of crystals (so called grains or crystallites). The sizes are on a millimetre to centimetre scale. On an average the orientation of the single crystals are randomly distributed, so

(19)

2.2. Constitutive equations and rheology of ice 11

that an isotropic behaviour is assumed for glacier ice. If a block of ice is sheared with a constant stress (τ) a shear angle (γ) will develop as show in Fig. 2.1. The angle measured over time gives a creep curve (shown in Fig. 2.1). First the ice responds elastic on applied force followed by a phase called primary creep, where the shear rate (γ) decrease with˙ time. At one point the minimum shear rate is reached and stays constant (the secondary creep phase). Before the last phase the tertiary creep sets in, which is the case at high homologous temperatures and long time scales, an acceleration phase can be recognized (Greve and Blatter, 2009). These phases can be attributed to dynamic recrystallization of the crystallites.

τ

γ

Time, t

Shear angle , γ

t₀

elastic primary

secondary accelerating tertiary

t₁

a) b)

Figure 2.1.: a) a block of ice is sheared (in a simple shear form) showing the resulting shear angle γ,

b) Diagram of the evolution of the shear angle over time, modified after Greve and Blatter (2009).

2.2.3. Glen’s flow law

For glacier ice laboratory and field data have shown that the connection between strain rate and stress depends on a power law (Glen, 1955). The generalized form by Nye (1957) reads as:

˙

_{i j} = E A τ_e^Dⁿ⁻¹τ_{i j}^D, (2.17) where (τ_e^D) is the effective deviatoric stress, the second invariant of the deviatoric stress tensor. The factor A is called the flow rate factor and is dependent on temperature and water content (A = A(T , ω)). E is called the enhancement factor and parameterizes other physical contributions, like fractures or impurities. The power law exponent n is set to 3 by convention (Cuffey and Paterson, 2010). In the model A and E are kept constant for simplicity reasons. The inverse form of Glen’s law reads as:

τ^D = 2η˙ with η( ˙) = 1

2(E A)⁻¹ⁿ˙

1−n

en , (2.18)

(20)

with the effective strain rate ˙_e = q1

2tr˙². Ice behaves like a Non-Newtonian fluid as the relation between stress and strain is non linear and the viscosity is dependent on the effective strain rateη =η(T ,˙e). The viscosity of ice ranges around10¹³P a s(atT = 0^◦C andτ_e = 100k P a) to 10¹⁷P a s (at T = −20^◦C andτ_e ≈10k P a). In comparison motor oil has a viscosity of 0.1P a s and the Earth’s mantel 10²¹P a s (Greve and Blatter, 2009).

2.3. Boundary conditions

The previously derived balance equations are only applicable if the thermodynamic field is sufficiently smooth and thus continuously differentiable. This assumption is not valid on the surfaces of the model and the surroundings. These interfaces, which experience a discontinuity in a physical quantity, are called singular surfaces. Jump conditions or boundary conditions have to be formulated in order to close the equations. There are three types of boundary conditions. The Neumann boundary condition prescribes the derivative of a quantity at the interface. If the value of the field quantity is set to a certain value, it is called a Dirichlet boundary condition. There is also a third type of condition, which is a mixture of the both and is called Robin condition. Four boundary conditions are need for the description of the model. They are the boundary condition for the boundaries between ice/atmosphere, ice/water, water/lithosphere and ice/ice. They are described in the following sections and are schematically illustrated in Fig. 2.2.

Ice domain

Water domain

Boundary Ice water Top / Atmosphere

Bottom / Lithosphere

Side A / Ice Side B / Ice

α

a _b

h_b

H_w

X Y

L H

Input flow Output flow

Periodic flow condition Periodic flow condition

q

⁺

= q

⁻

E = Epmp

Free slip

No slip Free slip

Free slip

Figure 2.2.: Schematic drawing of the geometry setup with implemented boundary condition switches. Please not that the x-axis is pointing horizontally and y-axis vertically. This is against convention but results from Comsol interior labling.

(21)

2.3. Boundary conditions 13

2.3.1. Ice surface

The boundary between ice and atmosphere is regarded to be a singular surface and its implicit form can be written as followed:

F_s(x, t) =z −z_s = z −h(x , y , t) = 0, (2.19) where z_s = h(x , y , t) is the position of the ice surface. The kinematic boundary condition reads as:

∂ z_s

∂t +vx

∂z_s

∂x +vy

∂z_s

∂y +vz = Nsa_s^⊥, (2.20) where thea_s^⊥ is the accumulation perpendicular to the surface andN_s is the gradient norm:

N_s = |∇F_s|= 1 + ∂z_s

∂x 2!¹₂

. (2.21)

Since the stress applied by the atmosphere and wind on the surface are very small in comparison to normal stresses ice and thus negligible. The ice/atmosphere boundary can be seen as a traction free surface:

t·n = 0. (2.22)

This is called the dynamic boundary condition. The thermodynamic boundary condition is described by a Dirichlet condition. The enthalpy is prescribed by a mean annual value E_s(x , y , t):

E = E_s(x , y , t) (2.23)

2.3.2. Ice base

Analogous to the singular surface at the ice/atmosphere boundary the kinematic boundary condition for the ice base (F_s(x, t) = z −z_b = z −h(x , y , t) = 0) can be derived and reads:

∂ zb

∂t +v_x ∂zb

∂x +v_y ∂zb

∂y +v_z = N_ba^⊥_b, (2.24) where z_b is the height of the basal interface andN_b is the norm of the normal vector.

N_s = |gr ad Fb|= 1 + ∂z_b

∂x ²!¹₂

(2.25) In the case that water is present at the base, the geothermal heat flux (q_geo) enters into the water domain. In the case of ice at the base, the thermal conductivity parameter (Kj) changes to the one of ice Thermodynamic boundary at the boundary between water and lithosphere is therefore:

−K_,· ∇E·n_b = q_geo with j = [w , i]. (2.26)

(22)

2.3.3. Ice water interface

Since this thesis deals with subglacial drainage systems and in order to evolve them melting and refreezing plays an enormous role. The basal melt rate (a_b) is calculated by the difference of heat fluxes at the interface (see equ. 2.29).

In the following it has be distinguished between the different conditions, which can occur underneath an ice sheet. The boundary conditions are linked to the presence of water and if the ice is temperate or not.

For the thermodynamic boundary condition there is a decision chart give by Aschwanden et al. (2012), which is adapted for this model. These decisions read as followed:

Cold base (dry): The glacier is cold at the base and or the water layer is smaller or equal the minimum water thickness (plus a small value for numerical stability) ((H_w ≤hb+ ∧ E < E_pmp)∨(H_w ≤hb+ ∧ E ≥E_pmp)) , then a Neumann boundary condition holds and sets fluxes on both sides equal:

−n⁺·(k∇E)⁺ = −n⁻·(k∇E)⁻ (2.27)

E⁺ = E⁻, (2.28)

where the superscript ⁺ and⁻ define respectively one side of a boundary. Here it can be referred to ice side ⁺ and water side ⁻.

Cold or temperate base with water: If the water layer is larger than the original minimum water thickness(H_w > hb+ ∧ E < E_pmp)∨(H_w > hb+ ∧ E ≥E_pmp), the conditions at the interface changes to a Dirichlet condition T = T_pmp. Therefore the interface condition represents a third type condition, which is updated at every timestep. Since the temperature is kept constant at the boundary, all exceeding energy is converted into melting. The melt rate is calculated by difference of heat fluxes at the boundary between ice and water ∂Ω_i−w and reads as followed (the unit is mm a⁻¹ in ice equivalen ):

a_b = −(q_i −q_w)n_b

L·ρ_i , (2.29)

where L is the latent heat of fusion ρ_i the density of ice. The terms q_i and q_w are the conductive heat fluxes of ice and water, respectively. They result from the enthalpy version of Fourier’s law and are calculated by:

q_j =−Kj∇E, (2.30)

where j =i , w indexing ice or water and Kj =kj/cj is the thermal conductivity in enthalpy form.

(23)

2.3. Boundary conditions 15

Further the flow conditions have to change for the ice. The boundary condition switches form a non-slip condition to a slip condition. It is achieved by manipulating the friction coefficient at the boundary. As soon as the water thickness rises above the original limits the friction coefficient (β²) decreases significantly. In case of a frozen bed the friction coefficient is infinitely large and hence no sliding occurs. The friction parameter (β²[P a a m⁻¹]) is related to the basal drag τb (sum of all basal resistance at the base) and the basal velocity (v_b) (Pattyn, 2010) by

τb =β²·vb (2.31)

2.3.4. Lateral boundaries

Since the model is representing the interior of an ice sheet, the lateral boundaries have to represent these conditions. Therefore the boundary conditions are chosen to simulate the as if multiple copies of itself would surround the domain. Since the domain lateral boundaries are fixed in space, a kinematic boundary condition can be omitted. The dynamic boundary conditions are chosen to represent a periodic wall and read as:

v⁺ = v⁻ (2.32)

p⁺ = p⁻, (2.33)

where v⁺ and v⁻ are the velocities of the destination and source boundary, respectively and p⁺andp⁻ the pressure of the destination and source boundary, respectively.

For the thermodynamic boundary condition at the lateral boundaries an insulating condition is chosen, so that the there is no enthalpy gradient over the interface. That is in contradiction to the cyclic environment proposed above. Nevertheless it brings a better numerical stability into the model. Thus the thermal insulation reads as:

−n·(−k∇E) = 0 (2.34)

(24)

(25)

3. The numerical model

3.1. Implementation

The implementation of the model into the commercial software COMSOL Multiphysics^c is shown in this Chapter. This program offers the ability to solve Partial Differential Equations (PDE), like the heat transfer equation or the Navier-Stokes equation (see Chapter 2), on the numerical approach of the Finite Element Method. COMSOL Multiphysics^c can be used via a Graphical User Interface (GUI), a live-link via Matlab^c or via JAVA API. In the following the use of the GUI will be described in detail. The GUI enables a quick and easy handling of models and their implementation. Since, COMSOL Multiphysics^c is a commercial software is not possible to see into the source code and this is a disadvantage in some cases for detailed manipulation in the model.

The GUI (see fig 3.1) is divided into three main parts theModel Builder, theNode properties and Graphics. The model is created by successively adding nodes (branches) in the model builder tree in the Model Builder section. By default the branches Global Definitions, Model, Study and Results are pre-set. The branch Global Definitions is used to define global valid parameters, shown in Table 3.1. If not mentioned over wise, these parameters and variables are used as default values for the model setups.

The Model builder tree is subdivided into several subbranches/subnodes. Here the nodes Definitions, Geometry and Mesh are set as default. It is a strength of Comsol to add pre-described PDE (the so called Physics) to the model at this point. The Physics nodes Heat transfer in Fluids, Laminar Flow 1 and 2 andMoving Mesh were added to the model tree. In the node Definitions variables can be defined and further selections of certain areas of the model geometry can be set, which improves the later handing.

3.1.1. Geometry

The geometry of the model defines the domain, where the PDEs are solved. The geometry node offers the ability to design all possible shape. This can be achieved by assembling different geometry object like circles, rectangle, lines or spline curves. These objects can be transformed in a more complex geometry by operators like union, split, difference and intersect. It also possible to load externally built geometries by a Computer Aided Design Program into Comsol (only possible with extended licence).

For subglacial water models a rather simple geometry is chosen, in order to represent a segment of an ice sheet. The shape is approximated by a two dimensional square (’ice

(26)

18 Chapter 3. The numerical model

Figure 3.1.: Screenshot of Comsol Multiphysics^c GUI

Table 3.1.: Default physical global parameters name expression unit description

rho_ice 910 kg/m³ Density of ice rho_water 1000 kg/m³ Density of water

g 9.81 m/s² gravity acceleration

spy 31556926 seconds per year

L_f 339000 J/kg Latent heat of fusion q_geo 42 mW/m² geothermal heat flux

n 3 Exponent in Glen’s Flow law

k_ice 2.1 W/m/K thermal conductivity of ice cp_ice 2009 J/kg/K specific heat capacity of ice cp_water 4216.278 J/kg/K constant specific heat capacity

of water at T_0

k_water 0.556 W/m/K constant thermal conductivity of water at T_0

cube’) with the height H and length L (all geometry linked parameters and variables are found in Table 3.2). This block is divided into two subdomains like shown in schematic fig 2.2. The creation of two domains is necessary due to topological reasons. Since the water is non existing at the start of the simulation it still needs an object to evolve itself.

Otherwise the water domain would be generated out of nothing, which is not possible.

The two sub domains are built as followed. The ice domain is overlaying the water domain and has an original thickness of h_t. Whereas, the water domain has a initial thickness of hb. These values were chosen in order to have the lowest possible thickness for the water domain, which is still in the building tolerances. Moreover, it fulfils the topology

(27)

3.1. Implementation 19

requirements and has the least effect on the starting conditions of the numerical model.

During the simulation run the water domain may evolved and change geometry according to the melt rate (a_b) at the boundary between ice and water. Further the domains have a slope with the angle α. Creating the slope and the internal boundary between ice and water is achieved by inserting a line according the function z_top, z_{i nter}, z_bottom in order to represent the boundaries ice/atmosphere, internal ice/water and ice/lithosphere respectively. Since all parameters are globally defined, they can be changed according to the model setups. Further, please note that the axis labelling is against convention. The x-axis is pointing in horizontal and the y-axis in vertical direction.

Table 3.2.: Default geometry parameters and variables name expression unit description

L 1000 m Length of block

H 1000 m Height of block

α 0...0.1 ^◦ Angle of the slope slope tan(α) Gradient of the slope

h_b 3 mm Initial height of water domain h_t H−h_b m Initial height of ice domain z_top - slope·x +H m Function of top boundary z_{i nter} - slope·x +h_b m Function of internal boundary

zbottom - slope·x m Function of bottom boundary

3.1.2. Switch values

Figure 3.2.: Smoothed Heaviside function by the Com- sol build in function flc2hs of the liquid parameter smoothed over 10⁻⁵ m

As mentioned above the water domain is evolving in time. This leads for example to a change in the boundary condition at ∂Ω_i−w as shown in Chapter 2.3.3. Switching condition during the simulation requires an additional processing step during calcula- tions. This change is enabled by a parameter switch.

The so called liquid parameter is a logical parameter indicating if the water domain is larger than the initial height of h_b.

liquid=

( 1 , H_w > h_b

0 , H_w ≤h_b , (3.1) where H_w is the thickness of the water domain. This parameter enables to control changes of other physical parameters and initialisation of processes. It

(28)

triggers the switch of the Neumann to the Dirichlet condition at the internal boundary, the change from the physical properties of ice to water in the water domain and the start of the water flow. The explicit use is described in the sections below when theliquidswitch is of need. Since numerical models do not cope well with hash changes the liquid parameter has to be smoothed. This is achieved by applying a smoothed Heaviside function the liquid parameter. Comsol offers the so calledflc2hs function which is a built in smoothed Heaviside function with a continuous second derivative without overshoot (COMSOL Inc, 2015). In fig 3.2 the smoothing of liquid is shown with a smoothing over 10⁻⁵ m. The implementation is done via adding a variable branch under the Definition branch in to Comsol GUI.

3.1.3. Mesh

The mesh is the discretisation of the domain, which enables solving the PDEs at the nodes. Therefore the mesh has a huge importance for the model. It has influence on the stabilization, convergence and computational cost. The mesh is chosen to have a fine resolution in the water domain, since the gradients are larger here due to the fast flowing water.

The meshing created with the Comsol built-in meshing operator Mapped, which is for building rectangular mesh elements. In order to define the amount of elementsDistributions are defined on the domain borders. The water domain has 10 elements in y directions and 250 in x directions, creating 2500 elements in the water domain with a height of 0.3 mm and a width of 4 m. The ice domain has the same number of elements in x directions and 25 in y direction. This creates 6250 mesh elements. They are not equal sized like in the water domain but increase in size in y directions with a factor of 0.01, since smaller elements are needed at the boundary between ice and water. In the region near the boundary between ice and atmosphere no huge changes are expected and therefore not as many elements are needed.

3.1.4. Arbitrary Lagrangian Method (ALE), Moving Mesh

In order account for an evolving water column the internal boundary has to move according to the melt rate at the interface between water and ice. The melt ratea_b in ice equivalent is calculated according to equ. 2.29. The melt rate a_b is applied on the boundary and is constrained by the minimal thickness of the water layer (hb) and the direction of the melt rate. Therefore the moving mesh velocity v_a_b reads:

v_a_b =

( 0 , H_w < hb ∧ a_b≤0

a_b·n , H_w ≥h_b ∧ a_b >0 , (3.2)

(29)

If the∂Ω_i−w reaches the minimum thickness hb then ∂Ω_i−w is constrained by:

∂Ω_i−w(y) =

( ∂Ω_i−w(y) , Hw > hb ∨ ab>0

∂Ω_i−w(hb) , H_w ≤h_b ∨ a_b≤0 (3.3) This constraint avoids penetrations of the ice domain into water domain and the lithosphere.

Moreover, it triggers the liquid parameter, which controls several other switches as mentioned above.

The calculation of the new mesh coordinatesx , y is done with linear mesh elements. Since the input velocity (melt rate v_a_b) is not always smooth across the elements, the mesh is smoothed with a Laplace algorithm.

The deformation of the mesh results from the calculated melt rate at the boundary between ice and water (∂Ω_i−w). In Comsol the node Fixed MeshandPrescribed Mesh Displacement are default in theALE branch. In order to allow mesh movement theFree Deformationhas to be added to the branch. Since only the boundary ∂Ω_i−w shall move free, it is necessary to set the movement of the other walls to zero with aPredefined Mesh velocity node.

This includes the boundaries top, bottom and sides (see Fig. 2.2).

A second Predefined Mesh velocity node is added to the branch in order to apply the melt rate on the internal boundary ∂Ω_i−w according to equations 3.2. The no penetration condition (equ. 3.3), are implemented by a Pointwise constraint node. This prohibits the boundary ∂Ω_i_−w to lower itself beneath the minimum height of h_b.

3.1.5. Heat transfer in fluids

Figure 3.3.: Top branch of the Heat transfer node The implementation of the heat flux in the model is done via

the Physics nodeHeat transfer in fluids. Here the enthalpy E is solved in the domains. As described in the Chapter, 2.1.3 and 2.1.4, transformation rules are needed to convert the enthalpy into temperature, which is the more comprehensible physical quantity for glaciologist and marine geologists. These transformation rules are implemented under a variable node in theLocal Definitions and listed in Table 3.3. Furthermore, the parameters for the boundary conditions, the thermal conductivity of ice and water and the calculation of the temperature / enthalpy at the pressure melting point (T_pmp and E_pmp, respectively) are listed in the same table.

The general picture of the model regarding the heat transfer is that there is a heat source at the bottom and a constant temperature at the top. Furthermore, the two domains have a different parameters resulting from the different phases in

(30)

Table 3.3.: Heat transfer parameters and variables

name expression unit description

T_ref 223.15 [K] Reference Temperature

T_0 273.15 [K] Melting point of water

under normal pressure

T_top 268.15 [K] Temperature at the top

E_start cp_ice*(T_top-T_ref) [J/kg/K] starting Enthalpy

my_p rho_ice*(H-y)*g pressure

beta 7.9e-8 [K/Pa] Clausious-Clapeyron

constant

T_pmp T_0-beta*my_p Temperature at

pressure melting point E_pmp cp_ice*(T_pmp-T_ref) [J/kg/K] enthalpy at Pressure

melting point K_0 k_ice/cp_ice/10 [kg/J*K] moisture mass diffusivity, conductivity of temperate ice

K_c k_ice/cp_ice [kg/J*K] Conductivity of cold ice

K_water k_water/cp_water [kg/J*K] constant thermal conductivity of water T if(E<E_pmp,

E/cp_ice+T_ref,T_pmp)

[K] Temperature omega if(E>=E_pmp,

(E-E_pmp)/L_f,0)

water content

the model. In order to create such a model the following adjustments in theHeat transfer node have to be made:

• The topbranch Heat transfer in fluids is for general settings (compare Fig. 3.3 ).

The domains one and two (all domains) are the regions where the PDE of the heat transfer is solved. They are set in the window domain settings. Under equations the equation type is chosen to be study dependent. In our case this is a time dependent form. In contrast if the stationary is chosen, the time dependent part would be set to zero (c_pρ ∂E/∂t = 0). Furthermore, Streamline and Crosswind diffusion are enabled as stabilization techniques. They avoid numerical instability caused by temperature advection. The streamline stabilisations adds numerical diffusion in streamline direction (Codina, 1998), whereas thecrosswind diffusion acts in orthogonal direction to the streamline (Hauke and Hughes, 1998). Under the section Discretization linear elements are chosen for the FEM shape functions. The option Compute boundary fluxes and Apply smoothing of boundary fluxes is enabled,

(31)

which is of major importance. These fluxes contribute to the calculation of the melt rate. The section Dependent variables is for the naming the variable, hereE.

• The first subnode is the Heat transfer in fluids 1. The parameters of the ice domain (domain 2 in Comsol) are defined here. Comsol solves the following equation for the heat transfer:

ρ cp

∂E

∂t +ρcpu· ∇E =∇ ·(k∇E) +Q+Qv h+wp, (3.4) whereρ, c_p and kcan be set in the node window and are the density, the heat capacity at constant pressure and the thermal conductivity respectively of the fluid. The terms w_pQ_{v h} are Comsol interior source terms and are set to zero. The term Q is given by equ. 2.15. The equation has to be modified in order to resemble the enthalpy formulation. Therefore, the conductivity of cold ice (K_csee Table 3.3 (K_c =k_c/c_p)) is inserted in the place of the thermal conductivity, for the heat capacity is one and for the rho is set to be the density of ice (rho_ice). The model will simulate the change of enthalpy according to the flow. Hence, the heat transfer is coupled with the laminar flow (described below) under the rubricModel inputs. The corresponding laminar flow 1 is chosen as input variables for the velocity field u and the absolute pressurep.

• In the above setting the implementation of the heat transfer for the ice domain is described. Corresponding the heat transfer is implemented into the water domain.

For that reason, a secondHeat transfer in fluids (subbranch) is added to the physics branch. This one is only valid for the water domain (domain 1 in Comsol). Here the same equ. 3.4 is solved for the water domain with all interior source terms set to zero. Similar to the above described the conductivity of water (K_water), the density of water (rho_water) and the thermal conductivity of water equal one are introduced into model. Like above, the water flow is linked to the heat transfer model via the rubric Model inputs and here chose the laminar flow 2.

• Numerical models need a initial value as we solve an initial- and boundary value problem. This value should be a good guess, representing the logical physical field corresponding to the later values. Thus an incorrect choice can lead to non-convergence of a model. Here, in the subnode Initial values the enthalpy E_{star t} is chosen to be the starting value for the enthalpy calculation in the stationary solver.

Boundaries The boundary conditions which are described already in Chapter 2 are here implemented as followed in the model. It needs three boundary conditions: the ice surface (boundary ice/atmosphere), the water/ice base (ice or water/lithosphere) and the sides (ice/ice).

(32)

• Ice surface: The temperature is defined at the top of the model, which in this case the enthalpy, to be a constant value T_top, transformed into E_start (see Table 3.3). This resembles a Dirichlet-boundary condition and is implemented via the subnode Temperature.

• Base: The boundary between water domain and lithosphere is described as a Neumann boundary condition. The geothermal heat flux (q_geo) is prescribed here. This is implemented via the Comsol subnode Heat flux.

• Sides: The model is supposed to represent an interior of an ice sheet. Therefore, the sides should represent as if the model would be surrounded by an infinite number of itself. This would be accomplished by periodic boundary condition. Hence, everything, which flows out of the model, flows in on the opposite side. This type of boundary condition leads unfortunately to non-convergence of the model. Therefore, a thermal insulation boundary condition is implemented, which is a default setting by Comsol.

The heat flux is set to zero at the sides, so that −n·(−k∇E) = 0.

• Internal Boundary: The boundary between ice and water experiences a switch during the calculation. During the period that liquid is equal to one, the boundary will be treated as a von Neumann boundary. The fluxes on both side are equal q_i =q_w. As soon as liquid increases to one the boundary condition flips to a Dirichlet condition and prescribes the enthalpy to be at pressure melting pointE =Epmp. This switch is implemented via the subnode Temperature. E_pmp is set in the field for temperature.

In order to introduce the switch the equation settings have to be manipulated. Here the constrain setting is multiplied by a combination of parameter as shown below in order to assure the Dirichlet condition to be full filled during an increased water thickness and limited only to the region where Epmp is already reached.

bc =

( 0 , i f liquid<1 &E < Epmp

1 , i f liquid= 1 &E ≥E_pmp ∨ liquid = 1 &E < E_pmp (3.5) 3.1.6. Laminar flow 1 + 2

The flow of the model is introduced via the physics node Laminar flow. As mentioned in the introduction, this thesis includes as well turbulent water flow, hence the title of seams to be a contradiction. This obstacle can be overcome, since theLaminar flow includes the Navier-Stokes equation and is therefore able to model turbulent fluxes. For better handling the Navier-Stokes equation is implemented via two separated Comsol physics nodes, the so calledLaminar flow 1 and2 (corresponding to the ice and the water domain, respectively).

The Navier-Stokes equation is able to model the high viscous flow of ice and turbulent flow of water. In the ice domain the high viscosity of ice (≈ 10¹⁴P a s), hence a low Reynolds number, will minimize the effect of turbulent part in Navier-Stokes equation (v∇v ≈0)

(33)

. Whereas, the low viscous water flow in the water domain will allow turbulent fluxes. It solves for the variablesu,v,p (v_xⁱ, v_yⁱ, pⁱ), the velocity and the pressure of the ice flow and u2,v2,p2 (v_x^w, v_y^w, p^w), the velocity and the pressure of the water flow. In order to avoid confusion, as mentioned above, the x and y subscript denotes the horizontal and vertical direction. All other variables and parameters necessary for the Laminar flow are set in Table 3.1 and Table 3.4 In the following the two Laminar flow nodes will be described.

3.1.6.1. Laminar flow 1, ice flow Fluid properties

• Laminar flow 1, ice flow: Similar to the physics nodeHeat transfer the first node is the top branch Laminar flow, ice flow. Here the Navier-Stokes equations is displayed, the PDE which is solved for. Further, the ice domain (domain 2) is chosen to be the region, where the PDE is valid. Under Physical model rubric the flow type is set to be incompressible since ice is assumed to be incompressible. Like in the Heat transfer branch the model needs a specify the discretisation of the PDE. In turbulent flow modelling it is convenient to discretize the velocity part with quadratic elements (P2) and the pressure with linear elements (P1). This arises from the Babuška- Brezzi condition, which states that for numerical stability the basis functions for velocity have to an order higher than the basis functions for pressure. If the same basis functions are applied in order to decrease computational cost, the usage of a stabilization technique is necessary. In Comsol the Streamline diffusion (Galerkin Last Square) is activated in order to achieve the necessary stabilization.

• Fluid properties 1: This subnode is used for the implementation of the ice density (rho_ice) and the dynamic ice viscosity (nu_ice) of ice (see Table 3.1 and 3.4).

The density of pure ice is 917 kg m³, which is a bit more than our value, due to the fact that our lower value accounts for possible gas bubble intrusions. Moreover, a small value of 10⁻²⁵ and 10⁻³⁰ is added to the dynamic viscosity of ice η_{i c e} and the effective deformation rate, respectively. These values are added to prevent that the initial effective deformation is no zero and thus the initial dynamic viscosity to be no zero. A starting dynamic viscosity would lead to unrealistic velocities and a non-convergence of the model.

• Initial values: Here, similar toHeat transfer, a initial velocity of zero and same for the pressure is implemented into the model via theInitial valuesubnode for the stationary solver.

• The Volume force, here the gravity, is given by the vectorF = (0, −ρi c e·g)^> ,with g being the gravity acceleration and therefore, pointing in downward (in negative y

(34)

direction see Fig. 2.2). It is implemented via the subnode Volume forces into the physics branch.

Table 3.4.: Laminar flow 1 + 2 variables and parameters

name expression unit description

A 10^-16 Pa⁻³ a⁻¹ Rate factor

of ice nu_ice 0.5*A^-(1/n)*(d+10^-30)^((1-n)/n) Pa s dynamic

viscosity of ice d sqrt(0.5*ux^2+0.25*(uy+vx)^2+10^-25) s⁻¹ effective

deformation rate

nu_water 1.7e-3 Pa s viscosity of

water

u_input 0.0001. . .10 m s⁻¹ input

velocity for water domain

beta2 50 and ∞ Pa a m⁻¹ friction

coefficient slip -beta2*(u*ny+v*nx)*test(u*ny+v*nx) slip value

which

implemented into slip equation

Boundary conditions

• The ice surface to the atmosphere is tension free, when neglecting wind stress. This boundary condition is implemented via the subnode Free Boundary. The stresses are set to zero, symbolizing the tension free surface.

• The boundary between ice and water: Since this boundary applies for the water and the ice domain this boundary has two conditions at the same time depending on the domain. Therefore, the boundary has to be discussed in two steps. Here, it is focused at the ice side, the water side will be described in the section 3.1.6.2, the water flow implementation. When the ice is frozen to the ground, the basal velocity is zero. In this case a non - slip boundary condition applies to the ice. This can be added to the model by choosing a Wall as subnode and defining that Wall with a no slip condition. If the water domain turns to liquid the boundary condition has to change to a slip boundary condition. In order to achieve such a switch, a slip wall condition is implemented as if the ice could slide freely over the water domain. In

(35)

the same time the friction parameter is controlled, which regulates the sliding. This is implemented via the subnode Wall, and changing to Slip in the drop down menu.

Further, in the equation settings of the wall the weak expression has to be edited and the friction depending on liquid is inserted into the equation settings. For the different experiments the friction coefficient (β²) is adapted to the necessary settings and is described in the experiments setups.

• The Sides of the ice domain are chosen to represent an interior of an ice sheet.

Therefore, periodic boundary conditions are applied on the sides. So that everything which flow out of the model flow into the model on the over sides. Hence, it seams, as the model would be surrounded by an infinite number of itself.

3.1.6.2. Laminar flow 2, water flow Fluid properties

• Laminar flow 2, water flow: Here in the top branch, similar to the above described physic nodes, the water domain (domain 1) is set to be the region, where the PDE is valid. Further the same proceeder is applied for the discretization to set the same order of basis function for velocity and pressure and enable the streamline diffusion.

The variables, which are solving for, are u2,v2,p2 in order to distinct them form the velocity and pressure variables in the ice domain.

• Fluid properties: Here the liquid switch is implemented as well. This leads to a switch from the density and viscosity of ice to density and viscosity of water (see Table 3.1). In the model a constant density and viscosity of water is applied for simplicity reason, since variance is small in the modelled temperature regime.

• Volume force: Corresponding to the ice domain the gravity as volume force would be F = (0,−g·ρ_{w ater}). Due numerical instabilities the volume force is set to zero and therefore, the input velocity only controls the water flow.

Boundary conditions

• Initial values: The boundary conditions applied on the water domain imply an input velocity and an output pressure. The value for the input velocity (u_input Table 3.4) is the initial value for the water domain.

• surface ice water and water lithosphere: The upper and lower wall surrounding the water domain have a slip boundary condition. This is implemented via adding the Wall subnode to the branch and switching to the wall condition to slip.

(36)

• The Sides have two different boundary condition, which require each other. On the left side of the inlet boundary condition is valid and on the right hand side operates the outlet boundary condition.

– The Inlet boundary conditions is added to the branch and enables to chose a inlet velocity or pressure at the boundary. Here, the inlet velocity is set to u_input parameter (see Table 3.4) in the field for the normal velocity.

– Outlet: In order that the inflowing water can exit the water domain, the outlet boundary is applied on the right side of the model. Here, the pressure at the outlet is put zero. This brings numerical stability into the model since defining the same outflow velocity brings a non-convergence. The implementation is done via the subnode Outflow.

3.1.7. Solver

In order to find a solution to the above implemented PDEs, Comsol Multiphysics ^c offers a variety of solvers. They approximated solution of the FEM problem is to a linearized problem and solve the equations. Since this thesis wants to study the evolution of the melt dynamics, we are interested in a prognostic solution (time dependent problems).

A diagnostic (stationary) solver is used in order to achieve a good first guess for the temperature and ice velocity field. The two solver are linked to each other like illustrated in Fig. 3.4. In Comsol the solvers are implemented in the Model tree under the branch Solver.

Diagnostic solver (Study step 1): The purpose of the stationary solver is to find initial values for the velocity and pressure variables (u,v, p and u2, v2, p2) and the enthalpy E. In the Stationary Solver branch is set to solve only the above variables. These variables are solved with the linear solver Direct, which offers the solver types Mumps, Paradiso or Spooles. In this thesis the Mumps type is used as direct solver. This solver can be linked either with theFully Coupled Solver or the Segregated Solver. The Fully Coupled Solver solves all variables at once, whereas theSegregated Solver solves groups of variables successively. It is called quasi Newton solver. Normally the Segregated Solver is not as memory intensive as the Fully Coupled Solver, as the later builds the Jacobian matrix for all unknowns at once, which requires computing power. The segregated solver builds one matrices for one or several unknowns. TheStationary Solver is linked to the Fully Coupled Solver as it converges faster. The solution of the stationary solver is saved in the Solver branch. It is later used as it serves as initial value for the time dependent solver.

Prognostic solver (Study step 2): The second step of the solution process is the time dependent solver. This solver technique enables us to keep track of a system’s evolution over