Shallow shelf approximation (SSA)¶

If the SSA stress balance is used, a choice of two solvers is available, namely -ssa_method fd (default) or -ssa_method fem. See Table 9, which describes additional controls on the numerical solution of the stress balance equations. If option -ssa_method fd is chosen then several more controls on numerics are available; see Table 10. If the ice sheet being modeled has any floating ice then the user is advised to read section PIK options for marine ice sheets on modeling marine ice sheets.

When using SSA as a “sliding law” one also needs to model the yield stress, or a pseudo-yield-stress in the case of power law sliding (section Controlling basal strength).

The basal yield stress is normally a function of the amount of water stored in the till and a (generally) spatially-varying till strength. The amount of stored basal water is modeled by the subglacial hydrology model (section Subglacial hydrology) based on the basal melt rate which is, primarily, thermodynamically-determined (see Modeling conservation of energy).

Table 9 Choice of, and controls on, the numerical SSA stress balance.¶
Option	Description
`-ssa_method` [ `fd \| fem` ]	Both finite difference (`fd`; the default) and finite element (`fem`) versions of the SSA numerical solver are implemented in PISM. The `fd` solver is the only one which allows PIK options (section PIK options for marine ice sheets). `fd` uses Picard iteration [10], while `fem` uses a Newton method. The `fem` solver has surface velocity inversion capability [27].
`-ssa_eps` (\(10^{13}\))	The numerical schemes for the SSA compute an effective viscosity \(\nu\) which depends on strain rates and ice hardness (thus temperature). The minimum value of the effective viscosity times the thickness (i.e. \(\nu H\)) largely determines the difficulty of solving the numerical SSA. This constant is added to keep \(\nu H\) bounded away from zero: \(\nu H \to \nu H + \epsilon_{\text{SSA}}\), where \(\epsilon_{\text{SSA}}\) is set using this option. Units of `ssa_eps` are \(\text{Pa}\,\text{m}\,\text{s}\). Set to zero to turn off this lower bound.
`-ssa_view_nuh`	View the product \(\nu H\) for your simulation as a runtime viewer (section Run-time diagnostic viewers). In a typical Greenland run we see a wide range of values for \(\nu H\) from \(\sim 10^{14}\) to \(\sim 10^{20}\) \(\text{Pa}\,\text{m}\,\text{s}\).

Table 10 Controls on the numerical iteration of the `-ssa_method fd` solver¶
Option	Description
`-ssafd_picard_maxi` (300)	Set the maximum allowed number of Picard (nonlinear) iterations in solving the shallow shelf approximation.
`-ssafd_picard_rtol` (\(10^{-4}\))	The Picard iteration computes a vertically-averaged effective viscosity which is used to solve the equations for horizontal velocity. Then the new velocities are used to recompute an effective viscosity, and so on. This option sets the relative change tolerance for the effective viscosity. The Picard iteration stops when successive values \(\nu^{(k)}\) of the vertically-averaged effective viscosity satisfy \[\\|(\nu^{(k)} - \nu^{(k-1)}) H\\|_1 \le Z \\|\nu^{(k)} H\\|_1\] where \(Z=\) `ssafd_picard_rtol`.
`-ssafd_ksp_rtol` (\(10^{-5}\))	Set the relative change tolerance for the iteration inside the Krylov linear solver used at each Picard iteration.
`-ssafd_max_speed` (\(50 km/yr\))	Limits computed SSA velocities: ice speed is capped at this limit after each Picard iteration of the SSAFD solver. This may allow PISM to take longer time steps by ignoring high velocities at a few troublesome locations.

Parameters¶

Prefix: stress_balance.ssa.

Glen_exponent (3) Glen exponent in ice flow law for SSA
compute_surface_gradient_inward (no) If yes then use inward first-order differencing in computing surface gradient in the SSA objects.
dirichlet_bc (no) apply SSA velocity Dirichlet boundary condition
enhancement_factor (1) Flow enhancement factor for SSA
epsilon (1e+13 Pascal second meter) Initial amount of regularization in computation of product of effective viscosity and thickness (\(\nu H\)). This default value for \(\nu H\) comes e.g. from a hardness for the Ross ice shelf (\(\bar B\)) = 1.9e8 Pa \(s^{1/3}\) [62] and a typical strain rate of 0.001 1/year for the Ross ice shelf, giving \(\nu = (\bar B) / (2 \cdot 0.001^{2/3})\) = 9.49e+14 Pa s ~ 30 MPa year, the value in [63], but with a tiny thickness \(H\) of about 1 cm.
fd.brutal_sliding (false) Enhance sliding speed brutally.
fd.brutal_sliding_scale (1) Brutal SSA Sliding Scale
fd.extrapolate_at_margins (true) Extrapolate SSA velocity to ice-free cells near the margin to improve SSA velocities at advancing ice fronts
fd.flow_line_mode (false) Set \(v\) (the \(y\) component of the ice velocity) to zero when assembling the system
fd.lateral_drag.enabled (false) Set viscosity at ice shelf margin next to ice free bedrock as friction parameterization
fd.lateral_drag.viscosity (5e+15 Pascal second) Staggered viscosity used as side friction parameterization.
fd.max_iterations (300) Maximum number of Picard iterations for the ice viscosity computation, in the SSAFD object
fd.max_speed (300000 km s-1) Upper bound for the ice speed computed by the SSAFD solver.
fd.nuH_iter_failure_underrelaxation (0.8) In event of “Effective viscosity not converged” failure, use outer iteration rule nuH <- nuH + f (nuH - nuH_old), where f is this parameter.
fd.relative_convergence (0.0001) Relative change tolerance for the effective viscosity in the SSAFD object
fd.replace_zero_diagonal_entries (yes) Replace zero diagonal entries in the SSAFD matrix with :config:’basal_resistance.beta_ice_free_bedrock’ to avoid solver failures.
flow_law (gpbld) The SSA flow law.
method (fd) Algorithm for computing the SSA solution.
read_initial_guess (yes) Read the initial guess from the input file when re-starting.
strength_extension.constant_nu (9.48681e+14 Pascal second) The SSA is made elliptic by use of a constant value for the product of viscosity (nu) and thickness (H). This value for nu comes from hardness (bar B)=1.9e8 \(Pa s^{1/3}\) [62] and a typical strain rate of 0.001 year-1: \(\nu = (\bar B) / (2 \cdot 0.001^{2/3})\). Compare the value of 9.45e14 Pa s = 30 MPa year in [63].
strength_extension.min_thickness (50 meters) The SSA is made elliptic by use of a constant value for the product of viscosity (nu) and thickness (H). At ice thicknesses below this value the product nu*H switches from the normal vertical integral to a constant value. The geometry itself is not affected by this value.

Technical remarks¶

Ice stress balance models ignoring inertia are “diagnostic” models that do not have a “state” evolving through time: the ice velocity is fully determined by ice geometry, basal boundary conditions, and the ice viscosity.

In addition to this, shallow stress balance models (other than the SIA) correspond to nonlinear systems of equations, which means that computing an estimate of the ice velocity is an iterative process starting with a particular initial guess.

The quality of this guess may affect the number of iterations needed and even whether the solver succeeds at all. It also has an influence on the computed solution: if an equation has a unique solution, estimates produced using different initial guesses should be close, but they need not be identical. If an equation has multiple solutions, even a “small” change in the initial guess may give a completely different result.

In the context of a evolutionary model we can usually assume that the change in the state of the model from one step to the next is “small” and we use (\(u, v\)) estimates from one time step as the initial guess for the next. To ensure that stopping and re-starting a simulation does not affect results we save these to the output file as variables u_ssa and v_ssa and read them in when re-starting a stopped simulation.

One could say that the continuum SSA model does not have a state, but its implementation does. Set stress_balance.ssa.read_initial_guess to “false” to ignore it during initialization and use the zero initial guess instead.