texactTestO.c - 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
---
texactTestO.c (5592B)
---
1 /*
2 Copyright (C) 2011, 2013, 2014, 2016 Ed Bueler
3
4 This file is part of PISM.
5
6 PISM is free software; you can redistribute it and/or modify it under the
7 terms of the GNU General Public License as published by the Free Software
8 Foundation; either version 3 of the License, or (at your option) any later
9 version.
10
11 PISM is distributed in the hope that it will be useful, but WITHOUT ANY
12 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 details.
15
16 You should have received a copy of the GNU General Public License
17 along with PISM; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21 #include "exactTestO.h"
22
23 static const double beta_CC = 7.9e-8; /* K Pa-1; Clausius-Clapeyron constant [@ref Luethi2002] */
24 static const double T_triple = 273.15; /* K; triple point of pure water */
25 static const double L = 3.34e5; /* J kg-1; latent heat of fusion for water [@ref AschwandenBlatter] */
26 static const double grav = 9.81; /* m/s^2; accel of gravity */
27 static const double rho_ICE = 910.0; /* kg/(m^3) density of ice */
28 static const double k_ICE = 2.10; /* J/(m K s) = W/(m K) thermal conductivity of ice */
29 static const double k_BEDROCK = 3.0; /* J/(m K s) = W/(m K) thermal conductivity of bedrock */
30 static const double H0 = 3000.0; /* m */
31 static const double Ts = 223.15; /* K */
32 static const double G = 0.042; /* W/(m^2) = J m-2 s-1 */
33
34 /*! \brief Implements an exact solution for basal melt rate. Utterly straightforward arithmetic. */
35 /*!
36 Assumes a steady-state temperature profile in ice and bedrock. This steady
37 profile is driven by geothermal flux `G` from the crust below the bedrock
38 thermal layer. The heat flux is everywhere upward in the bedrock thermal layer,
39 and it is equal by construction to `G`. The heat flux upward in the ice is
40 also a constant everywhere in the ice, but its value is smaller than `G`. This
41 imbalance is balanced by the basal melt rate `bmelt`.
42
43 Geometry and dynamics context: The top of the ice is flat so the ice
44 does not flow; the ice thickness has constant value `H0`. The ice surface has
45 temperature `Ts` which is below the melting point. The basal melt rate
46 `bmelt` does not contribute to the vertical velocity in the ice. The ice
47 pressure is hydrostatic: \f$p = \rho_i g (h-z)\f$.
48
49 The basic equation relates the fluxes in the basal ice, and in the top of the
50 bedrock layer, with the basal melt rate `bmelt` \f$= - M_b / \rho_i \f$. The
51 equation is from [\ref AschwandenBuelerKhroulevBlatter],
52 \f[ M_b H + (\mathbf{q} - \mathbf{q_{lith}}) \cdot \mathbf{n} = F_b + Q_b. \f]
53 Here \f$-M_b\f$ is the basal melt rate in units of mass per area per time.
54 In the above equation the basal friction heating
55 \f$F_b\f$ is zero and the subglacial aquifer enthalpy flux \f$Q_b\f$ includes no
56 horizontal velocity. (Note that \f$Q_b\f$ is the heat delivered by subglacial
57 water to the base of the ice.) We assume the subglacial water is at the ice
58 overburden pressure \f$p_0=\rho_i g H_0\f$, and we assume that the temperate
59 layer at the base of the ice has zero thickness, so \f$\omega = 0\f$. Thus
60 \f[ H_l(p_b) = H_l(p_0) = H_s(p_0) + L, \f]
61 \f[ H = H_s(p_0) + \omega L = H_s(p_0), \f]
62 \f[ Q_b = H_l(p_0) M_b. \f]
63
64 The basic equation therefore reduces to
65 \f[ \mathbf{q} \cdot \mathbf{n} - G = M_b L \f]
66 or
67 \f[ \text{\texttt{bmelt}} = -\frac{M_b}{\rho_i}
68 = \frac{G - \mathbf{q} \cdot \mathbf{n}}{L \rho_i}. \f]
69 On the other hand, the temperature in the ice is the steady-state result wherein
70 the upward flux is constant and the (absolute) temperature is a linear function
71 of the elevation \f$z\f$, so
72 \f[ \mathbf{q} \cdot \mathbf{n} = - k_i \frac{T_s - T_m(p)}{H_0}.\f]
73
74 The temperature in the ice (\f$0 \le z \le H_0\f$) is this linear function,
75 \f[ T(z) = T_m(p) + \frac{T_s - T_m(p)}{H_0} z \f]
76 and in the bedrock (\f$z \le 0\f$) is also linear,
77 \f[ T_b(z) = T_m(p) - \frac{G}{k_b} z. \f]
78
79 This method implements these formulas. It should be called both when setting-up
80 a verification test by setting temperature at different elevations within
81 the ice and bedrock, and when doing the verification itself by checking against
82 the exact `bmelt` value.
83 */
84 int exactO_old(double z, double *TT, double *Tm, double *qice, double *qbed, double *bmelt) {
85
86 double P_base;
87
88 P_base = rho_ICE * grav * H0; /* Pa; hydrostatic approximation to pressure
89 at base */
90
91 *Tm = T_triple - beta_CC * P_base;/* K; pressure-melting point at base */
92
93 *qice = - k_ICE * (Ts - *Tm) / H0;/* J m-1 K-1 s-1 K m-1 = J m-2 s-1 = W m-2;
94 equilibrium heat flux everywhere in ice */
95
96 *qbed = G; /* J m-2 s-1 = W m-2; equilibrium heat flux
97 everywhere in bedrock */
98
99 *bmelt = (G - *qice) / (L * rho_ICE);/* J m-2 s-1 / (J kg-1 kg m-3) = m s-1;
100 ice-equivalent basal melt rate */
101
102 if (z > H0) {
103 *TT = Ts; /* K; in air above ice */
104 } else if (z >= 0.0) {
105 *TT = *Tm + (Ts - *Tm) * (z / H0); /* in ice */
106 } else {
107 *TT = *Tm - (G / k_BEDROCK) * z; /* in bedrock */
108 }
109 return 0;
110 }
111
112 struct TestOParameters exactO(double z) {
113 struct TestOParameters result;
114
115 exactO_old(z, &result.TT, &result.Tm, &result.qice, &result.qbed, &result.bmelt);
116
117 return result;
118 }