Introduction to pyISSM
This notebook provides a first introduction to pyISSM, the Python interface to the Ice-sheet and Sea-level System Model (ISSM).
This notebook focuses on helping you understand how pyISSM is structured, how its syntax works, and how to interact with ISSM models in a Python-native way.
By the end of this notebook, you should be comfortable navigating pyISSM, inspecting model objects, and understanding where different pieces of functionality live.
Scope of this notebook
This notebook focuses on structure and syntax, not on glaciology or ice sheet model theory, or model development. If you are new to pyISSM, we strongly recommend working through this notebook before moving on.
Warning
MATLAB -> Python Differences
For those users who are familiar with the existing MATLAB API, pyISSM will feel familiar, but not identical. pyISSM preserves ISSMs core concepts (e.g. the model object and general workflows), while adopting Pythonic syntax and conventions.
Important differences include:
Adoption of keyword arguments, rather than positional argument pairs.
Function names use snake_case (underscores)
Import pyISSM
We’ll start by importing pyISSM and getting a feel for it’s structure. Instructions for installing pyISSM are provided here.
# Import pyISSM
import pyissm
# Check available submodules
import pkgutil
print(f'pyISSM contains the following submodules: {[m.name for m in pkgutil.iter_modules(pyissm.__path__)]}')
pyISSM contains the following submodules: ['analysis', 'data', 'model', 'plot', 'tools']
We can see that pyISSM contains five sub-modules.
analysiscontains tailored post-processing toolsdatacontains tailored tools for interacting with external datasets, including interpolation routines.modelcontains the central functionality, including the coreModel()class.plotconatains tailored tools for visualising ISSM models.toolscontains diverse tools for interacting with the ISSM Model Object.
This notebook focuses on the structure of the model submodule. Other modules are primarily used for parameterising, visualising, and analysing ISSM models and will be explored in other notebooks.
The ISSM Model Object (md)
pyISSM is developed around a central Model() class. Each model contains a series of sub-classes that define model components and parameters (e.g. mesh, geometry, etc.). Each model and sub-class contains detailed information that can be interrogated/explored interactively.
Across the ISSM ecosystem (within the core C++ code, in MATLAB, and here in pyISSM), ISSM models are commonly named md.
Let’s create an empty ISSM model!
# Create an empty model
md = pyissm.model.Model()
# View all model sub-classes by viewing the __repr__ of the model
md
ISSM Model Class
mesh: mesh properties
mask: defines grounded and floating elements
geometry: surface elevation, bedrock topography, ice thickness, ...
constants: physical constants
smb: surface mass balance
basalforcings: bed forcings
materials: material properties
damage: damage propagation laws
friction: basal friction / drag properties
flowequation: flow equations
timestepping: timestepping for transient models
initialization: initial guess / state
rifts: rifts properties
solidearth: solidearth inputs and settings
dsl: dynamic sea level
debug: debugging tools (valgrind, gprof
verbose: verbosity level in solve
settings: settings properties
toolkits: PETSc options for each solution
cluster: cluster parameters (number of CPUs...)
balancethickness: parameters for balancethickness solution
stressbalance: parameters for stressbalance solution
groundingline: parameters for groundingline solution
hydrology: parameters for hydrology solution
masstransport: parameters for masstransport solution
thermal: parameters for thermal solution
steadystate: parameters for steadystate solution
transient: parameters for transient solution
levelset: parameters for moving boundaries (level-set method)
calving: parameters for calving
frontalforcings: parameters for frontalforcings
esa: parameters for elastic adjustment solution
sampling: parameters for stochastic sampler
love: parameters for love solution
autodiff: automatic differentiation parameters
inversion: parameters for inverse methods
qmu: Dakota properties
amr: adaptive mesh refinement properties
outputdefinition: output definition
results: model results
radaroverlay: radar image for plot overlay
miscellaneous: miscellaneous fields
stochasticforcing: stochasticity applied to model forcings
We can access any model sub-class using standard Python dot notation. By default, interactively interrogating md will print a summary of the model, or specific model sub-class. When interrogating a model object, we see all available sub-classes (like the above example). when interrigating a model sub-class, we see all available fields/parameters within that sub-class. If a field/parameter is empty, N/A is displayed.
Let’s have a look at pre-defined parameters in the md.constants sub-class, and empty fields in the md.geometry sub-class.
md.constants
constants parameters:
g : 9.81 -- gravitational acceleration [m/s^2]
omega : 7.292e-05 -- angular velocity of Earth [rad/s]
yts : 31536000.0 -- number of seconds in a year [s/yr]
referencetemperature : 223.15 -- reference temperature used in the enthalpy model [K]
gravitational_constant : 6.67259e-11 -- Newtonian constant of gravitation [m^3/kg/s^2]
md.geometry
geometry parameters:
surface : N/A -- ice upper surface elevation [m]
thickness : N/A -- ice thickness [m]
base : N/A -- ice base elevation [m]
bed : N/A -- bed elevation [m]
hydrostatic_ratio : N/A -- hydrostatic ratio for floating ice
In the next cell, all available sub-classes are listed and commented out. Spend a few minutes interrogating each of them and familiarizing yourself with the structure of a pyISSM model. You can simply uncomment one line at a time and re-run the cell to view the fields/parameters included in each sub-class.
# md.amr
# md.autodiff
# md.balancethickness
# md.basalforcings
# md.calving
# md.cluster
# md.constants
# md.damage
# md.debris
# md.debug
# md.dsl
# md.esa
# md.flowequation
# md.friction
# md.frontalforcings
# md.geometry
# md.groundingline
# md.hydrology
# md.initialization
# md.inversion
# md.levelset
# md.love
# md.mask
# md.masstransport
# md.materials
# md.mesh
# md.miscellaneous
# md.outputdefinition
# md.private
# md.qmu
# md.radaroverlay
# md.results
# md.rifts,
# md.sampling
# md.settings
# md.smb
# md.solidearth
# md.steadystate
# md.stochasticforcing
# md.stressbalance
# md.thermal
# md.timestepping
# md.toolkits
# md.transient
# md.verbose
Alternative parameterisations
Within ISSM, there are multiple ways to represent various processes, depending on the use case. For example, modellers may wish to represent basal friction using the Budd friction law, or the Schoof friction law. Or, modellers may wish to use the Positive Degree Day (PDD) surface mass balance (SMB) approach, rather than a prescribed SMB. All of this flexibility is embedded in pyISSM and accessible via different model classes.
Let’s take a look at a few different friction classes, accessible via pyissm.model.classes.friction....
# By default, md.friction uses the default friction law, which is a Budd friction law
md.friction
Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b,
(effective stress Neff = rho_ice * g * thickness + rho_water * g * base, r = q / p and s = 1 / p)
coefficient : N/A -- friction coefficient [SI]
p : N/A -- p exponent
q : N/A -- q exponent
coupling : 0 -- Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)
linearize : 0 -- 0: not linearized, 1: interpolated linearly, 2: constant per element (default is 0)
effective_pressure : N/A -- Effective Pressure for the forcing if not coupled [Pa]
effective_pressure_l...: 0 -- Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)
# To change the friction law to a Weertman friction law, we can do:
md.friction = pyissm.model.classes.friction.weertman()
md.friction
Weertman sliding law parameters: Sigma_b = C^(- 1 / m) * |u_b|^(1 / m - 1) * u_b
C : N/A -- friction coefficient [SI]
m : N/A -- m exponent
linearize : 0 -- 0: not linearized, 1: interpolated linearly, 2: constant per element (default is 0)
The same appraoch allows users to adjust the SMB parameterisation (as well as other model sub-classes). Let’s look at a few different SMB options, accessible via pyissm.model.classes.smb....
# By default, md.smb uses the default surface mass balance model which uses a prescribed smb field
md.smb
# To change the smb model to a positive degree day smb model, we can do:
md.smb = pyissm.model.classes.smb.pdd()
md.smb
surface forcings parameters:
isdelta18o : 0 -- is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)
ismungsm : 0 -- is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)
issetpddfac : 0 -- is user passing in defined pdd factors (0 or 1, default is 0)
desfac : 0.5 -- desertification elevation factor (between 0 and 1, default is 0.5) [m]
s0p : N/A -- should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]
s0t : N/A -- should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]
rlaps : 6.5 -- present day lapse rate [degree/km]
rlapslgm : 6.5 -- LGM lapse rate [degree/km]
monthlytemperatures : N/A -- monthly surface temperatures [K], required if pdd is activated and delta18o not activated
precipitation : N/A -- monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated
delta18o : 0 -- delta18o [per mil], required if pdd is activated and delta18o activated
delta18o_surface : N/A -- surface elevation of the delta18o site, required if pdd is activated and delta18o activated
temperatures_presentday: N/A -- monthly present day surface temperatures [K], required if delta18o/mungsm is activated
temperatures_lgm : N/A -- monthly LGM surface temperatures [K], required if delta18o or mungsm is activated
precipitations_prese...: N/A -- monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated
precipitations_lgm : N/A -- monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated
Tdiff : N/A -- time interpolation parameter for temperature, 1D(year), required if mungsm is activated
sealev : N/A -- sea level [m], 1D(year), required if mungsm is activated
temperatures_presentday: N/A -- monthly present day surface temperatures [K], required if delta18o/mungsm is activated
temperatures_lgm : N/A -- monthly LGM surface temperatures [K], required if delta18o or mungsm is activated
precipitations_prese...: N/A -- monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated
precipitations_lgm : N/A -- monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated
Pfac : N/A -- time interpolation parameter for precipitation, 1D(year), required if mungsm is activated
Tdiff : N/A -- time interpolation parameter for temperature, 1D(year), required if mungsm is activated
sealev : N/A -- sea level [m], 1D(year), required if mungsm is activated
steps_per_step : 1 -- number of smb steps per time step
averaging : 0 -- averaging methods from short to long steps
0: Arithmetic (default)
1: Geometric
2: Harmonic
requested_outputs : ['default',] -- additional outputs requested