tssafem_notes.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
---
tssafem_notes.tex (18439B)
---
1 \documentclass{amsart}
2 \usepackage{amsmath}
3 \usepackage{amssymb}
4 \usepackage{hyperref}
5 \usepackage{underscore}
6 \usepackage[svgnames]{xcolor}
7 \usepackage[margin=1in]{geometry}
8 \usepackage{fancybox}
9
10 % disable indenting the first line of a paragraph
11 \parindent=0in
12 \parskip=0.5\baselineskip
13
14 % show sub-sub-sections in toc
15 \setcounter{tocdepth}{4}
16
17 % trace of a matrix
18 \DeclareMathOperator{\Tr}{tr}
19
20 % allow align to span page breaks
21 \allowdisplaybreaks[1]
22
23 % vectors
24 \newcommand{\N}{\mathbf{n}}
25 \newcommand{\U}{\mathbf{u}}
26 \newcommand{\T}{\boldsymbol{\tau}}
27 % area integral
28 \newcommand{\I}{\ensuremath{\int_{\Omega}}}
29 % boundary integral
30 \newcommand{\boundary}{\partial \Omega_N}
31 \newcommand{\bI}{\ensuremath{\int_{\boundary}}}
32 % effective stress tensor
33 \newcommand{\etaM}{\ensuremath{\eta M}}
34 % partial derivatives
35 \newcommand{\diff}[2]{\ensuremath{\frac{\partial #1}{\partial #2}}}
36 % sum over quadrature points
37 \newcommand{\sumQ}{\ensuremath{\sum_{q = 1}^{N_q}}}
38 % basis expansion of a partial derivative
39 \newcommand{\diffbasisexpansion}[3]{\ensuremath{\sum_{#3 = 1}^{N_k} {#1}_{#3} \diff{\phi_{#3}}{#2} }}
40 \newcommand{\UX}{\diffbasisexpansion{u}{x}{m}}
41 \newcommand{\UY}{\diffbasisexpansion{u}{y}{m}}
42 \newcommand{\VX}{\diffbasisexpansion{v}{x}{m}}
43 \newcommand{\VY}{\diffbasisexpansion{v}{y}{m}}
44 % drag coefficient as a function of velocity
45 \newcommand{\betaU}{\beta(\U)}
46
47 % basal shear stress
48 \newcommand{\basalshearstress}[1]{\T_{b#1}}
49 \newcommand{\taub}{\basalshearstress{}}
50 \newcommand{\taubx}{\basalshearstress{,x}}
51 \newcommand{\tauby}{\basalshearstress{,y}}
52
53 % basal driving stress
54 \newcommand{\drivingstress}[1]{\T_{d#1}}
55 \newcommand{\taud}{\drivingstress{}}
56 \newcommand{\taudx}{\drivingstress{,x}}
57 \newcommand{\taudy}{\drivingstress{,y}}
58
59 \newcommand{\highlight}[1]{{\color{red!80!black} \fbox{$ \displaystyle #1 $} }}
60
61 \newcommand{\R}{\mathbb{R}}
62
63 \begin{document}
64
65 \title{SSA FEM implementation notes}
66 \author{Constantine Khroulev}
67 \maketitle
68 \tableofcontents
69
70
71 \section{Notation}
72 \label{sec-1}
73
74 \begin{center}
75 \begin{tabular}{lp{0.5\textwidth}}
76 Symbol & Meaning\\
77 \hline
78 $B$ & vertically-averaged ice hardness\\
79 $g$ & acceleration due to gravity\\
80 $H$ & ice thickness\\
81 $h$ & ice top surface elevation\\
82 $n$ & Glen flow law exponent\\
83 $N_k$ & number of trial functions per element\\
84 $N_q$ & number of quadrature points per element\\
85 $q$ & sliding power law exponent\\
86 $w_q$ & weight corresponding to a quadrature point $q$\\
87 $J_q$ & Jacobian of the map from the reference element to the physical element, evaluated at a quadrature point $q$\\
88 $\U = (u,v)$ & horizontal ice velocity\\
89 $\epsilon_{\beta}$, $\epsilon_{\nu}$, $\epsilon_{\eta}$ & regularization parameters for $\betaU$, $\nu$, $\eta$\\
90 $\nu$ & effective viscosity of ice\\
91 $\phi$ & trial functions\\
92 $\psi$ & test functions\\
93 $\rho$ & ice density\\
94 $\taub$ & basal shear stress\\
95 $\taud$ & driving shear stress\\
96 \end{tabular}
97 \end{center}
98
99 Formulas that appear in the code are \highlight{highlighted.}
100
101 \section{The shallow shelf approximation}
102 \label{sec:ssa-strong}
103
104 \newcommand{\Mux}{4u_x + 2v_y}
105 \newcommand{\Muy}{u_y + v_x}
106 \newcommand{\Mvx}{u_y + v_x}
107 \newcommand{\Mvy}{2u_x + 4v_y}
108 Define the effective SSA strain rate tensor $M$ \cite{Dukowiczetal2010}:
109 \begin{equation}
110 \label{eq:1}
111 M =
112 \begin{pmatrix}
113 \Mux & \Muy \\
114 \Mvx & \Mvy \\
115 \end{pmatrix}.
116 \end{equation}
117 Then the strong form of the SSA system (without boundary conditions) is
118 \begin{align}
119 \label{eq:2}
120 - \nabla\cdot(\eta\, M) & = \taub + \taud,\\
121 \label{eq:3}
122 \eta &= \highlight{ \epsilon_{\eta} + \nu\, H }.
123 \end{align}
124 This is equivalent to the more familiar form (found in \cite{SchoofStream}, for example):
125 \begin{align*}
126 - \Big[ (\eta (\Mux))_x + (\eta (\Muy))_y \Big] & = \taubx + \taudx, \\
127 - \Big[ (\eta (\Mvx))_x + (\eta (\Mvy))_y \Big] & = \tauby + \taudy.\\
128 \end{align*}
129
130 Here $\taud = \highlight{ \rho g H \nabla h }$ is the gravitational driving shear stress; see subsections for definitions of $\taub$ and the ice viscosity $\nu$.
131
132
133 \subsection{Ice viscosity}
134 \label{sec:ice-viscosity}
135
136 Let $U = \{u,v,w\}$ and $X = \{x, y, z\}$.
137
138 The three-dimensional strain rate tensor $D$ (\cite{GreveBlatter2009}, equations 3.25 and 3.29) is defined by
139
140 \begin{equation*}
141 D_{i,j}(\U) = \frac 12 \left( \diff{U_i}{X_j} + \diff{U_j}{X_i} \right).
142 \end{equation*}
143
144 We assume that ice is incompressible, so $w_z = - (u_x + v_y)$. Moreover, in the shallow shelf approximation horizontal velocity components do not vary with depth, so $u_z = v_z = 0$.
145
146 With these assumptions $D$ becomes
147 \begin{equation*}
148 D =
149 \begin{pmatrix}
150 u_x & \frac{1}{2}\,(u_y + v_x) & 0\\
151 \frac{1}{2}\,(u_y + v_x) & v_y & 0\\
152 0 & 0 & - (u_x + v_y)\\
153 \end{pmatrix}
154 \end{equation*}
155
156 Now the second invariant of $D$ (\cite{GreveBlatter2009}, equation 2.42)
157 \begin{equation*}
158 \gamma = \Tr (D^2) - \left( \Tr D \right)^2
159 \end{equation*}
160 simplifies to
161 \begin{equation}
162 \label{eq:4}
163 \gamma = \highlight{ \frac{1}{2}\, \left( (u_x)^2 + (v_y)^2 + (u_x + v_y)^2 + \frac{1}{2}\,(u_y + v_x)^2 \right) }.
164 \end{equation}
165
166 We define the regularized effective viscosity of ice $\nu$ (\cite{SchoofStream}, equation 2.3):
167 \begin{equation}
168 \label{eq:5}
169 \nu = \highlight{ \frac{1}{2} B \left( \epsilon_{\nu} + \gamma \right)^{(1 - n) / (2n)} },\\
170 \end{equation}
171
172 \subsection{Basal shear stress}
173 \label{sec:beta}
174
175 The basal shear stress is defined by
176 \begin{equation}
177 \label{eq:6}
178 \taub = - \betaU \U,
179 \end{equation}
180
181 where $\beta = \betaU$ is a scalar-valued drag coefficient related to the yield stress.
182
183 In PISM, $\betaU$ is defined as follows (see \cite{SchoofHindmarsh}):
184 \begin{equation}
185 \label{eq:7}
186 \betaU = \highlight{ \frac{\tau_c}{u_{\text{threshold}}^q}\cdot (\epsilon_{\beta} + |\U|^2)^{(q - 1) / 2} }
187 \end{equation}
188
189
190 \section{The weak form of the SSA}
191
192 Multiplying (\ref{eq:2}) by a test function $\psi$ and integrating by parts, we get the weak form:
193
194 \begin{align}
195 - \nabla \cdot (\etaM) & = \taud + \taub,\notag \\
196 - \I \psi \nabla \cdot (\etaM) & = \I \psi (\taud + \taub),\notag \\
197 - \I \nabla \cdot (\psi\, \etaM) + \I \nabla \psi \cdot (\etaM) & = \I \psi (\taud + \taub),\notag \\
198 - \bI (\psi\, \etaM)\cdot \N\, ds + \I \nabla \psi \cdot (\etaM) & = \I \psi (\taud + \taub)\notag \\
199 \label{eq:8}
200 \I \Big[\nabla \psi \cdot (\etaM) - \psi (\taud + \taub)\Big] - \bI (\psi\, \etaM) \cdot \N\, ds &= 0.
201 \end{align}
202
203 Recall that $\psi$ are zero on the Dirichlet part of the boundary, so the boundary integral is equal to the integral over the \emph{Neumann} part of the boundary $\boundary$.
204
205 Let $\Phi$ be the line integral above:
206 \begin{equation}
207 \label{eq:9}
208 \Phi = \bI (\psi\, \etaM) \cdot \N\, ds
209 \end{equation}
210 We keep $\Phi$ in equations below to make these notes easier to digest. See section \ref{sec:boundary} for details of the residual contribution of $\Phi$.
211
212 Define $F(u,v)$
213 \begin{align}
214 \label{eq:10}
215 F_{1}(u,v) &= \I \left[ \diff{\psi}{x} \left( \eta (\Mux) \right) + \diff{\psi}{y} \left( \eta (\Muy) \right) - \psi (\taubx + \taudx) \right] - \Phi_{1}\\
216 \label{eq:11}
217 F_{2}(u,v) &= \I \left[ \diff{\psi}{x} \left( \eta (\Mvx) \right) + \diff{\psi}{y} \left( \eta (\Mvy) \right) - \psi (\tauby + \taudy) \right] - \Phi_{2}.
218 \end{align}
219
220 These notes are concerned with solving the system $F(u,v) = 0$.
221
222 \section{Solving the discretized system}
223 \label{sec-3}
224
225 In the following subscripts $x$ and $y$ denote partial derivatives, while subscripts $k$, $l$, $m$ denote nodal values of a particular quantity. Also, to simplify notation from here on $u$, $v$, etc stand for finite element approximations of corresponding continuum variables.
226
227 To build a Galerkin approximation of the SSA system, let $\phi$ be trial functions and $\psi$ be test functions\footnote{Note that in a Galerkin approximations $\phi$ and $\psi$ are the same.}. Then we have the following basis expansions:
228
229 \begin{equation}
230 \label{eq:12}
231 \begin{aligned}
232 u & = \sum_i \phi_i\, u_i, \\
233 v & = \sum_i \phi_i\, v_i,\\
234 \diff{u}{x_j} & = \sum_i \diff{\phi_i}{x_j}\ u_i,\\
235 \diff{v}{x_j} & = \sum_i \diff{\phi_i}{x_j}\ v_i.\\
236 \end{aligned}
237 \end{equation}
238
239 We use Newton's method to solve the system resulting from discretizing equations (\ref{eq:10}) and (\ref{eq:11}). This requires computing residuals and the Jacobian matrix.
240
241 \subsection{Residual evaluation}
242 \label{sec-3-1}
243
244 In this and following sections we focus on element contributions to the residual and the Jacobian. Basis functions used here are defined on the reference element, hence the added determinant of the Jacobian of the map from the reference element to a particular physical element $|J_q|$ appearing in all quadratures.
245
246 \begin{align}
247 \label{eq:13}
248 F_{k,1} & = \highlight{ \sumQ |J_q| \cdot w_q\cdot \left[ \eta \left( \diff{\psi_k}{x} (\Mux) + \diff{\psi_k}{y} (\Muy) \right) - \psi_k (\taubx + \taudx) \right]_{\text{evaluated at } q} - \Phi_{k,1} }\\
249 \label{eq:14}
250 F_{k,2} & = \highlight{ \sumQ |J_q| \cdot w_q\cdot \left[ \eta \left( \diff{\psi_k}{x} (\Mvx) + \diff{\psi_k}{y} (\Mvy) \right) - \psi_k (\tauby + \taudy) \right]_{\text{evaluated at } q} - \Phi_{k,2}}
251 \end{align}
252
253
254 \subsection{Jacobian evaluation}
255 \label{sec-3-2}
256
257 Equations (\ref{eq:13}) and (\ref{eq:14}) define a map from $\R^{2\times N}$ to $\R^{2\times N}$, where $N$ is the number of nodes in a FEM mesh. To use Newton's method, we need to be able to compute the Jacobian \emph{of this map}.
258
259 It is helpful to rewrite equations defining $F_{k,1}$ and $F_{k,2}$ using basis expansions for $u_x$, $u_y$, $v_x$, and $v_y$ (see (\ref{eq:12})), as follows:
260
261 \begin{align*}
262 F_{k,1} &= \highlight{
263 \begin{aligned}[t]
264 \sumQ |J_q| \cdot w_q\cdot \Bigg[ &\eta \left(\diff{\psi_k}{x} \left(4 \UX + 2 \VY\right)
265 + \diff{\psi_k}{y} \left(\UY + \VX\right)\right) \\
266 & - \psi_k \left(\taubx + \taudx\right) \Bigg]_{\text{evaluated at } q} - \Phi_{k,1}\\
267 \end{aligned} }\\
268 F_{k,2} &= \highlight{
269 \begin{aligned}[t]
270 \sumQ |J_q| \cdot w_q\cdot \Bigg[ &\eta \left(\diff{\psi_k}{x} \left(\UY + \VX\right)
271 + \diff{\psi_k}{y} \left(2 \UX + 4 \VY\right) \right) \\
272 & - \psi_k \left(\tauby + \taudy\right) \Bigg]_{\text{evaluated at } q} - \Phi_{k,2}\\
273 \end{aligned}}
274 \end{align*}
275
276 The Jacobian has elements
277 \begin{align*}
278 J_{k,l,1} & = \diff{F_{k,1}}{u_l}, & J_{k,l,2} & = \diff{F_{k,1}}{v_l},\\
279 J_{k,l,3} & = \diff{F_{k,2}}{u_l}, & J_{k,l,4} & = \diff{F_{k,2}}{v_l}.\\
280 \end{align*}
281
282 \begin{align}
283 J_{k,l,1} &= \highlight{
284 \begin{aligned}[t]
285 \sumQ |J_q| \cdot w_q\cdot \Bigg[&\diff{\eta}{u_l}\cdot
286 \left( \diff{\psi_k}{x} (4 u_x + 2 v_y) + \diff{\psi_k}{y} (u_y + v_x) \right)\\
287 & + \eta\cdot \left( \diff{\psi_k}{x}\cdot 4 \diff{\phi_l}{x} + \diff{\psi_k}{y}\cdot \diff{\phi_l}{y} \right)
288 - \psi_k\cdot \diff{\taubx}{u_l} \Bigg]_{\text{evaluated at } q} - \diff{\Phi_{k,1}}{u_{l}}\\
289 \end{aligned} }\\
290 J_{k,l,2} &= \highlight{
291 \begin{aligned}[t]
292 \sumQ |J_q| \cdot w_q\cdot \Bigg[&\diff{\eta}{v_l}\cdot
293 \left( \diff{\psi_k}{x} (4 u_x + 2 v_y) + \diff{\psi_k}{y} (u_y + v_x) \right)\\
294 & + \eta\cdot \left( \diff{\psi_k}{x}\cdot 2 \diff{\phi_l}{y} + \diff{\psi_k}{y}\cdot \diff{\phi_l}{x} \right)
295 - \psi_k\cdot \diff{\taubx}{v_l} \Bigg]_{\text{evaluated at } q} - \diff{\Phi_{k,1}}{v_{l}}\\
296 \end{aligned} }\\
297 J_{k,l,3} &= \highlight{
298 \begin{aligned}[t]
299 \sumQ |J_q| \cdot w_q\cdot \Bigg[&\diff{\eta}{u_l}\cdot
300 \left( \diff{\psi_k}{x} (u_y + v_x) + \diff{\psi_k}{y} (2 u_x + 4 v_y) \right)\\
301 & + \eta\cdot \left( \diff{\psi_k}{x}\cdot \diff{\phi_l}{y} + \diff{\psi_k}{y}\cdot 2 \diff{\phi_l}{x} \right)
302 - \psi_k\cdot \diff{\tauby}{u_l} \Bigg]_{\text{evaluated at } q} - \diff{\Phi_{k,2}}{u_{l}}\\
303 \end{aligned} }\\
304 J_{k,l,4} &= \highlight{
305 \begin{aligned}[t]
306 \sumQ |J_q| \cdot w_q\cdot \Bigg[&\diff{\eta}{v_l}\cdot
307 \left( \diff{\psi_k}{x} (u_y + v_x) + \diff{\psi_k}{y} (2 u_x + 4 v_y) \right)\\
308 & + \eta\cdot \left( \diff{\psi_k}{x}\cdot \diff{\phi_l}{x} + \diff{\psi_k}{y}\cdot 4 \diff{\phi_l}{y} \right)
309 - \psi_k\cdot \diff{\tauby}{v_l} \Bigg]_{\text{evaluated at } q} - \diff{\Phi_{k,2}}{v_{l}}\\
310 \end{aligned} }
311 \end{align}
312
313 In our case the number of trial functions $N_k$ is $4$ ($Q_1$ elements). Our test functions are the same as trial functions (a Galerkin method), i.e. we also have $4$ test functions per element. Moreover, each combination of test and trial functions corresponds to $4$ values in the Jacobian ($2$ equations, $2$ degrees of freedom). Overall, each element contributes to $4 \times 4 \times 4 = 64$ entries in the Jacobian matrix.
314
315 To evaluate $J_{\cdot,\cdot,\cdot}$, we need be able to compute the following\footnote{Also $\diff{\Phi_{\cdot,\cdot}}{u}$ and $\diff{\Phi_{\cdot,\cdot}}{v}$; see section \ref{sec:boundary}.}:
316 \begin{equation*}
317 \diff{\eta}{u_l}, \quad \diff{\eta}{v_l}, \quad \diff{\taubx}{u_l},
318 \quad \diff{\taubx}{v_l}, \quad \diff{\tauby}{u_l}, \quad \diff{\tauby}{v_l}.
319 \end{equation*}
320
321 Subsections that follow describe related implementation details.
322
323
324 \subsubsection{Ice viscosity}
325 \label{sec:viscosity-evaluation}
326
327 Recall (equation (\ref{eq:3})) that $\eta = \epsilon_{\eta} + \nu H$. We use chain rule to get
328 \begin{align*}
329 \diff{\eta}{u_l} &= \highlight{ H\,\diff{\nu}{\gamma}\cdot \diff{\gamma}{u_l}, } & \diff{\eta}{v_l} &= \highlight{ H\,\diff{\nu}{\gamma}\cdot \diff{\gamma}{v_l} }.
330 \end{align*}
331
332 The derivative of $\nu$ with respect to $\gamma$ (equation (\ref{eq:4})) can be written in terms of $\nu$ itself:
333
334 \begin{align*}
335 \diff{\nu}{\gamma} & = \frac{1}{2} B \cdot \frac{1 - n}{2n} \cdot \left(\epsilon_{\nu} + \gamma \right)^{(1 - n) / (2n) - 1}, \\
336 & = \frac{1 - n}{2n} \cdot \frac{1}{2} B \left( \epsilon_{\nu} + \gamma \right)^{(1 - n) / (2n)} \cdot \frac{1}{\epsilon_{\nu} + \gamma}, \\
337 & = \highlight{ \frac{1 - n}{2n} \cdot \frac{\nu}{\epsilon_{\nu} + \gamma} }.
338 \end{align*}
339
340 To compute $\diff{\gamma}{u_l}$ and $\diff{\gamma}{v_l}$ we need to re-write $\gamma$ (equation (\ref{eq:4})) using the basis expansion (\ref{eq:12}):
341 \begin{align*}
342 \gamma &=
343 \begin{aligned}[t]
344 \frac{1}{2}\, \Bigg(&\left(\UX\right)^2 + \left(\VY\right)^2 \\
345 & + \left(\UX + \VY\right)^2 + \frac{1}{2}\, \left(\UY + \VX\right)^2\Bigg). \\
346 \end{aligned}
347 \end{align*}
348
349 So,
350 \begin{align}
351 \diff{\gamma}{u_l} &= u_x \diff{\phi_l}{x} + (u_x + v_y)\diff{\phi_l}{x} + \frac{1}{2} (u_y + v_x)\diff{\phi_l}{y},\notag\\
352 &= \highlight{ (2 u_x + v_y) \diff{\phi_l}{x} + \frac{1}{2} (u_y + v_x)\diff{\phi_l}{y} },\\
353 \diff{\gamma}{v_l} &= v_y \diff{\phi_l}{y} + (u_x + v_y)\diff{\phi_l}{y} + \frac{1}{2} (u_y + v_x)\diff{\phi_l}{x},\notag\\
354 &= \highlight{ \frac{1}{2} (u_y + v_x)\diff{\phi_l}{x} + (u_x + 2 v_y)\diff{\phi_l}{y} }.
355 \end{align}
356
357
358 \subsubsection{Basal drag}
359 \label{sec:basal-drag}
360
361 The method \texttt{IceBasalResistancePlasticLaw::drag_with_derivative()} computes $\beta$ and the derivative of $\beta$ with respect to $\alpha = \frac12 |\U|^2 = \frac12 (u^2 + v^2)$.
362
363 So,
364 \begin{align*}
365 \diff{\alpha}{u_l} &= u\cdot \diff{u}{u_l} = u\cdot \phi_l,&
366 \diff{\alpha}{v_l} &= v\cdot \diff{v}{v_l} = v\cdot \phi_l.
367 \end{align*}
368
369 Recall from equation (\ref{eq:6})
370 \begin{align*}
371 \taubx &= \highlight{ - \betaU\cdot u }, & \tauby &= \highlight{ - \betaU\cdot v }.
372 \end{align*}
373
374 Using product and chain rules, we get
375 \begin{align*}
376 \diff{\taubx}{u_l} &= -\left(\betaU\cdot \diff{u}{u_l} + \diff{\betaU}{u_l}\cdot u\right)\\
377 &= -\left( \betaU\cdot \phi_l + \diff{\betaU}{\alpha}\cdot \diff{\alpha}{u_l}\cdot u \right)\\
378 &= \highlight{ -\left( \betaU\cdot \phi_l + \diff{\betaU}{\alpha}\cdot u^2 \phi_l \right) },\\
379 \diff{\taubx}{v_l} &= -\diff{\betaU}{v_l}\cdot u\\
380 &= -\diff{\betaU}{\alpha}\cdot \diff{\alpha}{v_l}\cdot u\\
381 &= \highlight{ -\diff{\betaU}{\alpha}\cdot u\cdot v \cdot \phi_l },\\
382 \diff{\tauby}{u_l} &= -\diff{\betaU}{u_l}\cdot v,\\
383 &= \highlight{ -\diff{\betaU}{\alpha}\cdot u\cdot v \cdot \phi_l },\\
384 \diff{\tauby}{v_l} &= -\left( \betaU\cdot \diff{v}{v_l} + \diff{\betaU}{u_l}\cdot v \right)\\
385 &= -\left( \betaU\cdot \phi_l + \diff{\betaU}{\alpha}\cdot \diff{\alpha}{v_l}\cdot v \right)\\
386 &= \highlight{ -\left( \betaU\cdot \phi_l + \diff{\betaU}{\alpha}\cdot v^2 \phi_l \right) }.
387 \end{align*}
388
389 \section{Boundary contributions}
390 \label{sec:boundary}
391
392 The boundary contribution consists of the integral
393 \begin{equation}
394 \label{eq:15}
395 \Phi = \bI (\psi\, \etaM) \cdot \N\, ds,
396 \end{equation}
397 where $M$ is the effective strain rate tensor for the shallow shelf
398 approximation and $\eta = \epsilon_{\eta} + \nu\, H$.
399
400 Define
401 \newcommand{\dP}{\Delta P}
402 \begin{equation}
403 \label{eq:17}
404 \dP = \int_{b}^{h} (p_{\mathrm{ice}} - p_{\mathrm{ocean}})\,dz,
405 \end{equation}
406 the ``pressure difference'' at the ice margin. Note that $\dP$ is a function of ice geometry and sea level, but it does \emph{not} depend on the ice velocity and therefore does \emph{not} contribute to the Jacobian.
407
408 In our particular application (see \cite{GreveBlatter2009} and \cite{KhroulevBueler2013}), the calving front boundary condition states
409
410 \begin{equation}
411 \label{eq:16}
412 \etaM \cdot \N = \dP\, \N.
413 \end{equation}
414
415 Combining equations (\ref{eq:16}) and (\ref{eq:17}) gives the following expression for $\Phi$:
416 \begin{equation}
417 \label{eq:18}
418 \Phi = \bI \psi\, \dP\,ds.
419 \end{equation}
420 Note that the right-hand side is a \emph{scalar}, so equation (\ref{eq:18}) gives the expression for \emph{both} $\Phi_{1}$ and $\Phi_{2}$.
421
422 \subsection{Residual evaluation}
423 \label{sec:boundary-residual}
424
425 To evaluate the contribution of (\ref{eq:18}) to the residual we need to rewrite this line integral as a regular integral by using a parameterization of $\boundary$.
426
427 As usual, we assemble it element-by-element.
428
429 \bibliography{ssafem-notes}
430 \bibliographystyle{siam}
431
432