Controlling basal strength ¶

Contents

Controlling basal strength

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 [18]. The user can control the parts of this model:

the so-called 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 poorly-observed 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 ISMIP-HEINO [86] and in EISMINT II experiment H [12], 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 [10], [47]. Within PISM, for historical reasons, there is an implementation of SIA-based sliding only for verification test E; see section Verification. PISM does not support this SIA-based sliding mode in other contexts.

Choosing the sliding law ¶

In PISM the sliding law can be chosen to be a purely-plastic (Coulomb) model, namely,

(15)¶\[|\boldsymbol{\tau}_b| \le \tau_c \quad \text{and} \quad \boldsymbol{\tau}_b = - \tau_c \frac{\mathbf{u}}{|\mathbf{u}|} \quad\text{if and only if}\quad |\mathbf{u}| > 0.\]

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

(16)¶\[\boldsymbol{\tau}_b = - \tau_c \frac{\mathbf{u}}{u_{\text{threshold}}^q |\mathbf{u}|^{1-q}},\]

where \(u_{\text{threshold}}\) is a parameter with units of velocity (see below). Condition (15) is studied in [21] and [87] in particular, while power laws for sliding are common across the glaciological literature (e.g. see [34], [55]). 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 “pseudo-plastic” 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.

Table 13 Sliding law command-line options¶
Option	Description
`-pseudo_plastic`	Enables the pseudo-plastic power law model. If this is not set the sliding law is purely-plastic, so `pseudo_plastic_q` and `pseudo_plastic_uthreshold` are inactive.
`-plastic_reg` (m/a)	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 `-pseudo_plastic` is set.
`-pseudo_plastic_q`	Set the exponent \(q\) in (16). The default is \(0.25\).
`-pseudo_plastic_uthreshold` (m/a)	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

(17)¶\[\boldsymbol{\tau}_b = - \beta \mathbf{u}\]

(The “\(\beta\)” coefficient is also called \(\beta^2\) in some sources (see [14], 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 power-law sliding relations in the form

\[\boldsymbol{\tau}_b = - C |\mathbf{u}|^{m-1} \mathbf{u},\]

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 [88] that combines processes of hard-bedded sliding (Coulomb behavior) and viscous bed deformation without required knowledge of the bed type. It is defined by

(18)¶\[\boldsymbol{\tau}_b = - \tau_c \frac{\uu}{\left(|\uu| + u_{\text{threshold}} \right)^{q} |\uu|^{1-q}},\]

Set -regularized_coulomb to select this sliding law.

The original equation (3) in [88] uses the exponent \(q=1/p\). Otherwise, same configuration parameters can be used as in the pseudo-plastic case, but they should have different values, namely \(q\)=0.2, \(u_{\text{threshold}}\)=40-80 m/year`:

-regularized_coulomb \
-pseudo_plastic_q 0.2 \
-pseudo_plastic_uthreshold 50.0

The model’s performance should be close to the pseudo-plastic 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 spatially-variable 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, or
-yield_stress constant allows the yield stress to be supplied as time-independent 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 [21]. 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 till-pore 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 Mohr-Coulomb criterion [34]:

(19)¶\[\tau_c = c_{0} + (\tan\phi)\,N_{till}.\]

Here \(c_0\) is called the “till cohesion”, whose default value in PISM is zero (see [21], 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.

Table 14 Command-line options controlling how yield stress is determined¶
Option	Description
`-yield_stress mohr_coulomb`	The default. Use equation (19) to determine \(\tau_c\). Only effective if `-stress_balance ssa` or `-stress_balance ssa+sia` is also set.
`-till_cohesion`	Set the value of the till cohesion (\(c_{0}\)) in the plastic till model. The value is a pressure, given in Pa.
`-tauc_slippery_grounding_lines`	If set, reduces the basal yield stress at grounded-below-sea-level grid points one cell away from floating ice or ocean. Specifically, it replaces the normally-computed \(\tau_c\) from the Mohr-Coulomb relation, which uses the effective pressure from the modeled amount of water in the till, by the minimum value of \(\tau_c\) from Mohr-Coulomb, i.e. using the effective pressure corresponding to the maximum amount of till-stored water. Does not alter the reported amount of till water, nor does this mechanism affect water conservation.
`-plastic_phi` (degrees)	Use a constant till friction angle. The default is \(30^{\circ}\).
`-topg_to_phi` (list of 4 numbers)	Compute \(\phi\) using equation (20).
`-yield_stress tillphi_opt`	Compute the till friction angle \(\phi\) in (19) iteratively in an equilibrium simulation using equation (21).
`-yield_stress constant`	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 `-stress_balance ssa` or `-stress_balance ssa+sia` is also set.
`-tauc`	Directly set the till yield stress \(\tau_c\), in units Pa, at all grounded locations and all times. Only effective if used with `-yield_stress constant`, because otherwise \(\tau_c\) is updated dynamically.

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 [2], [52], [17]. This heuristic is motivated by hypothesis that basal material with a marine history should be weak [45]. PISM has a mechanism setting \(\phi\) to be a piece-wise 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

(20)¶\[\begin{split}\phi(x,y) = \begin{cases} \phimin, & b(x,y) \le \bmin, \\ \phimin + (b(x,y) - \bmin) \,M, & \bmin < b(x,y) < \bmax, \\ \phimax, & \bmax \le b(x,y). \end{cases}\end{split}\]

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 option -topg_to_phi is set then this will be set to “yes”, and then MohrCoulombYieldStress will initialize the tillphi field using a piece-wise linear function of depth described by four parameters.
phi_max (15 degrees) upper value of the till friction angle; see the implementation of MohrCoulombYieldStress
phi_min (5 degrees) lower value of the till friction angle; see the implementation of MohrCoulombYieldStress
topg_max (1000 meters) the elevation at which the upper value of the till friction angle is used; see the implementation of MohrCoulombYieldStress
topg_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]\):

(21)¶\[\phi_{n+1} = \min(\max(\phimin, \phi_n + \Delta \phi), \phimax),\]

The value of \(\Delta \phi\) is proportional to the difference between the target (usually observed present day, e.g. [89]) and modeled surface elevations. It is similarly clipped to ensure \(\Delta \phi \in [\dphimin, \dphimax]\):

(22)¶\[ \begin{align}\begin{aligned}\Delta \phi &= \min (\max(\dphimin, \Delta\tilde\phi), \dphimax),\\\Delta \tilde \phi &= C \left( h_{\mathrm{observed}} - h_{\mathrm{modeled}} \right),\\\dphimin &= -2 \dphimax.\end{aligned}\end{align} \]

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 piece-wise linear function of the bed topography \(b\):

(23)¶\[\begin{split}\phi_0(b) = \begin{cases} \phi_{0,\mathrm{min}}, &b \le \bmin,\\ \phi_{0,\mathrm{min}} + (\phi_{0,\mathrm{max}} - \phi_{0,\mathrm{min}}) \frac{b - \bmin}{\bmax - \bmin}, & \bmin < b \le \bmax, \\ \phi_{0,\mathrm{max}} & \bmax < b. \end{cases}\end{split}\]

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 [90]).

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:

(24)¶\[\frac{\Delta h(x,y,T+\dt_{\phi}) - \Delta h(x,y,T) }{\dt_{\phi}} \leq D\]

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 year-1) rate of change in the surface elevation mismatch \(D\) used as a convergence criterion
dphi_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 mismatch
dt (250 365days) time step \(\dt_{\phi}\) for optimization of till friction angle
file Name of the file containing the time-independent variable usurf used as target surface elevation
phi0_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 which basal_yield_stress.mohr_coulomb.tillphi_opt.phi0_max is used
topg_min (-300 meters) the bed elevation \(\bmin\) below which basal_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 command-line option selecting a surface model (see Preventing grounding line retreat), and
fix 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 input-with-mask.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 fields
diff_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 conservation-of-energy (see Modeling conservation of energy) submodels in PISM.

Following [19], based on laboratory experiments with till extracted from an ice stream in Antarctica, [91] 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):

(25)¶\[N_{till} = \min\left\{P_o, N_0 \left(\frac{\delta P_o}{N_0}\right)^s \, 10^{(e_0/C_c) \left(1 - s\right).}\right\}\]

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 time-dependent variable mohr_coulomb_delta to use instead of basal_yield_stress.mohr_coulomb.till_effective_fraction_overburden.
delta.periodic (no) If true, interpret forcing data as periodic in time
tauc_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 [19]
till_effective_fraction_overburden (0.02) \(\delta\) in [91]; \(\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 runs
till_log_factor_transportable_water (0.1 meters) If basal_yield_stress.add_transportable_water is set then the water amount in the transport system is added to tillwat in determining tauc. 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 angle
till_reference_effective_pressure (1000 Pascal) reference effective pressure \(N_0\); value from [19]
till_reference_void_ratio (0.69) void ratio at reference effective pressure \(N_0\); value from [19]

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 [91])

(26)¶\[\delta P_o \le N_{till} \le P_o.\]

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 time-dependent 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 piece-wise linear interpolation in time for model times between records of \(\delta\).

Footnotes

1: This is similar to the simple inversion method described in [92] which optimizes the basal sliding coefficient \(\tau_c\) instead.