# Stress Field Equations for Pattyn 2003 Model

## Conservation Equations

To derive this model, we begin by stating the laws of conservation of mass, momentum, and energy as they apply to an incompressible fluid (Versteed). Consider an infinitely small "control volume" of ice. The law of conservation of mass in this context states that the rate of increase of the mass of the control volume must equal the net rate of flow into the control volume (Versteed). With an incompressible fluid such as ice, the mass of a control volume cannot change, so this law reduces to the fact that any flow into the control volume must be balanced by flow out of the control volume (Versteed). Stated formally:

$\nabla \cdot \mathbf{v} = 0$

Where $\mathbf{v}$ is velocity.

The conservation of momentum is a statement of Newton's second law as applied to fluid dynamics: the rate of change of momentum equals the sum the forces (Versteed). If we consider a control volume of ice, the rate of change of momentum (or the density times the acceleration) equals the sum of the forces due to stress on the control volume and the force of gravity (Pattyn). Formally,

$\rho_i \frac{d\mathbf{v}}{dt} = \nabla \cdot \mathbf{T} + \rho_i \mathbf{g}$

Where $\mathbf{T}$ is the total stress tensor, and $\rho_i$ and $g$ are the density of ice and gravitational acceleration respectively.

The conservation of energy is a statement of the first law of thermodynamics: the rate of change of energy is equal to the rate of heat addition plus the rate of work done (Versteed). Formally,

$\rho_i \frac{d(c_p \theta)}{dt} = \nabla(k_i \nabla \theta) + \Phi$

We now begin to develop the diagnostic model by applying our simplifying assumptions. We neglect the conservation of energy entirely. Although Pattyn's original model included temperature as a concern, in CISM temperature advection is handled by a different module and therefore does not concern this derivation (Glimmer doc). We expand the conservation of mass, using the common notation of referring to velocity components as u, v, and w (Pattyn):

$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0$

For the momentum equation, we neglect the acceleration term due the Froude number of ice being of the order 10-12 (the acceleration term is 1012 times smaller than the gravitational terms)

$0 = \nabla \cdot \sigma + \rho_i \mathbf{g}$

$-\rho_i \mathbf{g} = \nabla \cdot \sigma$

We define the coordinate system so that the top of the ice sheet is at z=0, with positive numbers denoting lower elevations (this introduces a rescaling term that will be discussed). This restricts gravitational influence to the z dimension, and allows us to expand the momentum equation (Pattyn):

$\frac{\partial T_{xx}}{\partial x} + \frac{\partial T_{xy}}{\partial y} + \frac{\partial T_{xz}}{\partial z} = 0$

$\frac{\partial T_{yx}}{\partial x} + \frac{\partial T_{yy}}{\partial y} + \frac{\partial T_{yz}}{\partial z} = 0$

$\frac{\partial T_{zx}}{\partial x} + \frac{\partial T_{zy}}{\partial y} + \frac{\partial T_{zz}}{\partial z} = \rho_i g$

Because vertical resistive stresses are neglected, we drop the $T_{zx}$ and $T_{zy}$ derivatives (Pattyn):

$\frac{\partial T_{zz}}{\partial z} = \rho_i g$

## Deviatoric Stresses

Some additional transformations are needed before directly relating the stresses and velocities. First, it has been established that the strain rates in ice do not largely depend on cryostatic pressure (Hooke). We can estimate the cryostatic pressure as the mean of the normal stresses:

$\frac{1}{3}(T_{xx} + T_{yy} + T_{zz})$

Thus, we define the deviatoric stress tensor as a stress tensor that neglects cryostatic pressure:

\begin{align} T'_{ij} & = T_{ij} - \frac{1}{3}\delta_{ij}(T_{xx} + T_{yy} + T_{zz}) \end{align}

Here, the Kroneker Delta, $\delta_{ij}$, is one if $i=j$, and zero otherwise.

In order to substitute deviatoric stresses for total stresses, we form the following system of equations:

$T_{xx} = T'_{xx} + \frac{1}{3}(T_{xx} + T_{yy} + T_{zz})$

$T_{yy} = T'_{yy} + \frac{1}{3}(T_{xx} + T_{yy} + T_{zz})$

Solving for $T_{xx}$ and $T_{yy}$:

$T_{xx} = 2T'_{xx} + T'_{yy} + T_{zz}$

$T_{yy} = 2T'_{yy} + T'_{xx} + T_{zz}$

We can remove the $\sigma_{zz}$ term by vertically integrating $\frac{\partial \sigma_{zz}}{\partial z} = \rho_i g$ from the surface $s$ to a height $z$ in the ice sheet (Pattyn):

$T_{zz}(z) = -\rho_i g (s - z)$

Thus:

$T_{xx} = 2T'_{xx} + T'_{yy} - \rho_i g (s - z)$

$T_{yy} = 2T'_{yy} + T'_{xx} - \rho_i g (s - z)$

## Stress, Strain, and Viscosity

We relate the deviatoric stress tensor to velocity gradients by equating each with the strain rates. Stress and strain rate are related nonlinearly by a Glen-type flow law (Paterson 1994):

$T_{ij} = 2 \mu \dot \epsilon_{ij}$

with the effective viscosity $\mu$ given by:

$\mu = \frac{A^{\frac{-1}{n}}}{2} \dot \epsilon^\frac{1-n}{n}$

where n is the flow law exponent, and defines the strength of the nonlinearity between stress and strain, usually taken to be 3. A is a thermomechanical coupling parameter, usually given by an Arrhenius relationship (Pattyn). It is not elaborated in this model because its computation is the responsibility of other CISM modules (Glimmer doc); for many basic experiments it is taken as a constant 10-16 (ISMIP-HOM). $\dot \epsilon$ is the second invariant of the strain rate tensor (Pattyn 2003) (Hooke):

$\dot \epsilon = \sqrt{ \dot \epsilon_{xy}^2 + \dot \epsilon_{yz}^2 + \dot \epsilon_{zx}^2 - \dot \epsilon_{xx}\dot \epsilon_{yy} - \dot \epsilon_{yy}\dot \epsilon_{zz} - \dot \epsilon_{zz}\dot \epsilon_{xx}}$

In practice, a small regularization is added to $\dot \epsilon$ in order to avoid a division by zero, particularly in the case of a frozen bed (Pattyn). This formulation of the second invariant will be returned to later.

We now formulate the viscosity term $\mu$ in terms of velocities rather than strain rates. In a full Stokes model, the relationship between strain rates and velocity are defined as follows (Hooke):

$\begin{pmatrix} \dot \epsilon_{xx} & \dot \epsilon_{xy} & \dot \epsilon_{xz} \\ \dot \epsilon_{yx} & \dot \epsilon_{yy} & \dot \epsilon_{yz} \\ \dot \epsilon_{zx} & \dot \epsilon_{zy} & \dot \epsilon_{zz} \\ \end{pmatrix} = \begin{pmatrix} \frac{\partial u}{\partial x} & \frac{1}{2}(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}) & \frac{1}{2}(\frac{\partial u}{\partial z} + \frac{\partial w}{\partial x}) \\ \frac{1}{2}(\frac{\partial v}{\partial x} + \frac{\partial u}{\partial y}) & \frac{\partial v}{\partial y} & \frac{1}{2}(\frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}) \\ \frac{1}{2}(\frac{\partial w}{\partial x} + \frac{\partial u}{\partial z}) & \frac{1}{2}(\frac{\partial w}{\partial y} + \frac{\partial v}{\partial z}) & \frac{\partial w}{\partial z} \\ \end{pmatrix}$

Because of the assumption made above that we can neglect the horizontal gradients of the vertical velocity as they are much smaller than the vertical gradients of the horizontal velocity (Pattyn), this becomes:

$\begin{pmatrix} \dot \epsilon_{xx} & \dot \epsilon_{xy} & \dot \epsilon_{xz} \\ \dot \epsilon_{yx} & \dot \epsilon_{yy} & \dot \epsilon_{yz} \\ \dot \epsilon_{zx} & \dot \epsilon_{zy} & \dot \epsilon_{zz} \\ \end{pmatrix} = \begin{pmatrix} \frac{\partial u}{\partial x} & \frac{1}{2}(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}) & \frac{1}{2}\frac{\partial u}{\partial z} \\ \frac{1}{2}(\frac{\partial v}{\partial x} + \frac{\partial u}{\partial y}) & \frac{\partial v}{\partial y} & \frac{1}{2}\frac{\partial v}{\partial z} \\ \frac{1}{2}\frac{\partial u}{\partial z} & \frac{1}{2}\frac{\partial v}{\partial z} & \frac{\partial w}{\partial z} \\ \end{pmatrix}$

From a modeling standpoint, this simplifying assumption has helped us because it allows us to neglect the vertical component of velocity when solving for the velocity fields, then reconstruct it in another module using the incompressibility condition. In order to remove our dependence on the vertical velocity entirely, we use the law of conservation of mass to remove all $\dot \epsilon_{zz}$ factors from the second invarant of the strain rate tensor (Pattyn). Rearranging terms from the conservation of mass equation:

$\frac{\partial w}{\partial z} = -(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y})$

From the strain rate tensor definition above, this is equivalent to:

$\dot \epsilon_{zz} = -(\dot \epsilon_{xx} + \dot \epsilon_{yy})$

Performing this substitution in the second invariant of the strain rate tensor given above and rearranging terms leads to an alternate form that has no dependence on the vertical velocity:

$\dot \epsilon = \sqrt{ \dot \epsilon_{xy}^2 + \dot \epsilon_{yz}^2 + \dot \epsilon_{xz}^2 + \dot \epsilon_{xx}\dot \epsilon_{yy} + \dot \epsilon_{yy}^2 + \dot \epsilon_{xx}^2}$

We can now finally expand the strain rate invariant factor in the viscosity. and perform substitutions such that viscosity is formulated in terms of velocity gradients rather than strain rates (Pattyn):

$\mu = \frac{A^{\frac{-1}{n}}}{2} \dot \epsilon^\frac{1-n}{n}$

$\mu = \frac{A^{\frac{-1}{n}}}{2} (\dot \epsilon_{xy}^2 + \dot \epsilon_{yz}^2 + \dot \epsilon_{xz}^2 + \dot \epsilon_{xx}\dot \epsilon_{yy} + \dot \epsilon_{yy}^2 + \dot \epsilon_{xx}^2)^\frac{1-n}{2n}$

$\mu = \frac{A^{\frac{-1}{n}}}{2} \left( (\frac{1}{4} (\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x})^2 + \frac{1}{4} (\frac{\partial u}{\partial z})^2 + \frac{1}{4} (\frac{\partial v}{\partial z})^2 + (\frac{\partial u}{\partial x})^2 + (\frac{\partial v}{\partial y})^2 + \frac{\partial u}{\partial z}\frac{\partial v}{\partial y} \right)$

## Stress and Velocity

With the relationship of both deviatoric stresses and velocities to strain rates, we can now reformulate the conservation of momentum purely in terms of velocity. This will lead us to a set of coupled elliptic partial differential equations that can be numerically approximated to find the u and v components of the velocity field. To derive the u component, let us return to the first equation in the system resulting from the conservation of momentum:

$\frac{\partial T_{xx}}{\partial x} + \frac{\partial T_{xy}}{\partial y} + \frac{\partial T_{xz}}{\partial z} = 0$

We want to express this equation in terms of the deviatoric stress rather than the total stress. Using the substitutions derived above (Pattyn):

$\frac{\partial}{\partial x}(2T'_{xx} + T'_{yy}) + \frac{\partial T'_{xy}}{\partial y} + \frac{\partial T'_{xz}}{\partial z} = \rho_i g \frac{\partial s}{\partial x}$

$\frac{\partial}{\partial x}(2T'_{xx} + T'_{yy}) + \frac{\partial T'_{xy}}{\partial y} + \frac{\partial T'_{yz}}{\partial z} = \rho_i g \frac{\partial s}{\partial x}$

Using Glen's flow law and the definition of strain rates given above, the stress tensor can be written as follows:

$\begin{pmatrix} T'_{xx} & T'_{xy} & T'_{xz} \\ T'_{yx} & T'_{yy} & T'_{yz} \\ T'_{zx} & T'_{zy} & T'_{zz} \\ \end{pmatrix} = \begin{pmatrix} 2 \mu \dot \epsilon_{xx} & 2 \mu \dot \epsilon_{xy} & 2 \mu \dot \epsilon_{xz} \\ 2 \mu \dot \epsilon_{yx} & 2 \mu \dot \epsilon_{yy} & 2 \mu \dot \epsilon_{yz} \\ 2 \mu \dot \epsilon_{zx} & 2 \mu \dot \epsilon_{zy} & 2 \mu \dot \epsilon_{zz} \\ \end{pmatrix} = \begin{pmatrix} 2 \mu \frac{\partial u}{\partial x} & \mu (\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}) & \mu \frac{\partial u}{\partial z} \\ \mu (\frac{\partial v}{\partial x} + \frac{\partial u}{\partial y}) & 2 \mu \frac{\partial v}{\partial y} & \mu \frac{\partial v}{\partial z} \\ \mu \frac{\partial u}{\partial z} & \mu \frac{\partial v}{\partial z} & 2 \mu \frac{\partial w}{\partial z} \\ \end{pmatrix}$

Directly substituting into the force balance equations gives us (Pattyn):

$\frac{\partial}{\partial x}(4 \mu \frac{\partial u}{\partial x} + 2 \mu \frac{\partial v}{\partial y}) + \frac{\partial}{\partial y}(\mu \frac{\partial u}{\partial y} + \mu \frac{\partial v}{\partial x}) + \frac{\partial}{\partial z}(\mu \frac{\partial u}{\partial z}) = \rho_i g \frac{\partial s}{\partial x}$

$\frac{\partial}{\partial y}(4 \mu \frac{\partial v}{\partial y} + 2 \mu \frac{\partial u}{\partial x}) + \frac{\partial}{\partial x}(\mu \frac{\partial u}{\partial y} + \mu \frac{\partial v}{\partial x}) + \frac{\partial}{\partial z}(\mu \frac{\partial v}{\partial z}) = \rho_i g \frac{\partial s}{\partial x}$

Finally, we expand the nested derivatives and rearrange terms so that each equation has the terms related to the variable we are solving for (u and v respectively) on the left-hand side, and terms related to the orthogonal velocity component (v and u respectively) on the right-hand side:

$4 \frac{\partial \mu}{\partial x} \frac{\partial u}{\partial x} + \frac{\partial \mu}{\partial y} \frac{\partial u}{\partial y} + \frac{\partial \mu}{\partial z} \frac{\partial u}{\partial z} + \mu(4 \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}) = \rho_i g \frac{\partial s }{\partial x} - 2 \frac{\partial \mu}{\partial x} \frac{\partial v}{\partial y} - \frac{\partial \mu}{\partial y} \frac{\partial v}{\partial x} - 3 \mu \frac{\partial^2 v}{\partial x \partial y}$

$4 \frac{\partial \mu}{\partial y} \frac{\partial v}{\partial y} + \frac{\partial \mu}{\partial x} \frac{\partial v}{\partial x} + \frac{\partial \mu}{\partial z} \frac{\partial v}{\partial z} + \mu(4 \frac{\partial^2 v}{\partial y^2} + \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial z^2}) = \rho_i g \frac{\partial s }{\partial x} - 2 \frac{\partial \mu}{\partial y} \frac{\partial u}{\partial x} - \frac{\partial \mu}{\partial x} \frac{\partial u}{\partial y} - 3 \mu \frac{\partial^2 u}{\partial x \partial y}$

## Rescaled Vertical Coordinate

Pattyn introduces a dimensionless vertical coordinate so that variations in ice thickness do not complicate the numerical treatment considerably (Pattyn 2003). This dimensionless vertical coordinate is referred to as $\zeta$ by Pattyn. I will refer to it as $\sigma$ for consistency with CISM's existing notation. The rescaled coordinate is defined as:

$\sigma = \frac{(s - z)}{H}$

This means that at the surface of the ice sheet $\sigma = 0$, and at the base $\sigma = 1$ regardless of the ice thickness (Pattyn 2003). As a result of this transformation, a coordinate $(x,y,z)$ is mapped to $(x',y',\sigma)$ (Pattyn 2003). This means that function derivatives must be re-written (using $\frac{\partial f}{\partial x}$ as an example) as:

$\frac{\partial f}{\partial x} = \frac{\partial f}{\partial x'} \frac{\partial x'}{\partial x} + \frac{\partial f}{\partial y'} \frac{\partial y'}{\partial x} + \frac{\partial f}{\partial \sigma} \frac{\partial \sigma}{\partial x}$

Similarly for $\frac{\partial f}{\partial y}$ and $\frac{\partial f}{\partial z}$. Pattyn simplifies this by assuming that

$\frac{\partial x'}{\partial x}, \frac{\partial y'}{\partial y} = 1$

and

$\frac{\partial x'}{\partial y}, \frac{\partial x'}{\partial z}, \frac{\partial y'}{\partial x}, \frac{\partial y'}{\partial z} = 0$.

This assumption is valid if the bed and surface gradients are not too large (Pattyn 2003). This simplifies the above to:

$\frac{\partial f}{\partial x} = \frac{\partial f}{\partial x'} + \frac{\partial f}{\partial \sigma}\frac{\partial \sigma}{\partial x}$

$\frac{\partial f}{\partial y} = \frac{\partial f}{\partial y'} + \frac{\partial f}{\partial \sigma}\frac{\partial \sigma}{\partial y}$

$\frac{\partial f}{\partial z} = \frac{\partial f}{\partial \sigma}\frac{\partial \sigma}{\partial z}$

Rescaling parameters $a_x$, $a_y$, $b_x$, $b_y$, and $c_{xy}$ are defined. Presenting the x derivative case, as the y derivative case is analogous,

$a_x = \frac{1}{H}(\frac{\partial s}{\partial x'} - \sigma \frac{\partial H}{\partial x'})$

$b_x = \frac{\partial a_x}{\partial x'} + a_x \frac{\partial a_x}{\partial \sigma} = \frac{1}{H} (\frac{\partial^2 s}{\partial x'^2} - \sigma \frac{\partial^2 H}{\partial x'^2} - 2a_x \frac{\partial H}{\partial x'})$

$c_{xy} = \frac{\partial a_y}{\partial x'} + a_x \frac{\partial a_y}{\partial \sigma} = \frac{\partial a_x}{\partial y'} + a_y \frac{\partial a_x}{\partial \sigma}$

Using these, expressions for the x derivatives become:

$\frac{\partial f}{\partial x} = \frac{\partial f}{\partial x'} + a_x \frac{\partial f}{\partial \sigma}$