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
EnthalpyModel.cc
Go to the documentation of this file.
1/* Copyright (C) 2016, 2017, 2018, 2022, 2023 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/EnthalpyModel.hh"
21#include "pism/energy/DrainageCalculator.hh"
22#include "pism/energy/enthSystem.hh"
23#include "pism/energy/utilities.hh"
24#include "pism/util/Context.hh"
25#include "pism/util/EnthalpyConverter.hh"
26#include "pism/util/array/CellType.hh"
27#include "pism/util/io/File.hh"
28
29namespace pism {
30namespace energy {
31
32EnthalpyModel::EnthalpyModel(std::shared_ptr<const Grid> grid,
33 std::shared_ptr<const stressbalance::StressBalance> stress_balance)
34 : EnergyModel(grid, stress_balance) {
35 // empty
36}
37
38void EnthalpyModel::restart_impl(const File &input_file, int record) {
39
40 m_log->message(2, "* Restarting the enthalpy-based energy balance model from %s...\n",
41 input_file.name().c_str());
42
43 m_basal_melt_rate.read(input_file, record);
44 init_enthalpy(input_file, false, record);
45
46 regrid("Energy balance model", m_basal_melt_rate, REGRID_WITHOUT_REGRID_VARS);
48}
49
50void EnthalpyModel::bootstrap_impl(const File &input_file,
51 const array::Scalar &ice_thickness,
52 const array::Scalar &surface_temperature,
53 const array::Scalar &climatic_mass_balance,
54 const array::Scalar &basal_heat_flux) {
55
56 m_log->message(2, "* Bootstrapping the enthalpy-based energy balance model from %s...\n",
57 input_file.name().c_str());
58
59 m_basal_melt_rate.regrid(input_file,
60 io::Default(m_config->get_number("bootstrapping.defaults.bmelt")));
61
62 regrid("Energy balance model", m_basal_melt_rate, REGRID_WITHOUT_REGRID_VARS);
63
64 int enthalpy_revision = m_ice_enthalpy.state_counter();
66
67 if (enthalpy_revision == m_ice_enthalpy.state_counter()) {
68 bootstrap_ice_enthalpy(ice_thickness, surface_temperature, climatic_mass_balance,
69 basal_heat_flux, m_ice_enthalpy);
70 }
71}
72
74 const array::Scalar &ice_thickness,
75 const array::Scalar &surface_temperature,
76 const array::Scalar &climatic_mass_balance,
77 const array::Scalar &basal_heat_flux) {
78
79 m_log->message(2, "* Bootstrapping the enthalpy-based energy balance model...\n");
80
82
83 regrid("Energy balance model", m_basal_melt_rate, REGRID_WITHOUT_REGRID_VARS);
84
85 int enthalpy_revision = m_ice_enthalpy.state_counter();
87
88 if (enthalpy_revision == m_ice_enthalpy.state_counter()) {
89 bootstrap_ice_enthalpy(ice_thickness, surface_temperature, climatic_mass_balance,
90 basal_heat_flux, m_ice_enthalpy);
91 }
92}
93
94//! Update ice enthalpy field based on conservation of energy.
95/*!
96This method is documented by the page \ref bombproofenth and by [\ref
97AschwandenBuelerKhroulevBlatter].
98
99This method updates array::Array3D m_work and array::Scalar basal_melt_rate.
100No communication of ghosts is done for any of these fields.
101
102We use an instance of enthSystemCtx.
103
104Regarding drainage, see [\ref AschwandenBuelerKhroulevBlatter] and references therein.
105 */
106
107void EnthalpyModel::update_impl(double t, double dt, const Inputs &inputs) {
108 // current time does not matter here
109 (void) t;
110
111 EnthalpyConverter::Ptr EC = m_grid->ctx()->enthalpy_converter();
112
113 const double
114 ice_density = m_config->get_number("constants.ice.density"), // kg m-3
115 bulgeEnthMax = m_config->get_number("energy.enthalpy.cold_bulge_max"), // J kg-1
116 target_water_fraction = m_config->get_number("energy.drainage_target_water_fraction");
117
119
120 inputs.check();
121
122 // give them names that are a bit shorter...
123 const array::Array3D
124 &strain_heating3 = *inputs.volumetric_heating_rate,
125 &u3 = *inputs.u3,
126 &v3 = *inputs.v3,
127 &w3 = *inputs.w3;
128
129 const auto &cell_type = *inputs.cell_type;
130
131 const array::Scalar
132 &basal_frictional_heating = *inputs.basal_frictional_heating,
133 &basal_heat_flux = *inputs.basal_heat_flux,
134 &surface_liquid_fraction = *inputs.surface_liquid_fraction,
135 &shelf_base_temp = *inputs.shelf_base_temp,
136 &ice_surface_temp = *inputs.surface_temp,
137 &till_water_thickness = *inputs.till_water_thickness;
138
139 const array::Scalar1 &ice_thickness = *inputs.ice_thickness;
140
141 energy::enthSystemCtx system(m_grid->z(), "energy.enthalpy", m_grid->dx(), m_grid->dy(), dt,
142 *m_config, m_ice_enthalpy, u3, v3, w3, strain_heating3, EC);
143
144 const size_t Mz_fine = system.z().size();
145 const double dz = system.dz();
146 std::vector<double> Enthnew(Mz_fine); // new enthalpy in column
147
148 array::AccessScope list{&ice_surface_temp, &shelf_base_temp, &surface_liquid_fraction,
149 &ice_thickness, &basal_frictional_heating, &basal_heat_flux, &till_water_thickness,
150 &cell_type, &u3, &v3, &w3, &strain_heating3, &m_basal_melt_rate, &m_ice_enthalpy,
151 &m_work};
152
153 double margin_threshold = m_config->get_number("energy.margin_ice_thickness_limit");
154
155 unsigned int liquifiedCount = 0;
156
157 ParallelSection loop(m_grid->com);
158 try {
159 for (auto pt = m_grid->points(); pt; pt.next()) {
160 const int i = pt.i(), j = pt.j();
161
162 const double H = ice_thickness(i, j);
163
164 system.init(i, j,
165 marginal(ice_thickness, i, j, margin_threshold),
166 H);
167
168 // enthalpy and pressures at top of ice
169 const double
170 depth_ks = H - system.ks() * dz,
171 p_ks = EC->pressure(depth_ks); // FIXME issue #15
172
173 const double Enth_ks = EC->enthalpy_permissive(ice_surface_temp(i, j),
174 surface_liquid_fraction(i, j), p_ks);
175
176 const bool ice_free_column = (system.ks() == 0);
177
178 // deal completely with columns with no ice; enthalpy and basal_melt_rate need setting
179 if (ice_free_column) {
180 m_work.set_column(i, j, Enth_ks);
181 // The floating basal melt rate will be set later; cover this
182 // case and set to zero for now. Also, there is no basal melt
183 // rate on ice free land and ice free ocean
184 m_basal_melt_rate(i, j) = 0.0;
185 continue;
186 } // end of if (ice_free_column)
187
188 if (system.lambda() < 1.0) {
189 m_stats.reduced_accuracy_counter += 1; // count columns with lambda < 1
190 }
191
192 const bool
193 is_floating = cell_type.ocean(i, j),
194 base_is_warm = system.Enth(0) >= system.Enth_s(0),
195 above_base_is_warm = system.Enth(1) >= system.Enth_s(1);
196
197 // set boundary conditions and update enthalpy
198 {
199 system.set_surface_dirichlet_bc(Enth_ks);
200
201 // determine lowest-level equation at bottom of ice; see
202 // decision chart in the source code browser and page
203 // documenting BOMBPROOF
204 if (is_floating) {
205 // floating base: Dirichlet application of known temperature from ocean
206 // coupler; assumes base of ice shelf has zero liquid fraction
207 double Enth0 = EC->enthalpy_permissive(shelf_base_temp(i, j), 0.0, EC->pressure(H));
208
209 system.set_basal_dirichlet_bc(Enth0);
210 } else {
211 // grounded ice warm and wet
212 if (base_is_warm && (till_water_thickness(i, j) > 0.0)) {
213 if (above_base_is_warm) {
214 // temperate layer at base (Neumann) case: q . n = 0 (K0 grad E . n = 0)
215 system.set_basal_heat_flux(0.0);
216 } else {
217 // only the base is warm: E = E_s(p) (Dirichlet)
218 // ( Assumes ice has zero liquid fraction. Is this a valid assumption here?
219 system.set_basal_dirichlet_bc(system.Enth_s(0));
220 }
221 } else {
222 // (Neumann) case: q . n = q_lith . n + F_b
223 // a) cold and dry base, or
224 // b) base that is still warm from the last time step, but without basal water
225 system.set_basal_heat_flux(basal_heat_flux(i, j) + basal_frictional_heating(i, j));
226 }
227 }
228
229 // solve the system
230 system.solve(Enthnew);
231
232 }
233
234 // post-process (drainage and bulge-limiting)
235 double Hdrainedtotal = 0.0;
236 {
237 // drain ice segments by mechanism in [\ref AschwandenBuelerKhroulevBlatter],
238 // using DrainageCalculator dc
239 for (unsigned int k=0; k < system.ks(); k++) {
240 if (Enthnew[k] > system.Enth_s(k)) { // avoid doing any more work if cold
241
242 const double
243 depth = H - k * dz,
244 p = EC->pressure(depth), // FIXME issue #15
245 T_m = EC->melting_temperature(p),
246 L = EC->L(T_m);
247
248 if (Enthnew[k] >= system.Enth_s(k) + 0.5 * L) {
249 liquifiedCount++; // count these rare events...
250 Enthnew[k] = system.Enth_s(k) + 0.5 * L; // but lose the energy
251 }
252
253 double omega = EC->water_fraction(Enthnew[k], p);
254
255 if (omega > target_water_fraction) {
256 double fractiondrained = dc.get_drainage_rate(omega) * dt; // pure number
257
258 fractiondrained = std::min(fractiondrained,
259 omega - target_water_fraction);
260 Hdrainedtotal += fractiondrained * dz; // always a positive contribution
261 Enthnew[k] -= fractiondrained * L;
262 }
263 }
264 }
265
266 // apply bulge limiter
267 const double lowerEnthLimit = Enth_ks - bulgeEnthMax;
268 for (unsigned int k=0; k < system.ks(); k++) {
269 if (Enthnew[k] < lowerEnthLimit) {
270 // Count grid points which have very large cold limit advection bulge... enthalpy not
271 // too low.
273 Enthnew[k] = lowerEnthLimit;
274 }
275 }
276
277 // if there is subglacial water, don't allow ice base enthalpy to be below
278 // pressure-melting; that is, assume subglacial water is at the pressure-
279 // melting temperature and enforce continuity of temperature
280 if (till_water_thickness(i, j) > 0.0) {
281 Enthnew[0] = std::max(Enthnew[0], system.Enth_s(0));
282 }
283 } // end of post-processing
284
285 // compute basal melt rate
286 {
287 bool base_is_cold = (Enthnew[0] < system.Enth_s(0)) && (till_water_thickness(i,j) == 0.0);
288 // Determine melt rate, but only preliminarily because of
289 // drainage, from heat flux out of bedrock, heat flux into
290 // ice, and frictional heating
291 if (is_floating) {
292 // The floating basal melt rate will be set later; cover
293 // this case and set to zero for now. Note that
294 // Hdrainedtotal is discarded (the ocean model determines
295 // the basal melt).
296 m_basal_melt_rate(i, j) = 0.0;
297 } else {
298 if (base_is_cold) {
299 m_basal_melt_rate(i, j) = 0.0; // zero melt rate if cold base
300 } else {
301 const double
302 p_0 = EC->pressure(H),
303 p_1 = EC->pressure(H - dz), // FIXME issue #15
304 Tpmp_0 = EC->melting_temperature(p_0);
305
306 const bool k1_istemperate = EC->is_temperate(Enthnew[1], p_1); // level z = + \Delta z
307 double hf_up = 0.0;
308 if (k1_istemperate) {
309 const double
310 Tpmp_1 = EC->melting_temperature(p_1);
311
312 hf_up = -system.k_from_T(Tpmp_0) * (Tpmp_1 - Tpmp_0) / dz;
313 } else {
314 double T_0 = EC->temperature(Enthnew[0], p_0);
315 const double K_0 = system.k_from_T(T_0) / EC->c();
316
317 hf_up = -K_0 * (Enthnew[1] - Enthnew[0]) / dz;
318 }
319
320 // compute basal melt rate from flux balance:
321 //
322 // basal_melt_rate = - Mb / rho in [\ref AschwandenBuelerKhroulevBlatter];
323 //
324 // after we compute it we make sure there is no refreeze if
325 // there is no available basal water
326 m_basal_melt_rate(i, j) = (basal_frictional_heating(i, j) + basal_heat_flux(i, j) - hf_up) / (ice_density * EC->L(Tpmp_0));
327
328 if (till_water_thickness(i, j) <= 0 && m_basal_melt_rate(i, j) < 0) {
329 m_basal_melt_rate(i, j) = 0.0;
330 }
331 }
332
333 // Add drained water from the column to basal melt rate.
334 m_basal_melt_rate(i, j) += Hdrainedtotal / dt;
335 } // end of the grounded case
336 } // end of the basal melt rate computation
337
338 system.fine_to_coarse(Enthnew, i, j, m_work);
339 }
340 } catch (...) {
341 loop.failed();
342 }
343 loop.check();
344
345 m_stats.liquified_ice_volume = ((double) liquifiedCount) * dz * m_grid->cell_area();
346}
347
352
354 m_ice_enthalpy.write(output);
355 m_basal_melt_rate.write(output);
356}
357
358} // end of namespace energy
359} // end of namespace pism
const Config::ConstPtr m_config
configuration database used by this component
Definition Component.hh:158
const Logger::ConstPtr m_log
logger (for easy access)
Definition Component.hh:162
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition Component.hh:156
void regrid(const std::string &module_name, array::Array &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS)
Definition Component.cc:159
std::shared_ptr< EnthalpyConverter > Ptr
std::string name() const
Definition File.cc:274
High-level PISM I/O class.
Definition File.hh:55
void failed()
Indicates a failure of a parallel section.
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition Array.hh:64
void copy_from(const Array2D< T > &source)
Definition Array2D.hh:73
void set_column(int i, int j, double c)
Set all values of scalar quantity to given a single value in a particular column.
Definition Array3D.cc:50
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition Array3D.hh:33
void read(const std::string &filename, unsigned int time)
Definition Array.cc:731
void define(const File &file, io::Type default_type) const
Define variables corresponding to an Array in a file opened using file.
Definition Array.cc:463
void write(const std::string &filename) const
Definition Array.cc:722
int state_counter() const
Get the object state counter.
Definition Array.cc:127
void regrid(const std::string &filename, io::Default default_value)
Definition Array.cc:736
const std::vector< double > & z() const
void fine_to_coarse(const std::vector< double > &input, int i, int j, array::Array3D &output) const
unsigned int ks() const
virtual double get_drainage_rate(double omega)
Return D(omega), as in figure in [AschwandenBuelerKhroulevBlatter].
Compute the rate of drainage D(omega) for temperate ice.
void init_enthalpy(const File &input_file, bool regrid, int record)
Initialize enthalpy by reading it from a file, or by reading temperature and liquid water fraction,...
EnergyModelStats m_stats
const array::Scalar & basal_melt_rate() const
Basal melt rate in grounded areas. (It is set to zero elsewhere.)
array::Array3D m_ice_enthalpy
array::Scalar m_basal_melt_rate
void regrid_enthalpy()
Regrid enthalpy from the -regrid_file.
virtual void initialize_impl(const array::Scalar &basal_melt_rate, const array::Scalar &ice_thickness, const array::Scalar &surface_temperature, const array::Scalar &climatic_mass_balance, const array::Scalar &basal_heat_flux)
virtual void restart_impl(const File &input_file, int record)
virtual void bootstrap_impl(const File &input_file, const array::Scalar &ice_thickness, const array::Scalar &surface_temperature, const array::Scalar &climatic_mass_balance, const array::Scalar &basal_heat_flux)
virtual void update_impl(double t, double dt, const Inputs &inputs)
Update ice enthalpy field based on conservation of energy.
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
EnthalpyModel(std::shared_ptr< const Grid > grid, std::shared_ptr< const stressbalance::StressBalance > stress_balance)
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
const array::Scalar * surface_liquid_fraction
const array::Scalar1 * ice_thickness
const array::Array3D * w3
const array::Array3D * v3
const array::Scalar * shelf_base_temp
const array::Scalar * basal_heat_flux
const array::Scalar * basal_frictional_heating
const array::Scalar * till_water_thickness
const array::Array3D * u3
const array::Scalar * surface_temp
const array::CellType * cell_type
const array::Array3D * volumetric_heating_rate
double Enth(size_t i) const
Definition enthSystem.hh:73
double Enth_s(size_t i) const
Definition enthSystem.hh:77
void init(int i, int j, bool ismarginal, double ice_thickness)
double k_from_T(double T) const
void set_basal_heat_flux(double heat_flux)
Set coefficients in discrete equation for Neumann condition at base of ice.
void solve(std::vector< double > &result)
Solve the tridiagonal system, in a single column, which determines the new values of the ice enthalpy...
void set_basal_dirichlet_bc(double E_basal)
Set coefficients in discrete equation for at base of ice.
void set_surface_dirichlet_bc(double E_surface)
Tridiagonal linear system for conservation of energy in vertical column of ice enthalpy.
Definition enthSystem.hh:38
static const double L
Definition exactTestL.cc:40
bool marginal(const array::Scalar1 &thickness, int i, int j, double threshold)
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
@ PISM_DOUBLE
Definition IO_Flags.hh:52
static const double k
Definition exactTestP.cc:42