Mass Transport Solution

Physical basis

Conservation of mass

The mass transport equation is derived from the depth-integrated form of the mass conservation equation and reads:

where:

• is the depth-averaged velocity vector
• is the ice thickness
• is the surface accumulation (in m/yr of ice equivalent, positive for accumulation)
• is the basal melting (in m/yr of ice equivalent, positive for melting)

For full-Stokes models, free surface equations are solved for the upper surface and the base of floating ice:

and:

where:

• is the elevation of the ice upper surface
• is the elevation of the floating ice lower surface
• are the ice velocity components at the upper surface
• are the ice velocity components at the base

Boundary conditions

Ice thickness is imposed at the inflow boundary:

For free surfaces models, both and are constrained at the inflow boundary.

Numerical implementation

Mass transport is solved using finite elements in space, and implicit finite difference in time. To stabilize the equation, a stabilization term might be added to the left hand side, for example:

where is the artificial diffusivity. We take:

There are other stabilization schemes available in ISSM: (1) Artificial Diffusion, (2) Streamline Upwinding, (3) Discontinuous Galerkin (DG), (4) Flux Corrected Transport (FCT), and (5) Streamline Upwind Petrov-Galerlin (SUPG). They can used by setting:

>> md.masstransport.stabilization = 1;

Model parameters

The parameters relevant to the mass transport solution can be displayed by running:

>> md.masstransport

• md.masstransport.spcthickness: thickness constraints (NaN means no constraint)
• md.masstransport.hydrostatic_adjustment: adjustment of ice shelves upper and lower surfaces: 'Incremental' or 'Absolute'
• md.masstransport.stabilization: 0: no stabilization, 1: artificial diffusion, 2: streamline upwinding 3: discontinuous Galerkin, 4: flux corrected transport (FCT), 5: streamline upwind Petrov-Galerkin (SUPG)
• md.masstransport.penalty_factor: offset used by penalties

• md.masstransport.vertex_pairing: pairs of vertices that are penalized (for periodic boundary conditions only)

The solution will also use the following model fields:

• md.smb.ablation_rate: surface ablation rate (in meters)
• md.smb.mass_balance: surface mass balance (in meters)
• md.initialization.vx: x component of velocity
• md.initialization.vy: y component of velocity
• md.basalforcings.groundedice_melting_rate: basal melting rate applied on grounded ice (positive if melting)
• md.basalforcings.floatingice_melting_rate: basal melting rate applied on floating ice (positive if melting)
• md.smb.mass_balance: surface mass balance (in meters/year ice equivalent)
• md.timestepping.time_step: length of time steps (in years)

Running a simulation

To run a simulation, use the following command:

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

The first argument is the model, the second is the nature of the simulation one wants to run. This will compute one time step of the mass transport equation; use the transient solution for multiple time steps.