Solution of the Blatter-Pattyn model
The final form of the equations we'd like to solve is:
For ice sheet modeling, it is convenient to recast the governing equations using a dimensionless, stretched vertical coordinate (often called a sigma coordinates). The stretched vertical coordinate is defined as:
This means that at the surface of the ice sheet , and at the base regardless of the ice thickness. As a result of this transformation, a coordinate is mapped to . This means that function derivatives must be re-written (using as an example) as:
Similarly for and . We can simplify this by assuming that
This assumption is valid if the bed and surface gradients are not too large. This simplifies the above to:
Rescaling parameters , , , , and are defined. For the x derivative case (the y derivative case is analogous) we have
Using these, expressions for the x derivatives become:
where hatted values refer to the coordinate directions in sigma coordinates. Similarly, the first cross-stress term on the RHS is given by
One term has become five terms and each one of those is pretty ugly looking on its own. Luckily, there is a lot of symmetry here. Notice that if we wanted to design subroutines to discretize the terms on the RHS, we could re-use a lot of them by either applying them to the correct velocity component (to either the U or the V discretization) or by passing the appropriate arguments (by passing either the grid spacing in the X direction or the Y direction, where appropriate).
A similar transform is applied to each of the terms in the governing equations given above. At any point within the grid, the grid spacing, coordinate transform, and viscosity information associated with the unknown velocity components (U and V) is made discrete using finite differences. This information ultimately equates to coefficients on the unknown velocities, allowing the governing equations over the entire grid (with appropriate discretizations for boundary conditions) to be recast as a system of n equations in n unknowns. In turn, this system is solved using standard linear algebraic methods for large, [sparse] systems of equations.
In the governing equations given above, note that for the x equation we have moved all terms containing gradients in v to the right-hand side (RHS) (and vice-versa for the y equation).
This allows us to solve the equations using an operator splitting approach; for the x equation, we treat v as known (where we take the values of v from the previous iteration, as discussed further below) and solve for u, and vice versa when we solve the y equation for v. The "splitting" refers to the fact that we are breaking the multi-dimensional divergence operation into two steps; rather than solving one big matrix equation for u and v simultaneously, we solve two smaller matrix equations in sequence with one of the unknowns treated as a known "source" term. This procedure was probably more common and important years ago when it was desirable to keep the matrix equations as small as possible for memory management issues. On today's machines, with fewer memory limitations (in particular when dealing with codes designed to run on parallel, distributed memory architectures) this splitting is not necessary and may even lead to some undesirable numerical side effects (i.e. a slow-down in the convergence of iterations used to treat nonlinearity in the governing equations).
General Matrix Form
A general matrix form of the equations, where coefficients on the u and v velocity components (i.e. viscosity, grid spacing, scalars) are contained in the block matrices A, is given by
where the uu subscript denotes block matrices containing coefficients for gradients on u in the equation for the x component of velocity (i.e. u). The subscript uv denotes block matrices containing coefficients for gradients on v in the equation for the x component of velocity (and similarly for the vv and vu subscripts). On the right-hand side, the subscripts denote the geometric source terms for the x and y components of velocity (subscripts u and v, respectively).
Solution of the Non-linear System Through a Fixed Point Iteration
The non-linearity in the equations - the fact that the coefficients on the velocity components (the viscosity) are dependent on the velocity (or more specifically, the velocity gradients) - is handled through a fixed-point iteration.
Go to Strang for this ..
Final Matrix Form
The final matrix form of the equations, accounting for the Picard iteration on the viscosity, is given by
where the index k denotes an unknown value being solved for during the current non-linear iteration and the index k-1 denotes a lagged value taken from solution at the end of the previous non-linear iteration.
Solution of the Linear System
Newton-based Methods for Solutions of the Non-linear System
Go to Blatter-Pattyn model.