Shallow ice approximation (SIA)

Note

The explicit time stepping of the mass continuity equation in the case of the SIA flow comes with a severe restriction on time step length:

(4)\[\dt \le \frac{2 R}{D\left( 1/\dx^2 + 1/\dy^2 \right)}\]

Here \(D\) is the maximum diffusivity of the SIA flow and \(R\) is time_stepping­.adaptive_ratio, a tuning parameter that further reduces the maximum allowed time step length.

The maximum diffusivity \(D\) may be achieved at an isolated grid point near the ice margin. In this case it might make sense to limit the diffusivity of the SIA flow, sacrificing accuracy at a few grid points to increase time step length and reduce the computational cost. Set stress_balance­.sia­.limit_diffusivity to enable this mechanism.

When stress_balance­.sia­.limit_diffusivity is false PISM stops as soon as the SIA diffusivity at any grid point exceeds stress_balance­.sia­.max_diffusivity. We do this to make it easier to detect problematic model configurations: in many cases it does not make sense to continue a simulation if \(D\) is very large.

Surface gradient method

PISM computes surface gradients to determine the “driving stress”

\[(\tau_{d,x},\tau_{d,y}) = - \rho g H \nabla h,\]

where \(H\) is the ice thickness, and \(h\) is the ice surface elevation.

In the SIA model surface gradients at staggered grid locations are computed using one of the following three finite-difference approximations (selected using stress_balance­.sia­.surface_gradient_method):

  1. mahaffy: This most “standard” way computes the surface slope onto the staggered grid for the SIA [64]. It makes \(O(\dx^2,\dy^2)\) errors.

  2. haseloff: This is the default method. It only differs from mahaffy at ice-margin locations, where the slope is approximated using one-sided finite differences in cases where an adjacent ice-free bedrock surface elevation is above the ice elevation.

  3. eta: In this method we first transform the thickness \(H\) by \(\eta = H^{(2n+2)/n}\) and then differentiate the sum of the thickness and the bed using centered differences:

    \[\nabla h = \nabla H + \nabla b = \frac{n}{(2n+2)} \eta^{(-n-2)/(2n+2)} \nabla \eta + \nabla b.\]

    Here \(b\) is the bed elevation, \(h\) is the surface elevation, and \(n\) is the Glen exponent. This transformation sometimes has the benefits that the surface values of the horizontal velocity and vertical velocity, and the driving stress, are better behaved near the margin. See [58] for technical explanation of this transformation and compare [65]. The actual finite difference schemes applied to compute the surface slope are similar to option mahaffy.

    Note

    This method may improve the model performance near grounded margins but should not be used in simulations of marine ice sheets.

To the best of our knowledge there is no theoretical advice on the best, most robust mechanism.

Parameterization of bed roughness

Schoof [66] describes how to alter the SIA stress balance to model ice flow over bumpy bedrock topgraphy. One computes the amount by which bumpy topography lowers the SIA diffusivity. An internal quantity used in this method is a smoothed version of the bedrock topography. As a practical matter for PISM, this theory improves the SIA’s ability to handle bed roughness because it parameterizes the effects of “higher-order” stresses which act on the ice as it flows over bumps. For additional technical description of PISM’s implementation, see Using Schoof’s parameterized bed roughness technique in PISM.

There are two associated parameters:

This mechanism is turned off by default in pismv.

Under the default output­.size (medium), PISM writes fields topgsmooth and schoofs_theta from this mechanism. The thickness relative to the smoothed bedrock elevation, namely topgsmooth, is the difference between the unsmoothed surface elevation and the smoothed bedrock elevation. It is only used internally by this mechanism, to compute a modified value of the diffusivity; the rest of PISM does not use this or any other smoothed bed. The field schoofs_theta is a number \(\theta\) between stress_balance­.sia­.bed_smoother­.theta_min (usually zero) and \(1\), with values significantly below one indicating a reduction in diffusivity, essentially a drag coefficient, from bumpy bed.

Coupling to the age of the ice

The age of the ice can be used in two parameterizations in the SIA stress balance model:

  1. Ice grain size parameterization based on data from [67] and [68] (Vostok core data). In PISM, only the Goldsby-Kohlstedt flow law (see Ice rheology) uses the grain size.

    Set stress_balance­.sia­.grain_size_age_coupling to enable this parameterization.

  2. The flow enhancement factor can be coupled to the age of the ice as in [69]: \(e\) in (12) is set to

    See time­.eemian_start, time­.eemian_end, and time­.holocene_start.

    Set stress_balance­.sia­.e_age_coupling to enable this parameterization.

Parameters

Prefix: stress_balance.sia.

  1. Glen_exponent (3) Glen exponent in ice flow law for SIA

  2. bed_smoother­.range (5000 meters) half-width of smoothing domain in the bed roughness parameterization for SIA [66]; set to zero to disable

  3. bed_smoother­.theta_min (0) minimum value of \(\theta\) in the bed roughness parameterization for SIA [66]

  4. e_age_coupling (no) Couple the SIA enhancement factor to age as in [44].

  5. enhancement_factor (1) Flow enhancement factor for SIA

  6. enhancement_factor_interglacial (1) Flow enhancement factor for SIA; used for ice accumulated during interglacial periods.

  7. flow_law (gpbld) The SIA flow law.

  8. grain_size_age_coupling (no) Use age of the ice to compute grain size to use with the Goldsby-Kohlstedt [70] flow law

  9. limit_diffusivity (no) Limit SIA diffusivity by stress_balance­.sia­.max_diffusivity.

  10. max_diffusivity (100 m2 s-1) Maximum allowed diffusivity of the SIA flow. PISM stops with an error message if the SIA diffusivity exceeds this limit.

  11. surface_gradient_method (haseloff) method used for surface gradient calculation at staggered grid points

Previous Up Next