tdiscretization.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
---
tdiscretization.tex (6864B)
---
1 \documentclass{article}
2 \usepackage{amsmath}
3 \usepackage{hyperref}
4 \usepackage[left=0.5in,right=0.5in,bottom=1in]{geometry}
5 \usepackage{fancyhdr}
6 \parindent=0in \parskip=0.5\baselineskip
7
8 \include{formulas}
9
10 \newcommand{\diff}[2]{\frac{\partial #1}{\partial #2}}
11 \newcommand{\wcase}[3]{\begin{cases} #1, & w^{n}_{#3} < 0,\\#2, & w^{n}_{#3} \ge 0.\end{cases}}
12 \newcommand{\wcaseeq}[4]{#1 &= \wcase{#2}{#3}{#4}}
13
14 \begin{document}
15 \pagestyle{fancy}
16 \fancyhead[L]{Discretization of the vertical problem for the conservation of energy}
17 \fancyhead[R]{Constantine Khroulev}
18
19 \section{Introduction}
20 \label{sec:introduction}
21
22 These notes dress up \textsc{Maxima}-generated formulas needed to
23 implement (and check the implementation of) the energy balance code in
24 PISM, including the implementation of boundary conditions.
25
26 \begin{center}
27 \begin{tabular}{ll}
28 Symbol or formula & Comment \\
29 \hline
30 $\rho_{i}$ & ice density\\
31 $c_{i}$ & specific heat capacity of ice\\
32 $k_{i}$ & thermal conductivity of ice\\
33 $\phi$ & heat flux at a boundary\\
34 $G = \phi / K$ & enthalpy flux at a boundary (see also (\ref{eq:2}))\\
35 $\Delta(U_{k}) = U_{k+1} - U_{k-1}$ & centered finite difference\\
36 $\delta_{+}(U_{k}) = U_{k+1} - U_{k}$ & forward finite difference\\
37 $\delta_{-}(U_{k}) = U_{k} - U_{k-1}$ & backward finite difference\\
38 $\operatorname{Up}(U_{k}, \alpha) =
39 \begin{cases}
40 \delta_{+}(U_{k}), & \alpha < 0,\\
41 \delta_{-}(U_{k}), & \alpha \ge 0.
42 \end{cases}$ & first order upwinding\\
43 \end{tabular}
44 \end{center}
45
46 \section{Energy balance in a column of ice}
47 \label{sec:energy-in-a-column}
48
49 The conservation of energy equation is
50
51 \begin{equation}
52 \label{eq:1}
53 \rho_{i} \cdot \frac{dE}{dt} = -\nabla\cdot\mathbf{q} + Q,
54 \end{equation}
55 where
56 $\mathbf{q}$, the enthalpy flux, is defined by
57 \begin{align}
58 \label{eq:2}
59 \mathbf{q} &= -K \nabla E,\\
60 K(E) &=
61 \begin{cases}
62 k_{i}/c_{i}, & E \le E_{s}\quad\text{(cold ice)},\\
63 K_{0}, & E > E_{s}\quad\text{(temperate ice)}.\\
64 \end{cases}
65 \end{align}
66 Here $E_{s}$ is the pressure-dependent enthalpy of the cold-temperate
67 transition and $Q$ is the source term consisting of the volumetric
68 strain heating and basal frictional heating.
69
70 Note that $K$ depends on $E$, so the system defined by equations
71 (\ref{eq:1}) and (\ref{eq:2}) with appropriate boundary conditions is
72 nonlinear. The PISM implementation linearizes it by treating the time
73 derivative implicitly and using the \emph{previous} (or
74 initial) distribution of $E$ to evaluate $K(E)$.
75
76 PISM uses a shallow approximation of this equation without horizontal
77 conduction terms. Here I focus on the vertical direction, so
78 (\ref{eq:1}) becomes
79 \begin{equation}
80 \label{eq:3}
81 \rho_{i} \left( \diff{E}{t} + w\,\diff{E}{z} \right) - \diff{}{z}\left( K \diff{E}{z} \right) = \Phi.
82 \end{equation}
83 The term $\Phi$ above combines horizontal advection and the source:
84 \begin{equation}
85 \label{eq:4}
86 \Phi = Q - \rho_{i} \left( u\,\diff{E}{x} + v\,\diff{E}{y} \right).
87 \end{equation}
88
89 \section{Interior of the column}
90 \label{sec:interior}
91
92 The discretization scheme for equation (\ref{eq:3}) used in PISM is
93 \begin{equation}
94 \label{eq:5}
95 \discretization.
96 \end{equation}
97 Note that (\ref{eq:5}) uses a convex combination of centered finite
98 difference and first-order upwinding approximations of vertical
99 advection with the weight $\lambda$.
100
101 For $w \ge 0$ this can be rewritten as
102 \begin{equation}
103 \label{eq:6}
104 (\Lp)\, \El + (\Dp)\, \E_{k} + (\Up)\, \Eu = \Bp,
105 \end{equation}
106 and for $w < 0$
107 \begin{equation}
108 \label{eq:7}
109 (\Lm)\, \El + (\Dm)\, \E_{k} + (\Um)\, \Eu = \Bm.
110 \end{equation}
111
112 Here
113 \begin{equation}
114 R^{n}_{k} = \R,\quad\quad \mu = \mufactor.
115 \end{equation}
116
117 This corresponds to the following lower-diagonal ($L_{k}$), diagonal
118 ($D_{k}$) and upper-diagonal ($U_{k}$) entries in the tridiagonal
119 matrix corresponding to the discretization (\ref{eq:5}).
120
121 \begin{align}
122 \wcaseeq{L_{k}}{\Lm}{\Lp}{k}\\
123 \notag\\
124 \wcaseeq{D_{k}}{\Dm}{\Dp}{k}\\
125 \notag\\
126 \wcaseeq{U_{k}}{\Um}{\Up}{k}
127 \end{align}
128 The right hand side $b$ has elements
129 \begin{equation}
130 \label{eq:8}
131 b_{k} = \Bp
132 \end{equation}
133 in both cases.
134
135 \newcommand{\Rm}{R^{n}_{k-\frac{1}{2}}}
136 \newcommand{\Rp}{R^{n}_{k+\frac{1}{2}}}
137 \newcommand{\W}{w^{n}_{k}}
138 \newcommand{\wstack}[2]{\left\{\begin{matrix}#1\\#2\end{matrix} \right\}}
139 This can be rewritten as
140 \begin{align}
141 L_{k} &= -\Rm + \mu\W \wstack{-\lambda}{\lambda - 2},\\
142 D_{k} &= 1 + \Rm + \Rp + 2\mu\W \wstack{\lambda - 1}{1 - \lambda},\\
143 U_{k} &= -\Rp + \mu\W \wstack{2 - \lambda}{\lambda}
144 \end{align}
145 with
146 \begin{equation}
147 \label{eq:12}
148 \wstack{a}{b} = \wcase{a}{b}{}.
149 \end{equation}
150 \section{Neumann B. C.}
151 \label{sec:neumann-b-c}
152
153 In the case of the initial boundary value system consisting of
154 equation (\ref{eq:1}) combined with
155 \begin{align}
156 \label{eq:9}
157 \left.\frac{\partial E}{\partial z}\right|_{z = 0} &= G_{0}(t),\\
158 E(z_{\text{surface}}) &= f(t)
159 \end{align}
160 the Neumann boundary conditions at the base are implemented by
161 combining the generic discretization equation (\ref{eq:5}) with the
162 equation (\ref{eq:10}) approximating (\ref{eq:9}) and eliminating
163 $E_{-1}$.
164 \begin{equation}
165 \label{eq:10}
166 \neumannb.
167 \end{equation}
168 This gives
169 \begin{align}
170 \wcaseeq{D_{0}}{\Dmb}{\Dpb}{0}\\
171 \notag\\
172 \wcaseeq{U_{0}}{\Umb}{\Upb}{0}\\
173 \notag\\
174 \wcaseeq{b_{0}}{\Bmb}{\Bpb}{0}
175 \end{align}
176
177 \renewcommand{\Rp}{R^{n}_{\frac{1}{2}}}
178 \renewcommand{\Rm}{R^{n}_{-\frac{1}{2}}}
179 \renewcommand{\W}{w^{n}_{0}}
180
181 \newcommand{\rhs}[1]{E^{n}_{#1} + \frac{\Delta t \Phi^{n}_{#1}}{\rho_{i}}}
182
183 Alternatively,
184 \begin{align}
185 D_{0} &= 1 + \Rm + \Rp + 2\mu\W \wstack{\lambda - 1}{1 - \lambda}\\
186 U_{0} &= -\Rp -\Rm + 2\mu\W \wstack{1 - \lambda}{\lambda - 1}\\
187 b_{0} &= \rhs{0} + 2G_{0}\Delta z \left( -\Rm + \mu\W\wstack{-\lambda}{\lambda - 2} \right).
188 \end{align}
189
190 When assembling the basal equation we approximate $R_{-\frac{1}{2}}$ by $R_{0}$.
191
192 \newcommand{\ks}{k_{s}}
193
194 Similarly, when the Neumann B.~C. is imposed at the surface ($k =
195 \ks$), we combine (\ref{eq:5}) with (\ref{eq:11}) and eliminate $E_{\ks+1}$.
196 \begin{equation}
197 \label{eq:11}
198 \neumanns.
199 \end{equation}
200 This gives
201 \begin{align}
202 \wcaseeq{L_{\ks}}{\Lms}{\Lps}{\ks}\\
203 \notag\\
204 \wcaseeq{D_{\ks}}{\Dms}{\Dps}{\ks}\\
205 \notag\\
206 \wcaseeq{b_{\ks}}{\Bms}{\Bps}{\ks}
207 \end{align}
208
209 \renewcommand{\Rp}{R^{n}_{\ks+\frac{1}{2}}}
210 \renewcommand{\Rm}{R^{n}_{\ks-\frac{1}{2}}}
211 \renewcommand{\W}{w^{n}_{\ks}}
212
213 Alternatively,
214 \begin{align}
215 L_{\ks} &= -\Rp -\Rm + 2\mu\W \wstack{1 - \lambda}{\lambda - 1}\\
216 D_{\ks} &= 1 + \Rm + \Rp + 2\mu\W \wstack{\lambda - 1}{1 - \lambda}\\
217 b_{\ks} &= \rhs{\ks} + 2G_{\ks}\Delta z\left( \Rp + \mu\W\wstack{\lambda-2}{-\lambda} \right)
218 \end{align}
219
220 When assembling the surface equation we approximate $R_{\ks + \frac{1}{2}}$ by $R_{\ks}$.
221
222 \end{document}