Thermal Solution

Physical basis

Thermal state

The heat transport equation is derived from the balance equation of internal energy $Equation 1$ combined with Fourier’s law of heat transfer and reads:

$Equation 2$

where radiative sources have been neglected, and:

$Equation 3$ is the velocity vector
$Equation 4$ is the strain rate tensor
$Equation 5$ is the internal energy density
$Equation 6$ is the specific heat conductivity, which can depend on the heat density
$Equation 7$ is the Cauchy stress tensor.

For constant heat conductivity and heat capacity $Equation 8$ , the previous equation becomes:

$Equation 9$

Boundary conditions

Dirichlet boundary conditions should be applied at the ice surface:

$Equation 10$

and Neumann boundary conditions at the ice base:

$Equation 11$

where:

$Equation 12$ is the elevation of the ice upper surface
$Equation 13$ is the elevation of the floating ice lower surface When using the enthalpy formulation, the basal boundary condition scheme from [Aschwanden2012], figure 5 is used instead of the previous equation.

NOTE: For regional model, make sure to set a Dirichlet condition on the inflow boundary (throughout the ice column) to avoid advection of noise.

Numerical implementation

The heat equation is solved using linear finite elements in space, and implicit finite difference in time (time stepping should satisfy the CFL condition). To stabilize the equation, we either add an isotropic artificial diffusion to the left hand side:

$Equation 14$

where $Equation 15$ is the artificial diffusivity. We take:

$Equation 16$

or rely on the Streamline upwind/Petrov-Galerkin formulation (SUPG) from [Franca2006].

Model parameters

The parameters relevant to the heat equation solution can be displayed by running:

>> md.thermal

md.thermal.spctemperature: temperature constraints (NaN means no constraint)
md.thermal.stabilization: type of stabilization (0: no stabilization; 1: artificial diffusion; 2: Streamline-Upwind Petrov-Galerkin)
md.thermal.maxiter: maximum number of iterations for thermal solver
md.thermal.penalty_lock: stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)
md.thermal.penalty_threshold: threshold to declare convergence of thermal solution (default is 0)
md.thermal.penalty_factor: offset used by penalties(default is 3):

$Equation 17$

md.thermal.isenthalpy: are we using the enthalpy formulation (Aschwanden et al., 2012)? (0: no, 1: yes)
md.thermal.isdynamicbasalspc: are we allowing changing basal boundary conditions for transient runs?
md.thermal.requested_outputs: specify further requested outputs here.

The solution will also use the following model fields:

md.initialization.temperature: temperature field (in K)
md.initialization.waterfraction: water fraction in ice (between $Equation 19$ and $Equation 18$ )
md.basalforcings.geothermalflux: geothermal heat flux (in $Equation 20$ )
md.basalforcings.meltingrate: basal melting rate (in m/yr w.e.)
md.timestepping.time_step: length of time steps (in yrs)

Running a simulation

To run a simulation solving only the thermal state, use the following command:

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

This will compute one time step only of the thermal equation; use the transient solution for multiple time steps.

To run a thermal steady-state simulation (i.e. $Equation 21$ ), you need to first set the time stepping as 0:

>> md.timestepping.time_step = 0
>> md = solve(md, 'Thermal');