Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.2-d6b3a29ca committed by Constantine Khrulev on 2025-03-28
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
utilities.cc
Go to the documentation of this file.
1/* Copyright (C) 2016, 2017, 2018, 2019, 2020, 2022, 2023, 2024 PISM Authors
2 *
3 * This file is part of PISM.
4 *
5 * PISM is free software; you can redistribute it and/or modify it under the
6 * terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 3 of the License, or (at your option) any later
8 * version.
9 *
10 * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 * details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with PISM; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20#include "pism/energy/utilities.hh"
21
22#include "pism/energy/bootstrapping.hh"
23#include "pism/util/ConfigInterface.hh"
24#include "pism/util/Context.hh"
25#include "pism/util/EnthalpyConverter.hh"
26#include "pism/util/Grid.hh"
27#include "pism/util/Logger.hh"
28#include "pism/util/VariableMetadata.hh"
29#include "pism/util/array/Array3D.hh"
30#include "pism/util/array/Scalar.hh"
31#include "pism/util/error_handling.hh"
32#include "pism/util/pism_utilities.hh"
33
34namespace pism {
35namespace energy {
36
37//! Compute ice enthalpy from temperature temperature by assuming the ice has zero liquid fraction.
38/*!
39First this method makes sure the temperatures is at most the pressure-melting
40value, before computing the enthalpy for that temperature, using zero liquid
41fraction.
42
43Because of how EnthalpyConverter::pressure() works, the energy
44content in the air is set to the value that ice would have if it a chunk of it
45occupied the air; the atmosphere actually has much lower energy content. It is
46done this way for regularity (i.e. dEnth/dz computations).
47*/
48void compute_enthalpy_cold(const array::Array3D &temperature, const array::Scalar &ice_thickness,
49 array::Array3D &result) {
50
51 auto grid = result.grid();
52 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
53
54 array::AccessScope list{ &temperature, &result, &ice_thickness };
55
56 const unsigned int Mz = grid->Mz();
57 const std::vector<double> &z = grid->z();
58
59 for (auto p = grid->points(); p; p.next()) {
60 const int i = p.i(), j = p.j();
61
62 const double *Tij = temperature.get_column(i, j);
63 double *Enthij = result.get_column(i, j);
64
65 for (unsigned int k = 0; k < Mz; ++k) {
66 const double depth = ice_thickness(i, j) - z[k]; // FIXME issue #15
67 Enthij[k] = EC->enthalpy_permissive(Tij[k], 0.0, EC->pressure(depth));
68 }
69 }
70
71 result.inc_state_counter();
72
73 result.update_ghosts();
74}
75
76void compute_temperature(const array::Array3D &enthalpy, const array::Scalar &ice_thickness,
77 array::Array3D &result) {
78
79 auto grid = result.grid();
80 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
81
82 array::AccessScope list{ &enthalpy, &ice_thickness, &result };
83
84 const unsigned int Mz = grid->Mz();
85 const std::vector<double> &z = grid->z();
86
87 for (auto p = grid->points(); p; p.next()) {
88 const int i = p.i(), j = p.j();
89
90 const double *E = enthalpy.get_column(i, j), H = ice_thickness(i, j);
91 double *T = result.get_column(i, j);
92
93 for (unsigned int k = 0; k < Mz; ++k) {
94 const double depth = H - z[k]; // FIXME issue #15
95 T[k] = EC->temperature(E[k], EC->pressure(depth));
96 }
97 }
98
99 result.inc_state_counter();
100
101 result.update_ghosts();
102}
103
104//! Compute `result` (enthalpy) from `temperature` and liquid fraction.
105void compute_enthalpy(const array::Array3D &temperature,
106 const array::Array3D &liquid_water_fraction,
107 const array::Scalar &ice_thickness, array::Array3D &result) {
108
109 auto grid = result.grid();
110 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
111
112 array::AccessScope list{ &temperature, &liquid_water_fraction, &ice_thickness, &result };
113
114 const unsigned int Mz = grid->Mz();
115 const std::vector<double> &z = grid->z();
116
117 for (auto p = grid->points(); p; p.next()) {
118 const int i = p.i(), j = p.j();
119
120 const double *T = temperature.get_column(i, j);
121 const double *omega = liquid_water_fraction.get_column(i, j);
122 double *E = result.get_column(i, j);
123
124 for (unsigned int k = 0; k < Mz; ++k) {
125 const double depth = ice_thickness(i, j) - z[k]; // FIXME issue #15
126 E[k] = EC->enthalpy_permissive(T[k], omega[k], EC->pressure(depth));
127 }
128 }
129
130 result.update_ghosts();
131
132 result.inc_state_counter();
133}
134
135//! Compute the liquid fraction corresponding to enthalpy and ice_thickness.
137 const array::Scalar &ice_thickness, array::Array3D &result) {
138
139 auto grid = result.grid();
140
141 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
142
143 result.set_name("liqfrac");
144 result.metadata(0).set_name("liqfrac");
145 result.metadata(0)
146 .long_name("liquid water fraction in ice (between 0 and 1)")
147 .units("1");
148
149 array::AccessScope list{ &result, &enthalpy, &ice_thickness };
150
151 ParallelSection loop(grid->com);
152 try {
153 for (auto p = grid->points(); p; p.next()) {
154 const int i = p.i(), j = p.j();
155
156 const double *Enthij = enthalpy.get_column(i, j);
157 double *omegaij = result.get_column(i, j);
158
159 for (unsigned int k = 0; k < grid->Mz(); ++k) {
160 const double depth = ice_thickness(i, j) - grid->z(k); // FIXME issue #15
161 omegaij[k] = EC->water_fraction(Enthij[k], EC->pressure(depth));
162 }
163 }
164 } catch (...) {
165 loop.failed();
166 }
167 loop.check();
168
169 result.inc_state_counter();
170}
171
172//! @brief Compute the CTS field, CTS = E/E_s(p), from `ice_enthalpy` and `ice_thickness`, and put
173//! in `result`.
174/*!
175 * The actual cold-temperate transition surface (CTS) is the level set CTS = 1.
176 *
177 * Does not communicate ghosts for array::Array3D result.
178 */
179void compute_cts(const array::Array3D &ice_enthalpy, const array::Scalar &ice_thickness,
180 array::Array3D &result) {
181
182 auto grid = result.grid();
183 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
184
185 result.set_name("cts");
186 result.metadata(0).set_name("cts");
187 result.metadata(0)
188 .long_name("cts = E/E_s(p), so cold-temperate transition surface is at cts = 1")
189 .units("1");
190
191 array::AccessScope list{&ice_enthalpy, &ice_thickness, &result};
192
193 const unsigned int Mz = grid->Mz();
194 const std::vector<double> &z = grid->z();
195
196 for (auto p = grid->points(); p; p.next()) {
197 const int i = p.i(), j = p.j();
198
199 double *CTS = result.get_column(i,j);
200 const double *enthalpy = ice_enthalpy.get_column(i,j);
201
202 for (unsigned int k = 0; k < Mz; ++k) {
203 const double depth = ice_thickness(i,j) - z[k]; // FIXME issue #15
204 CTS[k] = enthalpy[k] / EC->enthalpy_cts(EC->pressure(depth));
205 }
206 }
207
208 result.inc_state_counter();
209}
210
211//! Computes the total ice enthalpy in J.
212/*!
213 Units of the specific enthalpy field are J kg^-1. We integrate
214 \f$E(t,x,y,z)\f$ over the entire ice fluid region \f$\Omega(t)\f$, multiplying
215 by the density to get units of energy:
216 \f[ E_{\text{total}}(t) = \int_{\Omega(t)} E(t,x,y,z) \rho_i \,dx\,dy\,dz. \f]
217*/
218double total_ice_enthalpy(double thickness_threshold,
219 const array::Array3D &ice_enthalpy,
220 const array::Scalar &ice_thickness) {
221 double enthalpy_sum = 0.0;
222
223 auto grid = ice_enthalpy.grid();
224 Config::ConstPtr config = grid->ctx()->config();
225
226 auto cell_area = grid->cell_area();
227
228 const std::vector<double> &z = grid->z();
229
230 array::AccessScope list{&ice_enthalpy, &ice_thickness};
231 ParallelSection loop(grid->com);
232 try {
233 for (auto p = grid->points(); p; p.next()) {
234 const int i = p.i(), j = p.j();
235
236 const double H = ice_thickness(i, j);
237
238 if (H >= thickness_threshold) {
239 const int ks = grid->kBelowHeight(H);
240
241 const double
242 *E = ice_enthalpy.get_column(i, j);
243
244 for (int k = 0; k < ks; ++k) {
245 enthalpy_sum += cell_area * E[k] * (z[k+1] - z[k]);
246 }
247 enthalpy_sum += cell_area * E[ks] * (H - z[ks]);
248 }
249 }
250 } catch (...) {
251 loop.failed();
252 }
253 loop.check();
254
255 enthalpy_sum *= config->get_number("constants.ice.density");
256
257 return GlobalSum(grid->com, enthalpy_sum);
258}
259
260//! Create a temperature field within the ice from provided ice thickness, surface temperature, surface mass balance, and geothermal flux.
261/*!
262In bootstrapping we need to determine initial values for the temperature within
263the ice (and the bedrock). There are various data available at bootstrapping,
264but not the 3D temperature field needed as initial values for the temperature. Here
265we take a "guess" based on an assumption of steady state and a simple model of
266the vertical velocity in the column. The rule is certainly heuristic but it
267seems to work well anyway.
268
269The result is *not* the temperature field which is in steady state with the ice
270dynamics. Spinup is most-definitely needed in many applications. Such spinup
271usually starts from the temperature field computed by this procedure and then
272runs for a long time (e.g. \f$10^4\f$ to \f$10^6\f$ years), possibly with fixed
273geometry, to get closer to thermomechanically-coupled equilibrium.
274
275Consider a horizontal grid point. Suppose the surface temperature
276\f$T_s\f$, surface mass balance \f$m\f$, and geothermal flux \f$g\f$ are given at that location.
277Within the column denote the temperature by \f$T(z)\f$ at height \f$z\f$ above
278the base of the ice. Suppose the column of ice has height \f$H\f$, the ice
279thickness.
280
281There are two alternative bootstrap methods determined by the configuration parameter
282`config.get_number("bootstrapping.temperature_heuristic"))`. Allowed values are `"smb"` and
283`"quartic_guess"`.
284
2851. If the `smb` method is chosen, which is the default, and if \f$m>0\f$,
286then the method sets the ice
287temperature to the solution of the steady problem [\ref Paterson]
288 \f[\rho_i c w \frac{\partial T}{\partial z} = k_i \frac{\partial^2 T}{\partial z^2} \qquad \text{with boundary conditions} \qquad T(H) = T_s \quad \text{and} \quad \frac{\partial T}{\partial z}(0) = - \frac{g}{k_i}, \f]
289where the vertical velocity is linear between the surface value \f$w=-m\f$ and
290a velocity of zero at the base:
291 \f[w(z) = - m z / H.\f]
292(Note that because \f$m>0\f$, this vertical velocity is downward.)
293This is a two-point boundary value problem for a linear ODE. In fact, if
294\f$K = k_i / (\rho_i c)\f$ then we can write the ODE as
295 \f[K T'' + \frac{m z}{H} T' = 0.\f]
296Then let
297 \f[C_0 = \frac{g \sqrt{\pi H K}}{k_i \sqrt{2 m}}, \qquad \gamma_0 = \sqrt{\frac{mH}{2K}}.\f]
298(Note \f$\gamma_0\f$ is, up to a constant, the square root of the Peclet number
299[\ref Paterson]; compare [\ref vanderWeletal2013].) The solution to the
300two-point boundary value problem is then
301 \f[T(z) = T_s + C_0 \left(\operatorname{erf}(\gamma_0) - \operatorname{erf}\left(\gamma_0 \frac{z}{H}\right)\right).\f]
302If `usesmb` is true and \f$m \le 0\f$, then the velocity in the column, relative
303to the base, is taken to be zero. Thus the solution is
304 \f[ T(z) = \frac{g}{k_i} \left( H - z \right) + T_s, \f]
305a straight line whose slope is determined by the geothermal flux and whose value
306at the ice surface is the surface temperature, \f$T(H) = T_s\f$.
3072. If the `quartic_guess` method is chosen, the "quartic guess" formula which was in older
308versions of PISM is used. Namely, within the ice we set
309\f[T(z) = T_s + \alpha (H-z)^2 + \beta (H-z)^4\f]
310where \f$\alpha,\beta\f$ are chosen so that
311\f[\frac{\partial T}{\partial z}\Big|_{z=0} = - \frac{g}{k_i} \qquad \text{and} \qquad \frac{\partial T}{\partial z}\Big|_{z=H/4} = - \frac{g}{2 k_i}.\f]
312The purpose of the second condition is that when ice is advecting downward then
313the temperature gradient is much larger in roughly the bottom quarter of the
314ice column. However, without the surface mass balance, much less the solution
315of the stress balance equations, we cannot estimate the vertical velocity, so
316we make such a rough guess.
317
318In either case the temperature within the ice is not allowed to exceed the
319pressure-melting temperature.
320
321We set \f$T(z)=T_s\f$ above the top of the ice.
322
323This method determines \f$T(0)\f$, the ice temperature at the ice base. This
324temperature is used by BedThermalUnit::bootstrap() to determine a
325bootstrap temperature profile in the bedrock.
326*/
327void bootstrap_ice_temperature(const array::Scalar &ice_thickness,
328 const array::Scalar &ice_surface_temp,
329 const array::Scalar &surface_mass_balance,
330 const array::Scalar &basal_heat_flux,
331 array::Array3D &result) {
332
333 auto grid = result.grid();
334 auto ctx = grid->ctx();
335 auto config = ctx->config();
336 auto log = ctx->log();
337 auto EC = ctx->enthalpy_converter();
338
339 auto use_smb = config->get_string("bootstrapping.temperature_heuristic") == "smb";
340
341 if (use_smb) {
342 log->message(2,
343 " - Filling 3D ice temperatures using surface temperature"
344 " (and mass balance for velocity estimate)...\n");
345
346 } else {
347 log->message(2,
348 " - Filling 3D ice temperatures using surface temperature"
349 " (and a quartic guess without SMB)...\n");
350 }
351
352 const double
353 ice_k = config->get_number("constants.ice.thermal_conductivity"),
354 ice_density = config->get_number("constants.ice.density"),
355 ice_c = config->get_number("constants.ice.specific_heat_capacity"),
356 K = ice_k / (ice_density * ice_c),
357 T_min = config->get_number("energy.minimum_allowed_temperature"),
358 T_melting = config->get_number("constants.fresh_water.melting_point_temperature",
359 "kelvin");
360
361 array::AccessScope list{&ice_surface_temp, &surface_mass_balance,
362 &ice_thickness, &basal_heat_flux, &result};
363
364 ParallelSection loop(grid->com);
365 try {
366 for (auto p = grid->points(); p; p.next()) {
367 const int i = p.i(), j = p.j();
368
369 const double
370 T_surface = std::min(ice_surface_temp(i, j), T_melting),
371 H = ice_thickness(i, j),
372 G = basal_heat_flux(i, j);
373
374 const unsigned int ks = grid->kBelowHeight(H);
375
376 if (G < 0.0 and ks > 0) {
377 std::string units = basal_heat_flux.metadata()["units"];
378 int Mbz = config->get_number("grid.Mbz");
379 const char *quantity = (Mbz > 0 ?
380 "temperature of the bedrock thermal layer" :
381 "geothermal flux");
383 "Negative upward heat flux (%f %s) through the bottom of the ice column\n"
384 "is not allowed by PISM's ice temperature bootstrapping method.\n"
385 "Please check the %s at i=%d, j=%d.",
386 G, units.c_str(), quantity, i, j);
387 }
388
389 if (T_surface < T_min) {
391 "T_surface(%d,%d) = %f < T_min = %f kelvin",
392 i, j, T_surface, T_min);
393 }
394
395 double *T = result.get_column(i, j);
396
397 // within ice
398 if (use_smb) { // method 1: includes surface mass balance in estimate
399
400 // Convert SMB from "kg m-2 s-1" to "m second-1".
401 const double SMB = surface_mass_balance(i, j) / ice_density;
402
403 for (unsigned int k = 0; k < ks; k++) {
404 const double z = grid->z(k);
405 T[k] = ice_temperature_guess_smb(EC, H, z, T_surface, G, ice_k, K, SMB);
406 }
407
408 } else { // method 2: a quartic guess; does not use SMB
409
410 for (unsigned int k = 0; k < ks; k++) {
411 const double z = grid->z(k);
412 T[k] = ice_temperature_guess(EC, H, z, T_surface, G, ice_k);
413 }
414
415 }
416
417 // Make sure that resulting temperatures are not too low.
418 for (unsigned int k = 0; k < ks; k++) {
419 if (T[k] < T_min) {
421 "T(%d,%d,%d) = %f < T_min = %f kelvin",
422 i, j, k, T[k], T_min);
423 }
424 }
425
426 // above ice
427 for (unsigned int k = ks; k < grid->Mz(); k++) {
428 T[k] = T_surface;
429 }
430 }
431 } catch (...) {
432 loop.failed();
433 }
434 loop.check();
435
436 result.update_ghosts();
437}
438
439void bootstrap_ice_enthalpy(const array::Scalar &ice_thickness,
440 const array::Scalar &ice_surface_temp,
441 const array::Scalar &surface_mass_balance,
442 const array::Scalar &basal_heat_flux,
443 array::Array3D &result) {
444
445 bootstrap_ice_temperature(ice_thickness, ice_surface_temp,
446 surface_mass_balance, basal_heat_flux,
447 result);
448
449 compute_enthalpy_cold(result, ice_thickness, result);
450}
451
452} // end of namespace energy
453} // end of namespace pism
std::shared_ptr< const Config > ConstPtr
std::shared_ptr< EnthalpyConverter > Ptr
void failed()
Indicates a failure of a parallel section.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
VariableMetadata & set_name(const std::string &name)
double get_number(const std::string &name) const
Get a single-valued scalar attribute.
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition Array.hh:64
double * get_column(int i, int j)
Definition Array3D.cc:121
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition Array3D.hh:33
void set_name(const std::string &name)
Sets the variable name to name.
Definition Array.cc:342
std::shared_ptr< const Grid > grid() const
Definition Array.cc:131
void inc_state_counter()
Increment the object state counter.
Definition Array.cc:150
void update_ghosts()
Updates ghost points.
Definition Array.cc:615
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition Array.cc:476
#define PISM_ERROR_LOCATION
const double G
Definition exactTestK.c:40
void compute_liquid_water_fraction(const array::Array3D &enthalpy, const array::Scalar &ice_thickness, array::Array3D &result)
Compute the liquid fraction corresponding to enthalpy and ice_thickness.
Definition utilities.cc:136
void bootstrap_ice_temperature(const array::Scalar &ice_thickness, const array::Scalar &ice_surface_temp, const array::Scalar &surface_mass_balance, const array::Scalar &basal_heat_flux, array::Array3D &result)
Create a temperature field within the ice from provided ice thickness, surface temperature,...
Definition utilities.cc:327
void compute_cts(const array::Array3D &ice_enthalpy, const array::Scalar &ice_thickness, array::Array3D &result)
Compute the CTS field, CTS = E/E_s(p), from ice_enthalpy and ice_thickness, and put in result.
Definition utilities.cc:179
double ice_temperature_guess(EnthalpyConverter::Ptr EC, double H, double z, double T_surface, double G, double ice_k)
double ice_temperature_guess_smb(EnthalpyConverter::Ptr EC, double H, double z, double T_surface, double G, double ice_k, double K, double SMB)
void compute_enthalpy_cold(const array::Array3D &temperature, const array::Scalar &ice_thickness, array::Array3D &result)
Compute ice enthalpy from temperature temperature by assuming the ice has zero liquid fraction.
Definition utilities.cc:48
void compute_temperature(const array::Array3D &enthalpy, const array::Scalar &ice_thickness, array::Array3D &result)
Definition utilities.cc:76
void bootstrap_ice_enthalpy(const array::Scalar &ice_thickness, const array::Scalar &ice_surface_temp, const array::Scalar &surface_mass_balance, const array::Scalar &basal_heat_flux, array::Array3D &result)
Definition utilities.cc:439
void compute_enthalpy(const array::Array3D &temperature, const array::Array3D &liquid_water_fraction, const array::Scalar &ice_thickness, array::Array3D &result)
Compute result (enthalpy) from temperature and liquid fraction.
Definition utilities.cc:105
double total_ice_enthalpy(double thickness_threshold, const array::Array3D &ice_enthalpy, const array::Scalar &ice_thickness)
Computes the total ice enthalpy in J.
Definition utilities.cc:218
static const double k
Definition exactTestP.cc:42
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)