SIA — Shallow Ice Approximation

The SIA is the complementary limit to SSA: membrane stresses are neglected and the horizontal stress balance is dominated by the basal drag exerted by vertical shear. Yelmo’s SIA solver is implemented in src/physics/velocity_sia.f90 and is used to provide the vertical-shear contribution to the velocity field in the hybrid SIA + SSA mode and as a diagnostic on its own.

Crucially, Yelmo formulates SIA in stress form: the shear stress and the basal stress are computed from the driving stress \boldsymbol\tau_d exactly as in DIVA, and the velocity is then obtained by integrating Glen’s flow law vertically. This is mathematically equivalent to the familiar closed-form SIA expression in terms of \nabla s, but it keeps the link to the Stokes balance explicit, uses the same \boldsymbol\tau_d field as DIVA/SSA (including the same lateral and grounding-line treatments), and makes the SIA → DIVA limit transparent.

Continuum equations

The horizontal momentum balance, neglecting membrane stresses and vertical accelerations, reduces to a balance between the vertical gradient of horizontal shear stress and the horizontal pressure gradient of the gravitational driving:

\frac{\partial \tau_{xz}}{\partial z} \;=\; \rho_i\, g\, \frac{\partial s}{\partial x}, \qquad \frac{\partial \tau_{yz}}{\partial z} \;=\; \rho_i\, g\, \frac{\partial s}{\partial y}.

Integrating from z upward to the stress-free surface and writing the result in terms of the same driving stress as DIVA/SSA gives the linear-with-depth profile

\tau_{xz}(z) \;=\; -\,\frac{s-z}{H}\,\tau_{d,x}, \qquad \tau_{yz}(z) \;=\; -\,\frac{s-z}{H}\,\tau_{d,y},

with the driving stress

\tau_{d,x} \;=\; \rho_i\, g\, H\, \frac{\partial s}{\partial x}, \qquad \tau_{d,y} \;=\; \rho_i\, g\, H\, \frac{\partial s}{\partial y},

evaluated on ac-faces by the same calc_driving_stress routine used by DIVA and SSA (see DIVA). The basal stress is recovered at z = b: (\tau_{xz}, \tau_{yz})\big|_{z=b} = -\boldsymbol\tau_d, i.e. all of the driving stress is balanced at the bed, as expected when membrane stresses are dropped. The shear stress profile is built by calc_shear_stress_3D.

Velocity from Glen’s flow law

With the shear-stress field in hand, the horizontal shear strain rates follow directly from Glen’s law,

\dot\varepsilon_{xz} \;=\; A(T'(z))\, \tau_e^{\,n-1}\,\tau_{xz}, \qquad \dot\varepsilon_{yz} \;=\; A(T'(z))\, \tau_e^{\,n-1}\,\tau_{yz},

where the effective stress under SIA assumptions is \tau_e^{\,2} = \tau_{xz}^{\,2} + \tau_{yz}^{\,2}. Since \dot\varepsilon_{xz} = \frac{1}{2}\frac{\partial u}{\partial z} in the SIA limit, the horizontal velocity is obtained by integrating from the bed upward:

u(z) \;=\; u_b \,+\, 2\int_b^z A(T'(z'))\,\tau_e^{\,n-1}\,\tau_{xz}\,\mathrm dz',

and analogously for v. The integral is evaluated layer-by-layer with the trapezoidal rule (routine calc_uxy_sia_3D). The depth-averaged velocity \bar{\mathbf u}^{(\mathrm{SIA})} is obtained by trapezoidal integration over z.

In hybrid mode, Yelmo uses \bar{\mathbf u}^{(\mathrm{SIA})} to provide the shearing contribution that is added to the basal velocity produced by the membrane (DIVA/SSA) solve. In a pure-SIA configuration, the basal velocity is either zero (frozen bed) or comes from a sliding law of choice (not the default in current Yelmo versions).

Relationship to the textbook SIA

Substituting the closed-form \tau_{xz} = -\rho_i\,g\,(s-z)\,\frac{\partial s}{\partial x} into the integral above and assuming A constant with z recovers the familiar SIA velocity

u(z) \;=\; u_b \,-\, \frac{2\,(\rho_i g)^n\,A}{n+1}\, H^{n+1}\,|\nabla s|^{\,n-1}\,\frac{\partial s}{\partial x}\, \bigl[\,1 - (1 - (z-b)/H)^{n+1}\,\bigr],

but Yelmo does not collapse to that form internally. Keeping the stress-form integral is what allows the SIA solver to share \boldsymbol\tau_d, A(T'(z)) and the discretisation conventions with the DIVA solver and the rest of the model.