Glacial Isostatic Adjustment (GIA) Solution

Physical basis

The ISSM/GIA model assumes that the ice sheet rests on top of the solid Earth, which is considered to be a simple two-layered incompressible continuum with upper elastic lithosphere floating on the viscoelastic (Maxwell material) mantle half-space. Coordinate transformations allow simple axisymmetric solutions for the deformation of pre-stressed solid Earth (subject to a normal surface traction of ice/ocean) to retrieve semi-analytical solutions of vertical displacement at the lithosphere surface.

Vertical surface displacement

Vertical displacement at the lithosphere surface (i.e., ice/ocean-bedrock interface), , is the most relevant field variable for GIA assessment. For brevity, hereinafter, this is referred to as the GIA solution. Semi-analytical GIA solution is given by [Ivins1999]:

where:

• is the radial distance from the center of the cylindrical disc load
• is the evaluation time
• is the Hankel transform variable of (or wavenumber)
• is the radius of the cylindrical disc load
• is the shear modulus of elasticity of lithosphere
• is the lithosphere density
• is the vertical component of the gravity vector
• is the -th order Bessel function of the first kind
• accounts for the integrated influence of ice loading history (cf. Figure 1) at the evaluation time . (Note that is the -th order Hankel transform of function .)

Schematic of evolution of piecewise continuous load height, $h_0$, with $J$ linear segments (from [Ivins1999]). For $j$-th segment, we can compute $m_j$ and $b_j$ (cf. Eqs. 3–4) based on the ice load at time $t_{j-1}$ and $t_j$. At $t_j$, for example, ice load at the lithosphere surface is given by $\rho_0 g h_{0j}$, where $\rho_0$ is the ice density. Assuming , the term can be written as follows:

for :

and for (i.e. the last load segment):

where:

• is the slope of the linear -th load segment
• is the -intercept of the linear -th load segment
• is the inverse decay time
• is the amplitude factor

For , the inverse decay times are given by:

and the amplitude factors by:

Parameters appearing in Eqs. (5) and (6) are defined as follows:

where:

• is the lithosphere thickness
• is the Maxwell relaxation time
• is the effective viscosity of mantle
• is the shear modulus of elasticity of mantle

• parameters with primes, e.g. , are dimensionless (listed in Table 1) with the following dimensionless parameters:

• 
• 
• 
• 
• 

• where:

• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 

The following set of non-dimensionlized parameters are defined, as needed to express dimensionless terms listed in Table 2:

where:

• is the mantle density

Numerical implementation

In the Cartesian frame of ISSM, we treat the size of ice load as the property of mesh element and compute the GIA solution at each node of the element [Adhikari2014]. Individual 2-D (-plane) mesh elements are defined as the equivalence of footprint (i.e., projection onto the -plane) of cylindrical disc loads, ensuring that the corresponding element and disc both share the same origin and plan-form area. The height of ice load is then assigned to each element such that the average normal tractional force on the corresponding area of bedrock is conserved. At each node within the domain, the final GIA solutions are computed by integrating the solutions due to individual disc loads, defined as the property of mesh elements.

Model parameters

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

>> md.gia

• md.gia.mantle_viscosity: mantle viscosity (in Pa s)
• md.gia.lithosphere_thickness: lithosphere thickness (in km)
• md.gia.cross_section_shape: shape of the cylindrical disc load; 1: square-edged (default) 2: elliptical

The solution will also use the following model fields:

• md.materials.lithosphere_shear_modulus: shear modulus of lithosphere (in Pa)
• md.materials.lithosphere_density: lithosphere density (in g/cm)
• md.materials.mantle_shear_modulus: shear modulus of mantle (in Pa)
• md.materials.mantle_density: mantle density (in g/cm)
• md.timestepping.start_time: GIA evaluation time (in yr)
• md.timestepping.final_time: in Figure 1 (in yr).
• md.geometry.thickness: ice loading history in the matrix form; the -th row, for example, should be defined as (cf. Figure 1).

ISSM Configuration

To activate the GIA model, add the following in the configuration script and compile ISSM:

--with-math77-dir="${ISSM_DIR}/externalpackages/math77/install"

Running a simulation

To run a simulation, use the following command:

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

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

References