Boundary Conditions for Pattyn 2003 Model

From Interactive System for Ice sheet Simulation
Jump to: navigation, search


Stress-Free Surface

The stress-free boundary condition occurs at the surface of the ice sheet and specifies that the shear stress parallel to the surface is zero (Ven Der Veen 1999). This condition can be formalized (MacAyael Lessons) as:

\mathbf{T} \cdot \mathbf{\hat n_s} = 0

Where \mathbf{T} is the total (not deviatoric) stress tensor and \mathbf{\hat n_s} is a unit vector normal to the surface of the ice sheet that we derive by normalizing:

\mathbf{n_s} = \begin{pmatrix}-\frac{\partial s}{\partial x} \\ -\frac{\partial s}{\partial y} \\ 1 \end{pmatrix}

Expanding the tensor dot product and ignoring the vertical resistive stresses gives us the force balance for the upper boundary condition:

T_{xx} \frac{\partial s}{\partial x} + T_{xy} \frac{\partial s}{\partial y} - T_{xz} = 0

T_{yx} \frac{\partial s}{\partial x} + T_{yy} \frac{\partial s}{\partial y} - T_{yz} = 0

-T_{zz} = 0

Using the equations for moving from total stress to deviatoric stress, and observing that the term \rho_i g (s-z) is zero since s=z, we finally arrive at the stress-free surface in terms of deviatoric stresses (Pattyn 2003):

(2T'_{xx} + T'_{yy})\frac{\partial s}{\partial x} + T'_{xy} \frac{\partial s}{\partial y} - T'_{xz} = 0

(2T'_{xx} + T'_{yy})\frac{\partial s}{\partial y} + T'_{yx} \frac{\partial s}{\partial x} - T'_{yz} = 0

Applying the constitutive relationship (and noting that constant factors, including viscosity, can be dropped because of the equation with zero) (Pattyn 2003):

(4\frac{\partial u}{\partial x} + 2\frac{\partial v}{\partial y})\frac{\partial s}{\partial x} + (\frac{\partial u}{\partial u} + \frac{\partial v}{\partial x})\frac{\partial s}{\partial y} - \frac{\partial u}{\partial z} = 0

(4\frac{\partial v}{\partial y} + 2\frac{\partial v}{\partial y})\frac{\partial s}{\partial y} + (\frac{\partial u}{\partial u} + \frac{\partial v}{\partial x})\frac{\partial s}{\partial x} - \frac{\partial v}{\partial z} = 0

Finally, applying the vertical coordinate rescaling (Pattyn 2003):

4 \frac{\partial s}{\partial x'} \frac{\partial u}{\partial x'} 
+ \frac{\partial s}{\partial y'} \frac{\partial u}{\partial y'} 
+ (4 a_x \frac{\partial s}{\partial x'} + a_y \frac{\partial s}{\partial y'} + \frac{1}{H})\frac{\partial u}{\partial \sigma} 
= -\frac{\partial s}{\partial y'} \frac{\partial v}{\partial x'} 
- 2 \frac {\partial s}{\partial x'} \frac{\partial v}{\partial y'} 
- (2 a_y \frac{\partial s}{\partial x'} + a_x \frac{\partial s}{\partial y'}) \frac{\partial v}{\partial \sigma}

4 \frac{\partial s}{\partial y'} \frac{\partial v}{\partial y'} 
+ \frac{\partial s}{\partial x'} \frac{\partial v}{\partial x'} 
+ (4 a_y \frac{\partial s}{\partial y'} + a_x \frac{\partial s}{\partial x'} + \frac{1}{H})\frac{\partial v}{\partial \sigma} 
= -\frac{\partial s}{\partial x'} \frac{\partial u}{\partial y'} 
- 2 \frac {\partial s}{\partial y'} \frac{\partial u}{\partial x'} 
- (2 a_x \frac{\partial s}{\partial y'} + a_y \frac{\partial s}{\partial x'}) \frac{\partial u}{\partial \sigma}

Basal Boundary Condition

Significantly different outcomes can result from the model for bed strength that is used. The two most commonly used models are for plastic and power law materials. Plastic bed models are believed to be appropriate for ice that rests on marine sediments, whereas a power law relation for bed strength may be more appropriate for ice resting on a layer of exposed bedrock.

Hard Bed

A model for basal traction in cases where the bed consists of exposed bedrock is formulated as a power law material.

\tau_{bx} = \beta^2 u^{\left(\frac{1}{m} - 1\right)} |\mathbf{v}|

\tau_{by} = \beta^2 v^{\left(\frac{1}{m} - 1\right)} |\mathbf{v}|

At the moment, Pattyn's model allows only a linear bed, which is a special case with m = \frac{1}{2}.

To derive the force balance in this case, we divorce ourselves from a particular coordinate system and observe that given any unit vector tangential to the bedrock \mathbf{\hat t}, the traction can be expressed by taking the projection of basal drag onto that vector (Pattyn 2008):

\beta^2 \mathbf{\hat t} \cdot \mathbf{v} = \mathbf{\hat t} \cdot (\mathbf{T} \mathbf{\hat n_b}) = \tau_b

where \mathbf{\hat n_b} is a unit vector normal to the lower surface and pointing into the bedrock that we obtain by normalizing</math>

n_b = \begin{pmatrix}\frac{\partial b}{\partial x} \\ \frac{\partial b}{\partial y} \\ -1\end{pmatrix}

Observe that, because of the dot product, this produces only one equation. Since we are solving for the x and y components of velocity, we need two. But, because the equation is true for any tangent vector, we can produce a linearly independant set of equations by choosing linearly independant tangent vectors. Choosing vectors that will break \tau_b into components \tau_{bx} and \tau_{by}:

\mathbf{t_x} = \begin{pmatrix}1 \\ 0 \\ \frac{\partial b}{\partial x}\end{pmatrix}

\mathbf{t_y} = \begin{pmatrix}0 \\ 1 \\ \frac{\partial b}{\partial y}\end{pmatrix}

This gives us the system of equations:

\beta^2 \frac{\mathbf{t_x}}{||\mathbf{t_x}||} \cdot \mathbf{v} = \frac{\mathbf{t_x}}{||\mathbf{t_x}||} \cdot (\mathbf{T} \frac{\mathbf{n_b}}{||\mathbf{n_b}||}) = \tau_bx

\beta^2 \frac{\mathbf{t_y}}{||\mathbf{t_y}||} \cdot \mathbf{v} = \frac{\mathbf{t_y}}{||\mathbf{t_y}||} \cdot (\mathbf{T} \frac{\mathbf{n_b}}{||\mathbf{n_b}||}) = \tau_by

Expanding the vector and tensor products:

\beta^2 ||\mathbf{n_b}|| (u + w\frac{\partial b}{\partial x}) = T_{xx}\frac{\partial b}{\partial x} + T_{xy}\frac{\partial b}{\partial y} - T_{xz} - T_{zz}\frac{\partial b}{\partial x}

\beta^2 ||\mathbf{n_b}|| (u + w\frac{\partial b}{\partial y}) = T_{yx}\frac{\partial b}{\partial x} + T_{yy}\frac{\partial b}{\partial y} - T_{yz} - T_{zz}\frac{\partial b}{\partial y}

We make the assumption here that the vertical velocity component at the base is much larger than the horizontal component, which is valid as long as bedrock gradients are not too large (??) and melt rates are negligable. Hence, we drop the w term in the right-hand-side and substitute deviatoric stresses for total stresses, observing that the pressure at the base is \rho_i g H:

\beta^2 u ||\mathbf{n_b}|| = (2T'_{xx} + T'_{yy} - \rho_i g H)\frac{\partial b}{\partial x} + T'_{xy}\frac{\partial b}{\partial y} - T'_{xz} + \rho_i g H \frac{\partial b}{\partial x}

\beta^2 v ||\mathbf{n_b}|| = (2T'_{yy} + T'_{xx} - \rho_i g H)\frac{\partial b}{\partial y} + T'_{yx}\frac{\partial b}{\partial x} - T'_{yz} + \rho_i g H \frac{\partial b}{\partial y}

The pressure terms cancel and we apply the constitutive relationship:

\beta^2 u ||\mathbf{n_b}|| = 2 \mu (2 \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y})\frac{\partial b}{\partial x} + \mu (\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x})\frac{\partial b}{\partial y} - \mu \frac{\partial u}{\partial z}

\beta^2 v ||\mathbf{n_b}|| = 2 \mu (2 \frac{\partial v}{\partial y} + \frac{\partial u}{\partial x})\frac{\partial b}{\partial y} + \mu (\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x})\frac{\partial b}{\partial x} - \mu \frac{\partial u}{\partial z}

Finally, with rescaled vertical coordinates:

-\beta^2 u ||\mathbf{n_b}||
+ 4\mu \frac{\partial b}{\partial x'}\frac{\partial u}{\partial x'}
+ \mu \frac{\partial b}{\partial y'}\frac{\partial u}{\partial y'}
+ \mu (4 \frac{\partial b}{\partial x'} a_x + \frac{\partial b}{\partial y'}a_y - \frac{1}{H})\frac{\partial u}{\partial \sigma}
- 2\mu \frac{\partial b}{\partial x'} \frac{\partial v}{\partial y'} 
- \mu \frac{\partial b}{\partial y'} \frac{\partial v}{\partial x'}
- \mu (2\frac{\partial b}{\partial x'} a_y + \frac{\partial b}{\partial y'} a_x)\frac{\partial v}{\partial \sigma}

-\beta^2 v ||\mathbf{n_b}||
+ 4\mu \frac{\partial b}{\partial y'}\frac{\partial v}{\partial y'}
+ \mu \frac{\partial b}{\partial x'}\frac{\partial v}{\partial x'}
+ \mu (4 \frac{\partial b}{\partial y'} a_y + \frac{\partial b}{\partial x'}a_x - \frac{1}{H})\frac{\partial v}{\partial \sigma}
- 2\mu \frac{\partial b}{\partial y'} \frac{\partial u}{\partial x'} 
- \mu \frac{\partial b}{\partial x'} \frac{\partial u}{\partial y'}
- \mu (2\frac{\partial b}{\partial y'} a_x + \frac{\partial b}{\partial x'} a_y)\frac{\partial u}{\partial \sigma}

Soft Bed

The formulation for basal traction due to plastic bed is

\tau_{bx} = \tau_0 \frac{u}{ |\mathbf{v}|}

\tau_{by} = \tau_0 \frac{v}{|\mathbf{v}|}.

Where |\mathbf{v}| = (u^2 + v^2)^{\frac{1}{2}}

The rest of the formulation in the linear bed case is valid for this case as well, though with \beta^2 = \frac{\tau_0}{(u^2 + v^2)^{\frac{1}{2}}}. This is a much harder problem to solve because the strength of the nonlinearity is greater.

Frozen Bed

In addition to the dynamic boundary conditions, Pattyn provides the ability to specify that ice is frozen to the bed, as is the case in several of the ISMIP-HOM experiments (ISMIP-HOM). In this case, the finite differencing is set up so that velocity at the base is held at zero.

Lower Ice Shelf Boundary

Finally, we must consider the lower boundary condition of an ice shelf. In this case, we balance the outward-facing lower surface normal with hydrostatic pressure (MacAyeal Lessons). Assuming the ice shelf is in hydrostatic equilibrium, the lower surface will be at elevation -\frac{\rho_i}{\rho_w}H. Then,

\mathbf{T} \mathbf{\hat n_b} = -\rho_w g (\frac{\rho_i}{\rho_w}H)\mathbf{\hat n_b}

Where again, \hat n_b is the outward-facing normal unit vector to the ice surface that is a normalization of:

\mathbf{\hat n_b} = \begin{pmatrix}\frac{\partial b}{\partial x} \\ \frac{\partial b}{\partial y} \\ -1 \end{pmatrix}

Unlike the linear bed case above, we do not formulate this in terms of tangents to the lower surface because resistance to the velocity is not involved; the formulation of this boundary condition most closely resembles the stress-free boundary condition above. Expanding the tensor product and canceling the normalization factor gives us three equations:

T_{xx} \frac{\partial b}{\partial x} + T_{xy} \frac{\partial b}{\partial y} - T_{xz} = -\rho_i g H \frac{\partial b}{\partial x}

T_{xy} \frac{\partial b}{\partial x} + T_{yy} \frac{\partial b}{\partial y} - T_{yz} = -\rho_i g H \frac{\partial b}{\partial y}

- T_{zz} = \rho_i g H

Notice that we have simplified the right-hand side in a somewhat revealing way. Because we have assumed hydrostatic equilibrium, of course the ice shelf will float so as to equalize pressure inside and outside of it. This is exactly what the right-hand-side, which is the same as the cryostatic pressure at the base of the ice sheet, is telling us.

Next we substitute in deviatoric stresses. When we do this, the final equation becomes tautological (again reflecting hydrostatic equilibrium) so we remove it and are left with two equations:

(2T'_{xx} + T'_yy + \rho_i g H) \frac{\partial b}{\partial x} + T'_{xy} \frac{\partial b}{\partial y} - T'_{xz} = \rho_i g H \frac{\partial b}{\partial x}

T'_{xy} \frac{\partial b}{\partial x} + (2T'_{yy} + T'_{xx} + \rho_i g H) \frac{\partial b}{\partial y} - T'_{yz} = \rho_i g H \frac{\partial b}{\partial y}

Notice here that we can cancel \frac{\partial b}{\partial x}\rho_i g H and \frac{\partial b}{\partial y}\rho_i g H, leaving us with a stress-free base. From here on the basal boundary condition can be easily adapted from the stress-free surface by replacing surface gradients with bed gradients (Pattyn 2003).

Lateral Boundary Conditions

Two types of lateral boundary conditions have been introduced by us as extensions to Pattyn's original model to support ice shelf intercomparison experiments (e.g. Ross).

Kinematic Boundary

We have added the ability to specify a kinematic (Dirichlet) boundary condition at any point in the domain; this consists simply of specifying the velocity directly at a particular point. In particular, is used to specify ice stream inputs to the Ross ice shelf (MacAyeal et al 1996). Similarly to the frozen bed described above, implementation of this boundary condition consists simply of specifying in the finite difference scheme that a point is held at that value.

Dynamic Boundary

The second type of boundary condition handled is the dynamic (Neumann) boundary condition that occurs at the marine margin of an ice shelf. This is often formulated as a vertically integrated value assuming that the ice shelf is in hydrostatic equilibrium (MacAyeal et al 1996):

\int_{-\frac{\rho_i}{\rho_w}H}^{(1 - \frac{\rho_i}{\rho_w})H} \mathbf{T} \cdot \mathbf{n} dz 
-\frac{\rho_w g}{2}(\frac{\rho_i}{\rho_w}H)^2 \mathbf{n}

This integral assumes a coordinate system where z is strictly elevation (that is, z=0 at sea level, values increase as we ascend).

The boundary is often formulated this way because ice shelf models do not need to be vertically explicit, as ice shelves exhibit "plug flow". Hence, traditionally ice shelf models have used the "shallow shelf approximation" and are vertically integrated. However, as Pattyn's model is vertically explicit, we make the boundary vertically explicit by differentiating both sides of the equation with respect to z. Because the right-hand side is the integral of the hydrostatic pressure on the ice sheet, to "differentiate" it with respect to z we are really replacing it with the hydrostatic pressure equation. Below the water level, this becomes:

\mathbf{T}(z) \cdot \mathbf{n} = \rho_w g z \mathbf{n}, z < 0

Above, due to the approximation of atmospheric pressure as 0, it becomes:

\mathbf{T}(z) \cdot \mathbf{n} = 0, z \ge 0

For most of the following derivation, I will neglect the latter case and focus on the case where z < 0. Applying the tensor product and using the fact that n_z = 0, this leads to the system of equations:

T_{xx} n_x + T_{xy} n_y = \rho_w g z n_x

T_{yx} n_x + T_{yy} n_y = \rho_w g z n_y

T_{zx} n_x + T_{zy} n_y = 0

Using the definitions of deviatoric stresses:

(2T'_{xx} + T'_{yy} - \rho_i g (z_s - z)) n_x + T'_{xy} n_y = \rho_w g z n_x

T'_{yx} n_x + (2T'_{yy} + T'_{xx} - \rho_i g (z_s - z)) n_y = \rho_w g z n_y

T'_{zx} n_x + T'_{zy} n_y = 0

Applying the constitutive relationship:

  2 \mu(2 \dot \epsilon_{xx} + \dot \epsilon_{yy}) n_x + 2 \mu \dot \epsilon_{xy} n_y - \rho_i g (z_s-z) n_x = \rho_w g z n_x

  2 \mu(\dot \epsilon_{xx} + 2 \dot \epsilon_{yy}) n_y + 2 \mu \dot \epsilon_{yx} n_x - \rho_i g (z_s-z) n_y = \rho_w g z n_y

 2 \mu (\dot \epsilon_{xz} n_x + \dot \epsilon_{yz} n_y)= 0

Replacing strain rates with velocities:

 2 \mu [(2 \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}) n_x + \frac{1}{2}(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}) n_y]= n_x(\rho_i g (z_s-z) + \rho_w g z)

 2 \mu [\frac{1}{2}(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}) n_x + (\frac{\partial u}{\partial x} + 2\frac{\partial v}{\partial y}) n_y]= n_y(\rho_i g (z_s-z) + \rho_w g z)

 \mu (\frac{\partial u}{\partial z} n_x + \frac{\partial v}{\partial z} n_y) = 0

Above the surface, the water pressure term is 0 reducing the lateral boundary to a stress-free boundary condition:

 2 \mu [(2 \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}) n_x + \frac{1}{2}(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}) n_y]= n_x \rho_i g (z_s - z)

 2 \mu [\frac{1}{2}(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}) n_x + (\frac{\partial u}{\partial x} + 2\frac{\partial v}{\partial y}) n_y]= n_y \rho_i g (z_s - z)

 \mu (\frac{\partial u}{\partial z} n_x + \frac{\partial v}{\partial z} n_y) = 0

Additionally, here we assume plug flow on the ice shelf, meaning that the vertical derivatives of horizontal velocities are both zero. This removes the third equation from the system.

We must now recast the equations to use Pattyn's rescaled vertical coordinate. To do this, we must perform derivative substitutions so that variables are written in terms of (x', y', \sigma). We must also change the right-hand side as it has an explicit vertical dependence. In the sigma coordinate system, z can be substituted with:

z = z_s - H\sigma

Where z is the original vertical coordinate (z=0 at sea level, increasing upwards), H is the ice thickness, s is the ice surface, and \sigma is the rescaled vertical coordinate. Assuming hydrostatic equilibrium,

s = H(1 - \frac{\rho_i}{\rho_w})


z = H(1 - \frac{\rho_i}{\rho_w} - \sigma)

Above the surface, the pressure term (again assuming 0 atmospheric pressure) becomes:

\rho_i g (z_s - z) = \rho_i g H \sigma

Below the surface, the pressure term becomes:

\rho_i g (z_s-z) - \rho_w g z = \rho_i g H \sigma + \rho_w g H (1 - \frac{\rho_i}{\rho_w} - \sigma)

Thus, we re-write the above equations as:

 2 \mu [(2 \frac{\partial u}{\partial x'} + \frac{\partial v}{\partial y'} + 2 a_x\frac{\partial u}{\partial \sigma} + a_y \frac{\partial v}{\partial \sigma}) n_x + \mu (\frac{\partial u}{\partial y'} + \frac{\partial v}{\partial x'} + a_y\frac{\partial u}{\partial \sigma} + a_x\frac{\partial v}{\partial \sigma}) n_y]= n_x(\rho_i g H \sigma + \rho_w g H (1 - \frac{\rho_i}{\rho_w} - \sigma))

 2 \mu [(\frac{\partial u}{\partial x'} + 2\frac{\partial v}{\partial y'} + a_x \frac{\partial u}{\partial \sigma} + 2 a_y\frac{\partial v}{\partial \sigma} ) n_y + \mu (\frac{\partial u}{\partial y'} + \frac{\partial v}{\partial x'} + a_y\frac{\partial u}{\partial \sigma} + a_x\frac{\partial v}{\partial \sigma}) n_x]= n_y(\rho_i g H \sigma + \rho_w g H (1 - \frac{\rho_i}{\rho_w} - \sigma))

Rearranging terms to aid differencing:

4\mu n_x \frac{\partial u}{\partial x'} 
+ \mu n_y \frac{\partial u}{\partial y'} 
+ \mu(4 a_x n_x + a_y n_y)\frac{\partial u}{\partial \sigma}
= n_x(\rho_i g H \sigma + \rho_w g H (1 - \frac{\rho_i}{\rho_w} - \sigma)) 
- 2 \mu n_x \frac{\partial v}{\partial y'} 
- \mu n_y \frac{\partial v}{\partial x'}
- \mu(2 a_y n_x + a_x n_y) \frac{\partial v}{\partial \sigma}

4\mu n_y \frac{\partial v}{\partial y'} 
+ \mu n_x \frac{\partial v}{\partial x'} 
+ \mu(4 a_y n_y + a_x n_x)\frac{\partial v}{\partial \sigma}
= n_y(\rho_i g H \sigma + \rho_w g H (1 - \frac{\rho_i}{\rho_w} - \sigma)) 
- 2 \mu n_y \frac{\partial u}{\partial x'} 
- \mu n_x \frac{\partial u}{\partial y'}
- \mu (2 a_x n_y + a_y n_x) \frac{\partial u}{\partial \sigma}