tSSANeumBoundCond.tex - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tSSANeumBoundCond.tex (23380B)
---
1 \documentclass[a4paper,10pt]{article}
2 \usepackage{graphicx}
3 \usepackage{amsmath}
4 \usepackage{dsfont}
5 \usepackage{comment}
6 \usepackage{subfigure}
7 \usepackage[english]{babel}
8 \usepackage[top=3cm, bottom=3cm, left=3cm, right=3cm]{geometry}
9 \title{Implementation of Boundary Conditions at the Calving Front}
10 \author{Marianne Haseloff \& Torsten Albrecht,\\ Maria Martin, Ricarda Winkelmann, Anders Levermann}
11 \begin{document}
12 \maketitle
13 \section{Theoretics}
14 The force balance at calving fronts has been formulated by Morland~\cite{Morland87} in the following way:
15 \begin{equation}%MacAyeal1
16 \int_{-\frac{\rho}{\rho_w}H}^{(1-\frac{\rho}{\rho_w})H}\mathbf{\sigma}\cdot\mathbf{n}dz = -\frac{\rho_w}{2}g\left(\frac{\rho}{\rho_w}H \right)^2\mathbf{n}
17 \label{MacAyeal1}
18 \end{equation}
19 with $\mathbf{n}$ being the normal vector pointing from the ice
20 oceanward, $\mathbf{\sigma}$ the \emph{Cauchy} stress tensor, $H$ the ice thickness and $\rho$ and $\rho_{w}$ the densities of ice and seawater, respectively. A slightly different form allows for changing sealevel $z_s$:
21 \begin{equation}
22 \int_{z_s-\frac{\rho}{\rho_w}H}^{z_s+(1-\frac{\rho}{\rho_w})H}\mathbf{\sigma}\cdot\mathbf{n}dz = \int_{z_s-\frac{\rho}{\rho_w}H}^{z_s}\rho_w g (z-z_s) dz\mathbf{n}.
23 \label{MacAyeal2}
24 \end{equation}
25 The integration limits on the right hand side of Eq.~\eqref{MacAyeal2} account for the pressure exerted by the ocean on that part of the shelf, which is below sea level. For grounded ice the following calculations are similar but with different integration limits:
26 \begin{equation}
27 \int_{b}^{H+b}\mathbf{\sigma}\cdot\mathbf{n}dz = \int_{b}^{z_s}\rho_w g (z-z_s) dz\mathbf{n}.
28 \label{BC_sheet}
29 \end{equation}
30 Integration of the right-hand side of \eqref{MacAyeal2} yields:
31 \begin{eqnarray}
32 \int_{z_s-\frac{\rho}{\rho_w}H}^{z_s+(1-\frac{\rho}{\rho_w})H}\mathbf{\sigma}\cdot\mathbf{n}dz & = & \rho_w g \left(\frac{z^2}{2}-z_s z\right)\mid_{z_s-\frac{\rho}{\rho_w}H}^{z_s}\mathbf{n} \\
33 & = & -\frac{\rho_w}{2}g\left(\frac{\rho}{\rho_w}H \right)^2\mathbf{n} \label{MacAyeal3}
34 \end{eqnarray}
35 \noindent Using cartesian coordinates,
36 i.e. $\mathbf{n}=n_x\overrightarrow{e_x}+n_y\overrightarrow{e_y}$, equation
37 \eqref{MacAyeal3} can be rewritten in terms of tensor components:
38 \begin{eqnarray*}
39 \int_{h_m}^{h_s}(n_x\sigma_{xx}+n_y\sigma_{xy})dz & = & -\frac{\rho_w}{2}g\left(\frac{\rho}{\rho_w}H\right)^2n_x \\
40 \int_{h_m}^{h_s}(n_x\sigma_{yx}+n_y\sigma_{yy})dz & = & -\frac{\rho_w}{2}g\left(\frac{\rho}{\rho_w}H\right)^2n_y \\
41 \int_{h_m}^{h_s}(n_x\sigma_{zx}+n_y\sigma_{zy})dz & = & 0
42 \end{eqnarray*}
43 $h_s=z_s+(1-\frac{\rho}{\rho_w})H$ and
44 $h_m=z_s-\frac{\rho}{\rho_w}H$ denote the ice shelf's upper and lower boundary,
45 respectively. As PISM deals with velocities instead of stresses, the above equations have to be rewritten in terms of velocities. The stress is related to the strain rate $\dot{\epsilon_{ij}}$ by:
46 \begin{equation*}
47 \sigma_{ij}=\frac{1}{\left[A(T)\right]^{\frac{1}{n}}}d^{\frac{1-n}{n}} \dot{\epsilon_{ij}} - p\delta_{ij}
48 \end{equation*}
49 and
50 \begin{equation*}
51 \dot{\epsilon_{ij}}=\frac{1}{2}\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i} \right)
52 \end{equation*}
53 Hence, the velocities enter into the above equations via:
54 \begin{eqnarray*}
55 \int_{h_m}^{h_s}\sigma_{xx}dz & = & \int_{h_m}^{h_s}\left(\frac{1}{\left[A(T)\right]^{\frac{1}{n}}}d^{\frac{1-n}{n}}\frac{\partial u}{\partial x}-p \right) dz\\
56 \int_{h_m}^{h_s}\sigma_{yy}dz & = & \int_{h_m}^{h_s}\left(\frac{1}{\left[A(T)\right]^{\frac{1}{n}}}d^{\frac{1-n}{n}}\frac{\partial v}{\partial y}-p \right) dz\\
57 \int_{h_m}^{h_s}\sigma_{xy}dz & = & \int_{h_m}^{h_s}\frac{1}{2}\frac{1}{\left[A(T)\right]^{\frac{1}{n}}}d^{\frac{1-n}{n}}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right)dz \\
58 \end{eqnarray*}
59 It is an essential result of the SSA, that horiontal velocities do not depend on depth. Thus, the above equations can be further simplified:
60 \begin{eqnarray*}
61 \int_{h_m}^{h_s}\sigma_{xx}dz & = & \frac{\partial u}{\partial x} \int_{h_m}^{h_s}\frac{1}{\left[A(T)\right]^{\frac{1}{n}}}d^{\frac{1-n}{n}}dz-\int_{h_m}^{h_s}p dz \\
62 \int_{h_m}^{h_s}\sigma_{yy}dz & = & \frac{\partial v}{\partial y} \int_{h_m}^{h_s}\frac{1}{\left[A(T)\right]^{\frac{1}{n}}}d^{\frac{1-n}{n}}dz-\int_{h_m}^{h_s}p dz \\
63 \int_{h_m}^{h_s}\sigma_{xy}dz & = & \frac{1}{2}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right)\int_{h_m}^{h_s}\frac{1}{\left[A(T)\right]^{\frac{1}{n}}}d^{\frac{1-n}{n}}dz
64 \end{eqnarray*}
65 Defining the vertically integrated effective viscosity as
66 \begin{equation}%myNu: effective viscosity
67 \tilde{\nu}=\frac{1}{\rho g}d^{\frac{1-n}{n}}\int_{h_m}^{h_s}\frac{1}{\left[A(T)\right]^{\frac{1}{n}}}dz
68 \label{myNu}
69 \end{equation}
70 further simplification yields:
71 \begin{eqnarray}%tesnorX tensorY tensorZ
72 \int_{h_m}^{h_s}\sigma_{xx}dz & = &
73 \rho g\tilde{\nu} \frac{\partial u}{\partial x} - \int_{h_m}^{h_s}p dz \label{tensorX}\\
74 \int_{h_m}^{h_s}\sigma_{yy}dz & = &
75 \rho g\tilde{\nu} \frac{\partial v}{\partial y} - \int_{h_m}^{h_s}p dz \label{tensorY}\\
76 \int_{h_m}^{h_s}\sigma_{xy}dz & = &
77 \frac{1}{2} \rho g\tilde{\nu} \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right) \label{tensorXY}
78 \end{eqnarray}
79 By calculating the depth-integrated pressure, an expression for the derivative of the horizontal velocity at the calving front can be obtained. The pressure can be shown to satisfy (cp.~\cite{Weis01}):
80 \begin{eqnarray*}
81 -p(z) & = & \sigma_z - \tau_z \\
82 & = & \rho g\left[z-\left(z_s-\varrho H\right)\right] -\frac{1}{A(T)^\frac{1}{n}}d^\frac{1-n}{n}\frac{\partial v_z}{\partial z}
83 \end{eqnarray*}
84 with $\varrho=1-\frac{\rho}{\rho_w}$ the reduced density and $\mathbf{\tau}$ the deviatoric stress tensor. Furthermore, the equation of continuity is valid within the SSA limit, too and ice is assumed to be density preserving. Thus:
85 \begin{equation*}
86 \frac{\partial v_i}{\partial x_i} = 0,
87 \end{equation*}
88 and as horiziontal velocities do not depend on depth, $\partial v_z/ \partial
89 z$ may not depend on z, too. This allows to rewrite the left side of eq.\eqref{tensorX}:
90 \begin{eqnarray*}
91 \int_{h_m}^{h_s}\sigma_{xx}dz & = & \rho g\tilde{\nu} \frac{\partial u}{\partial x} + \int_{h_m}^{h_s} \left[ \rho g(z-z_s) - \varrho\rho gH+\frac{1}{A(T)^\frac{1}{n}}d^\frac{1-n}{n} \left( \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} \right) \right] dz \\
92 & = & \rho g\tilde{\nu} \frac{\partial u}{\partial x} + \frac{1}{2}\rho gH^2 - \rho g \frac{\rho}{\rho_w}H^2 -\varrho\rho gH^2+\rho g \tilde{\nu}\left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right) \\
93 & = & \rho g \left(2\tilde{\nu}\frac{\partial u}{\partial x} + \tilde{\nu}\frac{\partial v}{\partial y} - \frac{1}{2}H^2 \right)
94 \end{eqnarray*}
95 and of eq.\eqref{tensorY}, accordingly. Therefore, the boundary conditions at the calving front for an arbitrary normal vector $\mathbf{n}=n_x\overrightarrow{e}_x+n_y\overrightarrow{e}_y$ in terms of velocity gradients are:
96 \begin{eqnarray*}
97 n_x\left(2\tilde{\nu}\frac{\partial u}{\partial x}+\tilde{\nu}\frac{\partial v}{\partial y} \right) + \frac{n_y}{2} \left(\tilde{\nu}\frac{\partial u}{\partial y} + \tilde{\nu}\frac{\partial v}{\partial x} \right) & = & \frac{n_x}{2}\left(1-\frac{\rho}{\rho_w}\right)H^2\\
98 \frac{n_x}{2} \left(\tilde{\nu}\frac{\partial u}{\partial y} + \tilde{\nu}\frac{\partial v}{\partial x} \right) + n_y\left(\tilde{\nu}\frac{\partial u}{\partial x}+2\tilde{\nu}\frac{\partial v}{\partial y} \right) & = & \frac{n_y}{2}\left(1-\frac{\rho}{\rho_w}\right)H^2
99 \end{eqnarray*}
100 The effective viscosity $\tilde{\nu}$ defined by \eqref{myNu} and the viscosity used by PISM $\overline{\nu}$ are linked via:
101 \begin{equation*}%nuLink
102 2H\overline\nu\approx\rho g\tilde{\nu}
103 \label{nuLink}
104 \end{equation*}
105 Thus, the calving front boundary conditions at a \textbf{shelf} front can be formulated as:
106 \begin{eqnarray}%calv_BC1 and calv_BC2
107 n_x\left(2\overline{\nu}H\frac{\partial u}{\partial x}+\overline{\nu}H\frac{\partial v}{\partial y} \right) + \frac{n_y}{2} \left(\overline{\nu}H\frac{\partial u}{\partial y} + \overline{\nu}H\frac{\partial v}{\partial x} \right) & = & n_x\frac{\rho g}{4}\left(1-\frac{\rho}{\rho_w}\right)H^2 \label{calv_BC1}\\
108 n_y\left(\overline{\nu}H\frac{\partial u}{\partial x}+2\overline{\nu}H\frac{\partial v}{\partial y} \right) + \frac{n_x}{2} \left(\overline{\nu}H\frac{\partial u}{\partial y} + \overline{\nu}H\frac{\partial v}{\partial x} \right) & = & n_y\frac{\rho g}{4}\left(1-\frac{\rho}{\rho_w}\right)H^2 . \label{calv_BC2}
109 \end{eqnarray}
110 Similar calculations to rewrite \eqref{BC_sheet} yield the boundary condition on a \textbf{grounded} front:
111 \begin{eqnarray}%sheet_BC1 and sheet_BC2
112 n_x\left(2\overline{\nu}H\frac{\partial u}{\partial x}+\overline{\nu}H\frac{\partial v}{\partial y} \right) + \frac{n_y}{2} \left(\overline{\nu}H\frac{\partial u}{\partial y} + \overline{\nu}H\frac{\partial v}{\partial x} \right) & = & n_x\frac{\rho g}{4}\left(H^2-\frac{\rho_w}{\rho}( z_s-b)^2 \right) \label{sheet_BC1}\\
113 n_y\left(\overline{\nu}H\frac{\partial u}{\partial x}+2\overline{\nu}H\frac{\partial v}{\partial y} \right) + \frac{n_x}{2} \left(\overline{\nu}H\frac{\partial u}{\partial y} + \overline{\nu}H\frac{\partial v}{\partial x} \right) & = & n_y\frac{\rho g}{4}\left(H^2-\frac{\rho_w}{\rho}(z_s - b)^2 \right). \label{sheet_BC2}
114 \end{eqnarray}
115 As the only difference between shelf and sheet is the right hand side of the boundary condition (this is not surprising as the only difference is the area on which the pressure is exerted), both boundary conditions can be formulated as one:
116 \begin{eqnarray}%sheet_BC1 and sheet_BC2
117 n_x\left(2\overline{\nu}H\frac{\partial u}{\partial x}+\overline{\nu}H\frac{\partial v}{\partial y} \right) + \frac{n_y}{2} \left(\overline{\nu}H\frac{\partial u}{\partial y} + \overline{\nu}H\frac{\partial v}{\partial x} \right) & = & n_x\tau_{ocean} \label{BC1}\\
118 n_y\left(\overline{\nu}H\frac{\partial u}{\partial x}+2\overline{\nu}H\frac{\partial v}{\partial y} \right) + \frac{n_x}{2} \left(\overline{\nu}H\frac{\partial u}{\partial y} + \overline{\nu}H\frac{\partial v}{\partial x} \right) & = & n_y\tau_{ocean}. \label{BC2}
119 \end{eqnarray}
120 with
121 \begin{equation}
122 \tau_{ocean} = \left\{
123 \begin{array}{ll}
124 \frac{\rho g}{4}\left(1-\frac{\rho}{\rho_w}\right)H^2 & \text{if shelf} \\
125 \frac{\rho g}{4}\left(H^2-\frac{\rho_w}{\rho} ( z_s-b)^2 \right) & \text{if sheet}
126 \end{array} \right.
127 \end{equation}
128 Obviously, there is no need to normalize the `normal vector' $\mathbf{n}$.
129 %\newpage
130 \section{Discretization}
131
132 \noindent The calving front boundary condition determines the behavior of the velocity gradients at the calving front. To understand the implementation of the calving front boundary condition consider the SSA-equations:
133
134 \begin{align}
135 \frac{\partial}{\partial x}\left[ 2\bar\nu H\left( 2\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right)\right] + \frac{\partial }{\partial y}\left[\bar\nu H\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right) \right] &= \rho gH \frac{\partial h}{\partial x} \label{SSA1} \\
136 \frac{\partial}{\partial x}\left[ \bar\nu H\left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)\right] + \frac{\partial }{\partial y}\left[2\bar\nu H\left(\frac{\partial u}{\partial x}+2\frac{\partial v}{\partial y} \right) \right] &= \rho gH \frac{\partial h}{\partial y} \label{SSA2}
137 \end{align}
138 \noindent In the first step of discretization the outer derivatives of \eqref{SSA1} yield:
139 \begin{align}
140 \frac{2}{\Delta x}(\bar{\nu} H)|_{i+\frac{1}{2}}^j\left( 2\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right)_{i+\frac{1}{2}}^j
141 - \frac{2}{\Delta x}(\bar{\nu} H)|_{i-\frac{1}{2}}^j\left( 2\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right)_{i-\frac{1}{2}}^j \nonumber \\
142 + \frac{1}{\Delta y}(\bar{\nu} H)|_i^{j+\frac{1}{2}}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)_i^{j+\frac{1}{2}}
143 - \frac{1}{\Delta y}(\bar{\nu} H)|_i^{j-\frac{1}{2}}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)_i^{j-\frac{1}{2}}
144 &= \rho gH \frac{\partial h}{\partial x} \label{SSA1_dis1}
145 \end{align}
146 If one or more of the four nearest neighbors of the shelf box under consideration $[i,j]$ is an partially filled grid cell or ice-free ocean (\cite{Albrecht_Martin10}), the boundary between these two boxes will be a calving front and evaluating the boundary condition at these points allows to substitute the whole block which is calculated at this point:
147 \begin{itemize}
148 \item calving front at $[i+\frac{1}{2},j]$: $\hat{n}=(1,0)$:
149 \begin{align}
150 (\bar{\nu} H)|_{i+\frac{1}{2}}^j\left(2\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right)_{i+\frac{1}{2}}^j &= \tau_{ocean} \nonumber \\
151 (\bar{\nu} H)|_{i+\frac{1}{2}}^j\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)_{i+\frac{1}{2}}^j &= 0
152 \end{align}
153 \item calving front at $[i-\frac{1}{2},j]$: $\hat{n}=(-1,0)$:
154 \begin{align}
155 (\bar{\nu} H)|_{i-\frac{1}{2}}^j\left(2\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right)_{i-\frac{1}{2}}^j &= \tau_{ocean} \nonumber \\
156 (\bar{\nu} H)|_{i-\frac{1}{2}}^j\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)_{i-\frac{1}{2}}^j &= 0
157 \end{align}
158 \item calving front at $[i,j+\frac{1}{2}]$: $\hat{n}=(0,1)$:
159 \begin{align}
160 (\bar{\nu} H)|_{i}^{j+\frac{1}{2}}\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)_{i}^{j+\frac{1}{2}} &= 0 \nonumber \\
161 (\bar{\nu} H)|_{i}^{j+\frac{1}{2}}\left(\frac{\partial u}{\partial x} + 2\frac{\partial v}{\partial y} \right)_{i}^{j+\frac{1}{2}} &= \tau_{ocean}
162 \end{align}
163 \item calving front at $[i,j-\frac{1}{2}]$: $\hat{n}=(0,-1)$:
164 \begin{align}
165 (\bar{\nu} H)|_{i}^{j-\frac{1}{2}}\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)_{i}^{j-\frac{1}{2}} &= 0 \nonumber \\
166 (\bar{\nu} H)|_{i}^{j-\frac{1}{2}}\left(\frac{\partial u}{\partial x} + 2\frac{\partial v}{\partial y} \right)_{i}^{j-\frac{1}{2}} &= \tau_{ocean}
167 \end{align}
168 \end{itemize}
169 \begin{figure}[htb]
170 \begin{center}
171 \includegraphics[height=80mm, width=8cm]{f07.pdf}
172 \caption{\emph{SSA stencil}}
173 \label{result3}
174 \end{center}
175 \end{figure}
176 By rewriting equation \eqref{SSA1_dis1} with booleans $a_+, a_-, b_+, b_-$ equal zero if the boundary is a calving front and one otherwise (see Fig.~\ref{result3}), these boundary conditions can be taken into account:
177 \begin{align}
178 &\qquad a_+\frac{2}{\Delta x}(\bar{\nu} H)|_{i+\frac{1}{2}}^j\left( 2\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right)_{i+\frac{1}{2}}^j
179 - a_-\frac{2}{\Delta x}(\bar{\nu} H)|_{i-\frac{1}{2}}^j\left( 2\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right)_{i-\frac{1}{2}}^j \nonumber \\
180 + &\qquad b_+\frac{1}{\Delta y}(\bar{\nu} H)|_i^{j+\frac{1}{2}}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)_i^{j+\frac{1}{2}}
181 - b_-\frac{1}{\Delta y}(\bar{\nu} H)|_i^{j-\frac{1}{2}}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)_i^{j-\frac{1}{2}} \nonumber \\
182 = &\qquad \rho gH \frac{\partial h}{\partial x} + (a_+ - a_-)\frac{2}{\Delta x}\tau_{ocean} \label{SSA1_dis2}
183 \end{align}
184 Next, the calculation of the velocity gradients in \eqref{SSA1_dis2} is further specified by rewriting these in a way that allows to calculate the velocity derivatives as differential quotient between two neighboring points
185 \begin{align}
186 &\qquad a_+\frac{4}{\Delta x}(\bar{\nu}H)|_{i+\frac{1}{2}}^j\left[\frac{\partial u}{\partial x}_{i+\frac{1}{2}}^j + \frac{1}{8}\left( b_{+E}\frac{\partial v}{\partial y}_{i+1}^{j+\frac{1}{2}}+b_{-E}\frac{\partial v}{\partial y}_{i+1}^{j-\frac{1}{2}}+b_+\frac{\partial v}{\partial y}_{i}^{j+\frac{1}{2}}+b_-\frac{\partial v}{\partial y}_{i}^{j-\frac{1}{2}} \right) \right] \nonumber \\
187 - &\qquad a_-\frac{4}{\Delta x}(\bar{\nu}H)|_{i-\frac{1}{2}}^j\left[\frac{\partial u}{\partial x}_{i-\frac{1}{2}}^j + \frac{1}{8}\left( b_+\frac{\partial v}{\partial y}_{i}^{j+\frac{1}{2}} + b_-\frac{\partial v}{\partial y}_{i}^{j-\frac{1}{2}} + b_{+W}\frac{\partial v}{\partial y}_{i-1}^{j+\frac{1}{2}} + b_{-W}\frac{\partial v}{\partial y}_{i-1}^{j-\frac{1}{2}} \right) \right] \nonumber \\
188 + &\qquad b_+\frac{1}{\Delta y}(\bar{\nu}H)|_i^{j+\frac{1}{2}}\left[\frac{\partial u}{\partial y}_i^{j+\frac{1}{2}}+\frac{1}{4}\left( a_{+N}\frac{\partial v}{\partial x}_{i+\frac{1}{2}}^{j+1} + a_+\frac{\partial v}{\partial x}_{i+\frac{1}{2}}^{j} + a_{-N}\frac{\partial v}{\partial x}_{i-\frac{1}{2}}^{j+1} + a_-\frac{\partial v}{\partial x}_{i-\frac{1}{2}}^{j} \right) \right] \nonumber \\
189 - &\qquad b_-\frac{1}{\Delta y} (\bar{\nu}H)|_i^{j-\frac{1}{2}} \left[\frac{\partial u}{\partial y}_i^{j-\frac{1}{2}}+\frac{1}{4}\left( a_+\frac{\partial v}{\partial x}_{i+\frac{1}{2}}^{j} + a_{+S}\frac{\partial v}{\partial x}_{i+\frac{1}{2}}^{j-1} + a_-\frac{\partial v}{\partial x}_{i-\frac{1}{2}}^{j} + a_{-S}\frac{\partial v}{\partial x}_{i-\frac{1}{2}}^{j-1} \right) \right] \nonumber \\
190 %%
191 = &\qquad \rho gH \frac{\partial h}{\partial x} + (a_+ - a_-)\frac{2}{\Delta x}\tau_{ocean}. \label{SSA1_dis3}
192 \end{align}
193 Again, the coefficients $a_{\pm N,S}$, $b_{\pm E,W}$ etc. are determined by the nature of the boundary between the two points used to calculate the differential quotient: zero if it is a calving front boundary and one if it is an ice-ice boundary. Thus an ice-inward scheme is used to calculate gradients, justified by assuming that the only influence of the ocean on the shelf is by exerting the extra pressure introduced on the right hand side of equation \eqref{SSA1_dis3}.\newline
194 \newline
195 With $c_{01}=(\bar{\nu}H)|_{i+\frac{1}{2}}^j, c_{00}=(\bar{\nu}H)|_{i-\frac{1}{2}}^j, c_{11}=(\bar{\nu}H)|_{i}^{j+\frac{1}{2}}, c_{10}=(\bar{\nu}H)|_{i}^{j-\frac{1}{2}}$ and if not both boundary sides exist at the same time ($a_{\pm}=0$), we can reorganize \eqref{SSA1_dis3} for each velocity contributions
196 \begin{align}
197 & u_{i}^{j}\left(\frac{-4a_+c_{01}-4a_-c_{00}}{\Delta x^2}+\frac{-b_+c_{11}-b_-c_{10}}{\Delta y^2}\right)+ \nonumber \\
198 & u_{i+1}^{j}\left(\frac{4a_+c_{01}}{\Delta x^2}\right)+u_{i-1}^{j}\left(\frac{4a_-c_{00}}{\Delta x^2}\right)+u_{i}^{j+1}\left(\frac{b_+c_{11}}{\Delta y^2}\right)+u_{i}^{j-1}\left(\frac{b_-c_{10}}{\Delta y^2}\right)+ \nonumber \\
199 & v_{i}^{j}\left(\frac{-(b_+-b_-)(a_+ c_{01}-a_- c_{00})}{2 \Delta x \Delta y}+\frac{-(a_+-a_-)(b_+ c_{11}-b_- c_{10})}{4 \Delta x \Delta y}\right)+ \nonumber \\
200 & v_{i+1}^{j}\left(\frac{-a_+ c_{01} (b_{+E}-b_{-E})}{2 \Delta x \Delta y}+\frac{a_+ (b_{+}c_{11}-b_{-}c_{10})}{4 \Delta x \Delta y}\right)+ \nonumber \\
201 & v_{i-1}^{j}\left(\frac{a_- c_{00} (b_{+W}-b_{-W})}{2 \Delta x \Delta y}+\frac{-a_-(b_{+}c_{11}-b_{-}c_{10})}{4 \Delta x \Delta y}\right)+ \nonumber \\
202 & v_{i}^{j+1}\left(\frac{-b_+ c_{11} (a_{+N}-a_{-N})}{4 \Delta x \Delta y}+\frac{b_+(a_{+}c_{01}-a_{-}c_{00})}{2 \Delta x \Delta y}\right)+ \nonumber \\
203 & v_{i}^{j-1}\left(\frac{b_- c_{10} (a_{+S}-a_{-S})}{4 \Delta x \Delta y}+\frac{-b_-(a_{+}c_{01}-a_{-}c_{00})}{2 \Delta x \Delta y}\right)+ \nonumber \\
204 & v_{i+1}^{j+1}\left(\frac{a_+ b_{+E} c_{01}}{2 \Delta x \Delta y}+\frac{a_{+N} b_{+} c_{11}}{4 \Delta x \Delta y}\right)+ v_{i-1}^{j+1}\left(\frac{-a_- b_{+W} c_{00}}{2 \Delta x \Delta y}+\frac{-a_{-N} b_{+} c_{11}}{4 \Delta x \Delta y}\right)+ \nonumber \\
205 & v_{i+1}^{j-1}\left(\frac{-a_+ b_{+E} c_{01}}{2 \Delta x \Delta y}+\frac{-a_{+S} b_{-} c_{10}}{4 \Delta x \Delta y}\right)+ v_{i-1}^{j-1}\left(\frac{a_+ b_{-W} c_{00}}{2 \Delta x \Delta y}+\frac{a_{-S} b_{-} c_{10}}{4 \Delta x \Delta y}\right)+ \nonumber \\
206 %%
207 = &\qquad \rho gH_{i}^{j} \left(\frac{a_+(h_{i+1}^{j}-h_{i}^{j})+a_-(h_{i}^{j}-h_{i-1}^{j})}{(a_++a_-)\Delta x}\right) + (a_+ - a_-)\frac{2}{\Delta x}\tau_{ocean}. \label{SSA1_dis4}
208 \end{align}
209 The coefficients enter the matrix to calculate SSA-velocities for a given right-hand side of the equation. The second SSA-equation \eqref{SSA2} is discretized accordingly.
210 The beauty is, that if there is any boundary between the neighbor cells, the stencil is identical with the default stencil used for the interior of the SSA-domain.\newline
211 However, as mentioned in \cite{Winkelmann_Martin10}, there is an inconsistency regarding the mixed derivatives. For the example of Fig.~\ref{result3}, we get for the left hand side of the grid cell $a_-=b_{+W}=b_+=1, b_-=b_{-W}=0$ and hence the second term of Eq.~\eqref{SSA1_dis3} turns into
212 \begin{align}
213 -\frac{2}{\Delta x}(\bar{\nu}H)|_{i-\frac{1}{2}}^j\left[2\frac{\partial u}{\partial x}_{i-\frac{1}{2}}^j + \frac{1}{4}\left( \frac{\partial v}{\partial y}_{i}^{j+\frac{1}{2}} + 0 + \frac{\partial v}{\partial y}_{i-1}^{j+\frac{1}{2}} + 0 \right) \right].
214 \label{SSA1_dis3b}
215 \end{align}
216 %\begin{comment}
217 The part in round parantheses means an average over the adjacent derivatives $\frac{\partial v}{\partial y}$. If one of theses derivatives does not exist due to a boundary between adjacent neighbors, it will vanish and not substituted by a pressure term on the right-hand side. It would be actually more reasonable to calculate the average by dividing by $(b_{+W}+b_++b_-+b_{-W}>0)$ instead of dividing by a constant $4$, and we would get
218 \begin{align}
219 -\frac{2}{\Delta x}(\bar{\nu}H)|_{i-\frac{1}{2}}^j\left[2\frac{\partial u}{\partial x}_{i-\frac{1}{2}}^j + \frac{1}{2}\left( \frac{\partial v}{\partial y}_{i}^{j+\frac{1}{2}} + \frac{\partial v}{\partial y}_{i-1}^{j+\frac{1}{2}}\right) \right].
220 \label{SSA1_dis3c}
221 \end{align}
222
223 %\end{comment}
224
225 \begin{thebibliography}{9}
226
227 \bibitem{Bueler2009a} Bueler, E., and J. Brown. “Shallow shelf approximation as a ‘sliding law’ in a thermomechanically coupled ice sheet model.” Journal of Geophysical Research 114, no. 3 (2009): F03008.
228
229 \bibitem{Bueler2011}Bueler, Ed, Constantine Khroulev, Andy Aschwanden, Jed Brown, and Nathan Shemonski. “Documentation for PISM, a Parallel Ice Sheet Model”, 2011. http://www.pism-docs.org/wiki/doku.php.
230
231 \bibitem{Albrecht_Martin10} Albrecht, T., M. Martin, M. Haseloff, R. Winkelmann, and A. Levermann. “Parameterization for subgrid-scale motion of ice-shelf calving fronts.” The Cryosphere 5 (2011): 35-44.
232
233 \bibitem{Winkelmann_Martin10} R. Winkelmann, M. A. Martin, M. Haseloff, T. Albrecht, E. Bueler, C. Khroulev, A. Levermann
234 The Potsdam Parallel Ice Sheet Model (PISM-PIK), Part I: Model Description
235 (2010), The Cryosphere Discussions 4, 1277-1306, DOI:10.5194/tcd-4-1277-2010.
236
237 \bibitem{Martin_Winkelmann10} M. A. Martin, R. Winkelmann, M. Haseloff, T. Albrecht, E. Bueler, C. Khroulev, A. Levermann
238 The Potsdam Parallel Ice Sheet Model (PISM-PIK), Part II: Dynamical Equilibrium Simulation of the Antarctic Ice Sheet
239 (2010), The Cryosphere Discussions 4, 1307-1341, DOI:10.5194/tcd-4-1307-2010.
240
241 \bibitem{Morland87} Morland, L. W. “Unconfined ice-shelf flow.” In Dynamics of the West Antarctic Ice Sheet: Proceedings of a Workshop held in Utrecht, May 6–8, 1987, 99–116, 1987.
242
243 %\bibitem{MacAyeal96} MacAyeal, D. R., V. Rommelaere, P. Huybrechts, C. L. Hulbe, J. Determan, and C. Ritz. “An ice-shelf model test based on the Ross Ice Shell, Antarctica.” Annals of Glaciology 23 (1996): 46–51.
244
245 \bibitem{Weis01} Weis, M. “Theory and finite element analysis of shallow ice shelves” (2001).
246
247
248 \end{thebibliography}
249
250