Stress Balance Solution

Physical basis

Conservation of linear momentum

The conservation of momentum reads:

$Equation 1$

where:

$Equation 2$ is the ice density
$Equation 3$ is the velocity vector
$Equation 4$ is the Cauchy stress tensor
$Equation 5$ is a body force Now if we assume that:
The ice motion is a Stokes flow (acceleration negligible)
The only body force is due to gravity (Coriolis effect negligible) The equation of momentum conservation is reduced to:

$Equation 6$

Conservation of angular momentum

For a non-polar material body, the balance of angular momentum imposes the stress tensor to be symmetrical:

$Equation 7$

Ice constitutive equations

Ice is treated as a purely viscous incompressible material [Cuffey2010]. Its constitutive equation therefore only involves the deviatoric stress and the strain rate tensor:

$Equation 8$

where:

$Equation 10$ is the deviatoric stress tensor ( $Equation 9$ )
$Equation 11$ is the ice effective viscosity
$Equation 12$ is the strain rate tensor Ice is a non-Newtonian fluid, its viscosity follows the generalized Glen’s flow law [Glen1955]:

$Equation 13$

where:

$Equation 14$ is the ice hardness or rigidity
$Equation 15$ is Glen’s flow law exponent, generally taken as equal to 3
$Equation 16$ is the effective strain rate The effective strain rate is defined as:

$Equation 17$

where $Equation 18$ is the Frobenius norm.

Full-Stokes (FS) field equations

Without any further approximation, the previous system of equations are called the Full-Stokes model.

Higher-Order (HO) field equations

We make two assumptions:

Bridging effects are neglected
Horizontal gradient of vertical velocities are neglected compared to vertical gradients of horizontal velocities With these two assumptions, the Full-Stokes equations are reduced to a system of 2 equations with 2 unknowns [Blatter1995, Pattyn2003]:

$Equation 19$

with:

$Equation 20$

Shelfy-Stream Approximation (SSA) field equations

We make the following assumption:

Vertical shear is negligible With this assumption, we have a system of 2 equations with 2 unknowns in the horizontal plane [Morland1987a, MacAyeal1989]:

$Equation 21$

with:

$Equation 22$

where:

$Equation 23$ is the depth-averaged viscosity
$Equation 24$ is the ice thickness
$Equation 25$ is the basal friction coefficient

Boundary conditions

At the surface of the ice sheet, $Equation 28$ , we assume a stress-free boundary condition. A viscous friction law is applied at the base of the ice sheet, $Equation 27$ , and water pressure is applied at the ice/water interface $Equation 26$ . For FS, these boundary conditions are:

$Equation 29$

where

$Equation 30$ is the outward-pointing unit normal vector
$Equation 31$ is the water density
$Equation 32$ is the vertical coordinate with respect to sea level

For HO, these boundary conditions become:

$Equation 33$

where $Equation 34$ .

For SSA, these boundary conditions are:

$Equation 35$

$Equation 36$

Model parameters

The parameters relevant to the stress balance solution can be displayed by typing:

>> md.stressbalance

md.stressbalance.restol: mechanical equilibrium residue convergence criterion
md.stressbalance.reltol: velocity relative convergence criterion, (NaN if not applied)
md.stressbalance.abstol: velocity absolute convergence criterion, (NaN if not applied)
md.stressbalance.maxiter: maximum number of nonlinear iterations (default is 100)
md.stressbalance.spcvx: x-axis velocity constraint (NaN means no constraint)
md.stressbalance.spcvy: y-axis velocity constraint (NaN means no constraint)
md.stressbalance.spcvz: z-axis velocity constraint (NaN means no constraint)
md.stressbalance.rift_penalty_threshold: threshold for instability of mechanical constraints
md.stressbalance.rift_penalty_lock: number of iterations before rift penalties are locked
md.stressbalance.penalty_factor: offset used by penalties:

$Equation 37$

md.stressbalance.vertex_pairing: pairs of vertices that are penalized
md.stressbalance.shelf_dampening: use dampening for floating ice? Only for Stokes model
md.stressbalance.referential: local referential
md.stressbalance.requested_outputs: additional outputs requested

The solution will also use the following model fields:

md.flowequations: FS, HO or SSA
md.materials: material parameters
md.initialization.vx: x component of velocity (used as an initial guess)
md.initialization.vy: y component of velocity (used as an initial guess)
md.initialization.vz: y component of velocity (used as an initial guess)

Running a simulation

To run a simulation, use the following command:

>> md = solve(md, 'Stressbalance');

The first argument is the model, the second is the nature of the simulation one wants to run.