PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
exactTestO.c
Go to the documentation of this file.
1 /*
2  Copyright (C) 2011, 2013, 2014, 2016, 2023 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 "pism/verification/tests/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 }
static const double L
Definition: exactTestO.c:25
static const double rho_ICE
Definition: exactTestO.c:27
static const double T_triple
Definition: exactTestO.c:24
static const double grav
Definition: exactTestO.c:26
static const double k_ICE
Definition: exactTestO.c:28
static const double H0
Definition: exactTestO.c:30
static const double beta_CC
Definition: exactTestO.c:23
int exactO_old(double z, double *TT, double *Tm, double *qice, double *qbed, double *bmelt)
Implements an exact solution for basal melt rate. Utterly straightforward arithmetic.
Definition: exactTestO.c:84
struct TestOParameters exactO(double z)
Definition: exactTestO.c:112
static const double G
Definition: exactTestO.c:32
static const double Ts
Definition: exactTestO.c:31
static const double k_BEDROCK
Definition: exactTestO.c:29