Controlling basal strength¶
Contents
When using option stress_balance ssa+sia
, the SIA+SSA hybrid stress balance, a
model for basal resistance is required. This model for basal resistance is based, at least
conceptually, on the hypothesis that the ice sheet is underlain by a layer of till
[103]. The user can control the parts of this model:
the socalled sliding law, typically a power law, which relates the ice base (sliding) velocity to the basal shear stress, and which has a coefficient which is or has the units of a yield stress,
the model relating the effective pressure on the till layer to the yield stress of that layer, and
the model for relating the amount of water stored in the till to the effective pressure on the till.
This subsection explains the relevant options.
The primary example of stress_balance ssa+sia
usage is in section Getting started: a Greenland ice sheet example of
this Manual, but the option is also used in sections MISMIP,
MISMIP3d, and Example: A regional model of the Jakobshavn outlet glacier in Greenland.
In PISM the key coefficient in the sliding is always denoted as yield stress \(\tau_c\),
which is tauc
in PISM output files. This parameter represents the strength of the
aggregate material at the base of an ice sheet, a poorlyobserved mixture of pressurized
liquid water, ice, granular till, and bedrock bumps. The yield stress concept also extends
to the power law form, and thus most standard sliding laws can be chosen by user options
(below). One reason that the yield stress is a useful parameter is that it can be
compared, when looking at PISM output files, to the driving stress (taud_mag
in
PISM output files). Specifically, where tauc
\(<\) taud_mag
you are likely to
see sliding if option stress_balance ssa+sia
is used.
A historical note on modeling basal sliding is in order. Sliding can be added directly to a SIA stress balance model by making the sliding velocity a local function of the basal value of the driving stress. Such an SIA sliding mechanism appears in ISMIPHEINO [104] and in EISMINT II experiment H [22], among other places. This kind of sliding is not recommended, as it does not make sense to regard the driving stress as the local generator of flow if the bed is not holding all of that stress [29], [48]. Within PISM, for historical reasons, there is an implementation of SIAbased sliding only for verification test E; see section Verification. PISM does not support this SIAbased sliding mode in other contexts.
Choosing the sliding law¶
In PISM the sliding law can be chosen to be a purelyplastic (Coulomb) model, namely,
Equation (15) says that the (vector) basal shear stress \(\boldsymbol{\tau}_b\) is at most the yield stress \(\tau_c\), and that only when the shear stress reaches the yield value can there be sliding. The sliding law can, however, also be chosen to be the power law
where \(u_{\text{threshold}}\) is a parameter with units of velocity (see below). Condition (15) is studied in [45] and [105] in particular, while power laws for sliding are common across the glaciological literature (e.g. see [80], [56]). Notice that the coefficient \(\tau_c\) in (16) has units of stress, regardless of the power \(q\).
In both of the above equations (15) and (16) we call
\(\tau_c\) the yield stress. It corresponds to the variable tauc
in PISM output files.
We call the power law (16) a “pseudoplastic” law with power \(q\) and
threshold velocity \(u_{\text{threshold}}\). At the threshold velocity the basal shear
stress \(\boldsymbol{\tau}_b\) has exact magnitude \(\tau_c\). In equation
(16), \(q\) is the power controlled by pseudo_plastic_q
, and the
threshold velocity \(u_{\text{threshold}}\) is controlled by pseudo_plastic_uthreshold
.
The plastic model (15) is the \(q=0\) case of (16).
See Table 13 for options controlling the choice of sliding law. The
purely plastic case is the default; just use stress_balance ssa+sia
to turn it on.
(Or use stress_balance ssa
if a model with no vertical shear is desired.)
Warning
Options pseudo_plastic_q
and pseudo_plastic_uthreshold
have no effect if
pseudo_plastic
is not set.
Option 
Description 


Enables the pseudoplastic power law model. If this is not set the sliding law is
purelyplastic, so 

Set the value of \(\epsilon\) regularization of the plastic law, in the formula
\(\boldsymbol{\tau}_b =  \tau_c \mathbf{u}/\sqrt{\mathbf{u}^2 + \epsilon^2}\). The
default is \(0.01\) m/a. This parameter is inactive if 

Set the exponent \(q\) in (16). The default is \(0.25\). 

Set \(u_{\text{threshold}}\) in (16). The default is \(100\) m/a. 
Equation (16) is a very flexible power law form. For example, the linear case is \(q=1\), in which case if \(\beta=\tau_c/u_{\text{threshold}}\) then the law is of the form
(The “\(\beta\)” coefficient is also called \(\beta^2\) in some sources (see [44],
for example).) If you want such a linear sliding law, and you have a value
\(\beta=\) beta
in \(\text{Pa}\,\text{s}\,\text{m}^{1}\), then use this option
combination:
pseudo_plastic \
pseudo_plastic_q 1.0 \
pseudo_plastic_uthreshold 1m/s \
yield_stress constant tauc beta
More generally, it is common in the literature to see powerlaw sliding relations in the form
where \(C\) is a constant, as for example in sections MISMIP and MISMIP3d. In that case, use this option combination:
pseudo_plastic \
pseudo_plastic_q m \
pseudo_plastic_uthreshold 1m/s \
yield_stress constant \
tauc C
Another alternative is the slip law by [106] that combines processes of hardbedded sliding (Coulomb behavior) and viscous bed deformation without required knowledge of the bed type. It is defined by
Set regularized_coulomb
to select this sliding law.
The original equation (3) in [106] uses the exponent \(q=1/p\). Otherwise, same configuration parameters can be used as in the pseudoplastic case, but they should have different values, namely \(q\)=0.2, \(u_{\text{threshold}}\)=4080 m/year`:
regularized_coulomb \
pseudo_plastic_q 0.2 \
pseudo_plastic_uthreshold 50.0
The model’s performance should be close to the pseudoplastic implementation (Eq. (16)), although there ought to be slightly more fast sliding and a less diffuse onset of sliding.
Lateral drag¶
PISM prescribes lateral drag at ice margins next to ground with elevation above the ice.
(This is relevant in outlet glaciers flowing through fjords, valley glaciers, and next to
nunataks.) Set basal_resistance
.beta_lateral_margin
to control the amount of
additional drag at these margins.
Determining the yield stress¶
Other than setting it to a constant, which only applies in some special cases, the above discussion does not determine the yield stress \(\tau_c\). As shown in Table 14, there are two schemes for determining \(\tau_c\) in a spatiallyvariable manner:
yield_stress mohr_coulomb
(the default) determines the yields stress by models of till material property (the till friction angle) and of the effective pressure on the saturated till, oryield_stress constant
allows the yield stress to be supplied as timeindependent data, read from the input file.
In normal modelling cases, variations in yield stress are part of the explanation of the
locations of ice streams [45]. The default model yield_stress
mohr_coulomb
determines these variations in time and space. The value of \(\tau_c\) is
determined in part by a subglacial hydrology model, including the modeled tillpore water
amount (\(W_{till}\), NetCDF variable tillwat
; see Subglacial hydrology), which then
determines the effective pressure \(N_{till}\) (see below). The value of \(\tau_c\) is also
determined in part by a material property field \(\phi\), the “till friction angle” (NetCDF
variable tillphi
). These quantities are related by the MohrCoulomb criterion
[80]:
Here \(c_0\) is called the “till cohesion”, whose default value in PISM is zero (see
[45], formula (2.4)) but which can be set by option till_cohesion
.
Option combination yield_stress constant tauc X
can be used to fix the yield stress
to have value \(\tau_c = X\) at all grounded locations and all times if desired. This is
unlikely to be a good modelling choice for real ice sheets.
Option 
Description 


The default. Use equation (19) to determine \(\tau_c\). Only
effective if 

Set the value of the till cohesion (\(c_{0}\)) in the plastic till model. The value is a pressure, given in Pa. 

If set, reduces the basal yield stress at groundedbelowsealevel grid points one cell away from floating ice or ocean. Specifically, it replaces the normallycomputed \(\tau_c\) from the MohrCoulomb relation, which uses the effective pressure from the modeled amount of water in the till, by the minimum value of \(\tau_c\) from MohrCoulomb, i.e. using the effective pressure corresponding to the maximum amount of tillstored water. Does not alter the reported amount of till water, nor does this mechanism affect water conservation. 

Use a constant till friction angle. The default is \(30^{\circ}\). 

Compute \(\phi\) using equation (20). 

Compute the till friction angle \(\phi\) in (19) iteratively in an equilibrium simulation using equation (21). 

Keep the current values of the till yield stress \(\tau_c\). That is, do not update
them by the default model using the stored basal melt water. Only effective if


Directly set the till yield stress \(\tau_c\), in units Pa, at all grounded locations
and all times. Only effective if used with 
Till friction angle heuristic¶
We find that an effective, though heuristic, way to determine \(\phi\) in (19) is to make it a function of bed elevation [66], [3], [37]. This heuristic is motivated by hypothesis that basal material with a marine history should be weak [4]. PISM has a mechanism setting \(\phi\) to be a piecewise linear function of bed elevation. The option is
topg_to_phi phimin,phimax,bmin,bmax
Thus the user supplies 4 parameters: \(\phimin\), \(\phimax\), \(\bmin\), \(\bmax\), where \(b\) stands for the bed elevation. To explain these, we define \(M = (\phimax  \phimin) / (\bmax  \bmin)\). Then
It is worth noting that an earth deformation model (see section Earth deformation models) changes
\(b(x,y)\) (NetCDF variable topg
) used in (20), so that a sequence
of runs such as
pismr i foo.nc bed_def lc stress_balance ssa+sia \
topg_to_phi 10,30,50,0 ... o bar.nc
pismr i bar.nc bed_def lc stress_balance ssa+sia \
topg_to_phi 10,30,50,0 ... o baz.nc
will use different tillphi
fields in the first and second runs. PISM will print a
warning during initialization of the second run:
* Initializing the default basal yield stress model...
option topg_to_phi seen; creating tillphi map from bed elev ...
PISM WARNING: topg_to_phi computation will override the 'tillphi' field
present in the input file 'bar.nc'!
Omitting topg_to_phi
in the second run would make PISM continue with the
same tillphi
field which was set in the first run.
Parameters
Prefix: basal_yield_stress.mohr_coulomb.topg_to_phi.
enabled
(no) If the optiontopg_to_phi
is set then this will be set to “yes”, and thenMohrCoulombYieldStress
will initialize thetillphi
field using a piecewise linear function of depth described by four parameters.phi_max
(15 degrees) upper value of the till friction angle; see the implementation of MohrCoulombYieldStressphi_min
(5 degrees) lower value of the till friction angle; see the implementation of MohrCoulombYieldStresstopg_max
(1000 meters) the elevation at which the upper value of the till friction angle is used; see the implementation of MohrCoulombYieldStresstopg_min
(1000 meters) the elevation at which the lower value of the till friction angle is used; see the implementation of MohrCoulombYieldStress
Till friction angle optimization¶
Warning
This is a work in progress. Use at your own risk.
In grounded areas the distribution of till friction angle \(\phi\) (see (19)) can be iteratively optimized in a forward equilibrium simulation. 1
The iteration starts from a \(\phi\) distribution set using plastic_phi
,
topg_to_phi
, or read from an input file (variable tillphi
).
During each step, an adjustment \(\Delta \phi\) is added to the previous value of \(\phi\) and the result is clipped to ensure \(\phi \in [\phimin, \phimax]\):
The value of \(\Delta \phi\) is proportional to the difference between the target (usually observed present day, e.g. [107]) and modeled surface elevations. It is similarly clipped to ensure \(\Delta \phi \in [\dphimin, \dphimax]\):
Here \(C\) is the (positive) scaling factor (units: \({}^\circ / \mathrm{m}\)) set using
basal_yield_stress
.mohr_coulomb
.tillphi_opt
.dphi_scale
.
Note
The adjustment \(\Delta \phi\) is positive if the modeled surface elevation is below the reference value, and negative otherwise.
In other words, the basal resistance is increased if the ice thickness is too low and decreased otherwise.
The lower bound \(\phimin = \phi_0\) is a piecewise linear function of the bed topography \(b\):
Similarly to the till friction angle heuristic (20), we assume that “marine” sediments (below \(\bmin\)) can be much weaker than rather “continental” bedrock material (above \(\bmax\)). In sensitivity experiments we found a strong sensitivity of the Antarctic Ice Sheet’s ice volume in particular to the choice of \(\phimin\) (see [108]).
To allow ice geometry to respond to changes in the till friction angle the simulation goes on for \(\dt_{\phi}\) years between iterations.
Iterations at a particular location are considered “done” when the rate of change of the
surface elevation mismatch \(\Delta h = h_{\mathrm{observed}}  h_{\mathrm{modeled}}\)
approximated using subsequent steps falls below a threshold \(D\) set using
basal_yield_stress
.mohr_coulomb
.tillphi_opt
.dhdt_min
:
The diff_mask
diagnostic variable is set to 0 to indicate that \(\phi\) “converged”
at this location.
Parameters
Prefix: basal_yield_stress.mohr_coulomb.tillphi_opt.
dhdt_min
(0.001 meters year1) rate of change in the surface elevation mismatch \(D\) used as a convergence criteriondphi_max
(1 degrees) maximum till friction angle adjustment \(\dphimax\)dphi_scale
(0.002 degree / meters) scaling factor \(C\) used to compute the \(\Delta \phi\) adjustment, \(C\) degrees per meter of surface elevation mismatchdt
(250 365days) time step \(\dt_{\phi}\) for optimization of till friction anglefile
Name of the file containing the timeindependent variableusurf
used as target surface elevationphi0_max
(5 degrees) maximum value of the lower bound of the till friction angle, \(\phi_{0,\mathrm{max}}\)phi0_min
(2 degrees) minimum value of the lower bound of the till friction angle, \(\phi_{0,\mathrm{min}}\)phi_max
(70 degrees) upper bound of the till friction angle \(\phimax\)topg_max
(700 meters) the bed elevation \(\bmax\) above whichbasal_yield_stress
.mohr_coulomb
.tillphi_opt
.phi0_max
is usedtopg_min
(300 meters) the bed elevation \(\bmin\) below whichbasal_yield_stress
.mohr_coulomb
.tillphi_opt
.phi0_min
is used
When the domain contains a grounding line the mismatch between modeled and observed surface elevation is meaningful only if the ice is grounded in both data sets. A retreat of the grounding line would make it impossible to optimize the till friction angle in areas that are observed to contain grounded ice but are covered by water in a simulation.
To avoid this issue, we
disable the influence of the basal melt rate on geometry evolution by setting
geometry
.update
.use_basal_melt_rate
to “false”,modify the surface mass balance in grounded areas to disallow grounding line retreat by adding
no_gl_retreat
to the commandline option selecting a surface model (see Preventing grounding line retreat), andfix ice thickness where the bed elevation is below sea level and the ice (if present) is floating.
To fix ice thickness, we create a mask with ones where ice is floating or there is no ice and the bed is below sea level:
ncap2 O \
s "where(topg < 0 && thk*(910.0/1028.0) < 0  topg) thk_bc_mask=1;" \
input.nc inputwithmask.nc
Here \(910\) is the ice density (see constants
.ice
.density
), \(1028\) is the sea
water density (see constants
.sea_water
.density
), and \(0\) is the sea level
elevation.
Reported diagnostic quantities
usurf_difference
reports the mismatch between modeled and reference surface elevation fieldsdiff_mask
reports the area where the till friction angle is iteratively adjusted or where the convergence criterion is met.usurf_target
reports the reference (target) ice surface elevation in use.
Determining the effective pressure¶
When using the default option yield_stress mohr_coulomb
, the effective pressure on
the till \(N_{till}\) is determined by the modeled amount of water in the till. Lower
effective pressure means that more of the weight of the ice is carried by the pressurized
water in the till and thus the ice can slide more easily. That is, equation
(19) sets the value of \(\tau_c\) proportionately to \(N_{till}\). The amount
of water in the till is, however, a nontrivial output of the hydrology (see
Subglacial hydrology) and conservationofenergy (see Modeling conservation of energy) submodels in PISM.
Following [109], based on laboratory experiments with till extracted
from an ice stream in Antarctica, [110] propose the following
parameterization which is used in PISM. It is based on the ratio
\(s=W_{till}/W_{till}^{max}\) where \(W_{till}\) is the effective thickness of water in the
till and \(W_{till}^{max}\) (hydrology
.tillwat_max
) is the maximum amount of water
in the till (see section Subglacial hydrology):
Here \(P_o\) is the ice overburden pressure, which is determined entirely by the ice thickness and density, and the remaining parameters are listed below
Parameters
Prefix: basal_yield_stress.mohr_coulomb.
delta
.file
Name of the file containing space and timedependent variablemohr_coulomb_delta
to use instead ofbasal_yield_stress
.mohr_coulomb
.till_effective_fraction_overburden
.delta
.periodic
(no) If true, interpret forcing data as periodic in timetauc_to_phi
.file
File containing the basal yield stress field that should be used to recover the till friction angle distribution.till_cohesion
(0 Pascal) cohesion of till; = \(c_0\) in most references; note Schoof uses zero but Paterson pp 168–169 gives range 0–40 kPa; but Paterson notes that “… all the pairs \(c_0\) and \(\phi\) in the table would give a yield stress for Ice Stream B that exceeds the basal shear stress there…”till_compressibility_coefficient
(0.12) coefficient of compressiblity of till; value from [109]till_effective_fraction_overburden
(0.02) \(\delta\) in [110]; \(\delta P_o \le N_{\text{till}} \le P_o\) where \(P_o\) is overburden pressure and \(N_{\text{till}}\) is the effective pressure of the overlying ice on the saturated till; default value of \(\delta\) corresponds to Greenland and Antarctic model runstill_log_factor_transportable_water
(0.1 meters) Ifbasal_yield_stress
.add_transportable_water
is set then the water amount in the transport system is added totillwat
in determiningtauc
. Normally only the water in the till is used. This factor multiplies the logarithm in that formula.till_phi_default
(30 degrees) fill value for till friction angletill_reference_effective_pressure
(1000 Pascal) reference effective pressure \(N_0\); value from [109]till_reference_void_ratio
(0.69) void ratio at reference effective pressure \(N_0\); value from [109]
Note
While there is experimental support for the default values of \(C_c\), \(e_0\), and \(N_0\), the value of \(\delta\) should be regarded as uncertain, important, and subject to parameter studies to assess its effect.
Controlling minimum effective pressure¶
The effective pressure \(N_{till}\) above satisfies (see equation 20 in [110])
In other words, \(\delta\) controls the lower bound of the effective pressure. In addition
to setting it using a configuration parameter one can use a space and timedependent
field. Set basal_yield_stress
.mohr_coulomb
.delta
.file
to the name of the file
containing the variable mohr_coulomb_delta
(dimensionless, i.e. units of “1”).
Note
PISM restricts the time step to capture changes in \(\delta\). For example, if you provide monthly records of \(\delta\) PISM will make sure no time step spans more than one month.
PISM uses piecewise linear interpolation in time for model times between records of \(\delta\).
Footnotes
Previous  Up  Next 