Glacial Isostatic Adjustment (GIA)
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 twolayered incompressible continuum with upper elastic
lithosphere floating on the viscoelastic (Maxwell material) mantle halfspace. Coordinate
transformations allow simple axisymmetric solutions for the deformation of prestressed
solid Earth (subject to a normal surface traction of ice/ocean) to retrieve
semianalytical solutions of vertical displacement at the lithosphere surface.
Vertical surface displacement
Vertical displacement at the lithosphere surface (i.e., ice/oceanbedrock interface),
, is the most relevant
field variable for GIA assessment. For brevity, hereinafter, this is referred to as the GIA
solution. Semianalytical GIA solution is given by [Ivins1999]:
where:
 is the radial distance from the centre 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, , with linear segments (from [Ivins1999]). For th segment, we can compute and (cf. Eqs. 34) based on the ice load at time and . At , for example, ice load at the lithosphere surface is given by , where 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 nondimensionlized 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 2D
(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 planform 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 typing:
>> 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: squareedged (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:
withmath77dir="$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

S. Adhikari, E. Ivins, E. Larour, H. Seroussi, M. Morlighem, and S. Nowicki.
Future antarctic bed topography and its implications for ice sheet
dynamic.
Solid Earth, 5(1):569584, 2014.

Erik R. Ivins and Thomas S. James.
Simple models for late Holocene and presentday Patagonian
glacier fluctuations and predictions of a geodetically detectable isostatic
response.
Geophys. J. Int., 138(3):601624, 1999.