tlocalMassBalance.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tlocalMassBalance.cc (17051B)
---
1 // Copyright (C) 2009, 2010, 2011, 2013, 2014, 2015, 2016, 2017, 2018 Ed Bueler and Constantine Khroulev and Andy Aschwanden
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 #include <cassert>
20 #include <ctime> // for time(), used to initialize random number gen
21 #include <gsl/gsl_rng.h>
22 #include <gsl/gsl_randist.h>
23 #include <gsl/gsl_math.h> // M_PI
24 #include <cmath> // for erfc() in CalovGreveIntegrand()
25 #include <algorithm>
26
27 #include "pism/util/pism_utilities.hh"
28 #include "pism/util/ConfigInterface.hh"
29 #include "localMassBalance.hh"
30 #include "pism/util/IceGrid.hh"
31
32 namespace pism {
33 namespace surface {
34
35 LocalMassBalance::Changes::Changes() {
36 firn_depth = 0.0;
37 snow_depth = 0.0;
38 melt = 0.0;
39 runoff = 0.0;
40 smb = 0.0;
41 }
42
43 LocalMassBalance::LocalMassBalance(Config::ConstPtr myconfig, units::System::Ptr system)
44 : m_config(myconfig), m_unit_system(system),
45 m_seconds_per_day(86400) {
46 // empty
47 }
48
49 LocalMassBalance::~LocalMassBalance() {
50 // empty
51 }
52
53 std::string LocalMassBalance::method() const {
54 return m_method;
55 }
56
57 PDDMassBalance::PDDMassBalance(Config::ConstPtr config, units::System::Ptr system)
58 : LocalMassBalance(config, system) {
59 precip_as_snow = m_config->get_flag("surface.pdd.interpret_precip_as_snow");
60 Tmin = m_config->get_number("surface.pdd.air_temp_all_precip_as_snow");
61 Tmax = m_config->get_number("surface.pdd.air_temp_all_precip_as_rain");
62 pdd_threshold_temp = m_config->get_number("surface.pdd.positive_threshold_temp");
63 refreeze_ice_melt = m_config->get_flag("surface.pdd.refreeze_ice_melt");
64
65 m_method = "an expectation integral";
66 }
67
68
69 /*! \brief Compute the number of points for temperature and
70 precipitation time-series.
71 */
72 unsigned int PDDMassBalance::get_timeseries_length(double dt) {
73 const unsigned int NperYear = static_cast<unsigned int>(m_config->get_number("surface.pdd.max_evals_per_year"));
74 const double dt_years = units::convert(m_unit_system, dt, "seconds", "years");
75
76 return std::max(1U, static_cast<unsigned int>(ceil(NperYear * dt_years)));
77 }
78
79
80 //! Compute the integrand in integral (6) in [\ref CalovGreve05].
81 /*!
82 The integral is
83 \f[\mathrm{PDD} = \int_{t_0}^{t_0+\mathtt{dt}} dt\,
84 \bigg[\frac{\sigma}{\sqrt{2\pi}}\,\exp\left(-\frac{T_{ac}(t)^2}{2\sigma^2}\right)
85 + \frac{T_{ac}(t)}{2}\,\mathrm{erfc}
86 \left(-\frac{T_{ac}(t)}{\sqrt{2}\,\sigma}\right)\bigg] \f]
87 This procedure computes the quantity in square brackets. The value \f$T_{ac}(t)\f$
88 in the above integral is in degrees C. Here we think of the argument `TacC`
89 as temperature in Celsius, but really it is the temperature above a threshold
90 at which it is "positive".
91
92 This integral is used for the expected number of positive degree days. The user can choose
93 \f$\sigma\f$ by option `-pdd_std_dev`. Note that the integral is over a time interval of
94 length `dt` instead of a whole year as stated in \ref CalovGreve05 . If `sigma` is zero,
95 return the positive part of `TacC`.
96 */
97 double PDDMassBalance::CalovGreveIntegrand(double sigma, double TacC) {
98
99 if (sigma == 0) {
100 return std::max(TacC, 0.0);
101 } else {
102 const double Z = TacC / (sqrt(2.0) * sigma);
103 return (sigma / sqrt(2.0 * M_PI)) * exp(-Z*Z) + (TacC / 2.0) * erfc(-Z);
104 }
105 }
106
107
108 //! Compute the expected number of positive degree days from the input temperature time-series.
109 /**
110 * Use the rectangle method for simplicity.
111 *
112 * @param S standard deviation for air temperature excursions
113 * @param dt_series length of the step for the time-series
114 * @param T air temperature (array of length N)
115 * @param N length of the T array
116 * @param[out] PDDs pointer to a pre-allocated array with N-1 elements
117 */
118 void PDDMassBalance::get_PDDs(double dt_series,
119 const std::vector<double> &S,
120 const std::vector<double> &T,
121 std::vector<double> &PDDs) {
122 assert(S.size() == T.size() and T.size() == PDDs.size());
123 assert(dt_series > 0.0);
124
125 const double h_days = dt_series / m_seconds_per_day;
126 const size_t N = S.size();
127
128 for (unsigned int k = 0; k < N; ++k) {
129 PDDs[k] = h_days * CalovGreveIntegrand(S[k], T[k] - pdd_threshold_temp);
130 }
131 }
132
133
134 //! \brief Extract snow accumulation from mixed (snow and rain)
135 //! precipitation using the temperature time-series.
136 /** Uses the temperature time-series to determine whether the
137 * precipitation is snow or rain. Rain is removed entirely from the
138 * surface mass balance, and will not be included in the computed
139 * runoff, which is meltwater runoff. There is an allowed linear
140 * transition for Tmin below which all precipitation is interpreted as
141 * snow, and Tmax above which all precipitation is rain (see, e.g.
142 * [\ref Hock2005b]).
143 *
144 * Sets P[i] to the *solid* (snow) accumulation *rate*.
145 *
146 * @param[in,out] P precipitation rate (array of length N)
147 * @param[in] T air temperature (array of length N)
148 * @param[in] N array length
149 */
150 void PDDMassBalance::get_snow_accumulation(const std::vector<double> &T,
151 std::vector<double> &P) {
152
153 assert(T.size() == P.size());
154 const size_t N = T.size();
155
156 // Following \ref Hock2005b we employ a linear transition from Tmin to Tmax
157 for (unsigned int i = 0; i < N; i++) {
158 // do not allow negative precipitation
159 if (P[i] < 0.0) {
160 P[i] = 0.0;
161 continue;
162 }
163
164 if (precip_as_snow || T[i] <= Tmin) { // T <= Tmin, all precip is snow
165 // no change
166 } else if (T[i] < Tmax) { // linear transition from Tmin to Tmax
167 P[i] *= (Tmax - T[i]) / (Tmax - Tmin);
168 } else { // T >= Tmax, all precip is rain -- ignore it
169 P[i] = 0.0;
170 }
171 }
172
173 }
174
175
176 //! \brief Compute the surface mass balance at a location from the number of positive
177 //! degree days and the accumulation amount in a time interval.
178 /*!
179 * This is a PDD scheme. The input parameter `ddf.snow` is a rate of
180 * melting per positive degree day for snow.
181 *
182 * `accumulation` has units "meter / second".
183 *
184 * - a fraction of the melted snow and ice refreezes, conceptualized
185 * as superimposed ice, and this is controlled by parameter \c
186 * ddf.refreeze_fraction
187 *
188 * - the excess number of PDDs is used to melt both the ice that came
189 * from refreeze and then any ice which is already present.
190 *
191 * Ice melts at a constant rate per positive degree day, controlled by
192 * parameter `ddf.ice`.
193 *
194 * The scheme here came from EISMINT-Greenland [\ref RitzEISMINT], but
195 * is influenced by R. Hock (personal communication).
196 */
197 PDDMassBalance::Changes PDDMassBalance::step(const DegreeDayFactors &ddf,
198 double PDDs,
199 double thickness,
200 double old_firn_depth,
201 double old_snow_depth,
202 double accumulation) {
203 double
204 firn_depth = old_firn_depth,
205 snow_depth = old_snow_depth,
206 max_snow_melted = PDDs * ddf.snow,
207 firn_melted = 0.0,
208 snow_melted = 0.0,
209 excess_pdds = 0.0;
210
211 assert(thickness >= 0);
212
213 // snow depth cannot exceed total thickness
214 snow_depth = std::min(snow_depth, thickness);
215
216 assert(snow_depth >= 0);
217
218 // firn depth cannot exceed thickness - snow_depth
219 firn_depth = std::min(firn_depth, thickness - snow_depth);
220
221 assert(firn_depth >= 0);
222
223 double ice_thickness = thickness - snow_depth - firn_depth;
224
225 assert(ice_thickness >= 0);
226
227 snow_depth += accumulation;
228
229 if (PDDs <= 0.0) { // The "no melt" case.
230 snow_melted = 0.0;
231 firn_melted = 0.0,
232 excess_pdds = 0.0;
233 } else if (max_snow_melted <= snow_depth) {
234 // Some of the snow melted and some is left; in any case, all of
235 // the energy available for melt, namely all of the positive
236 // degree days (PDDs) were used up in melting snow.
237 snow_melted = max_snow_melted;
238 firn_melted = 0.0;
239 excess_pdds = 0.0;
240 } else if (max_snow_melted <= firn_depth + snow_depth) {
241 // All of the snow is melted but some firn is left; in any case, all of
242 // the energy available for melt, namely all of the positive
243 // degree days (PDDs) were used up in melting snow and firn.
244 snow_melted = snow_depth;
245 firn_melted = max_snow_melted - snow_melted;
246 excess_pdds = 0.0;
247 } else {
248 // All (firn_depth and snow_depth meters) of snow and firn melted. Excess_pdds is the
249 // positive degree days available to melt ice.
250 firn_melted = firn_depth;
251 snow_melted = snow_depth;
252 excess_pdds = PDDs - ((firn_melted + snow_melted) / ddf.snow); // units: K day
253 }
254
255 double
256 ice_melted = std::min(excess_pdds * ddf.ice, ice_thickness),
257 melt = snow_melted + firn_melted + ice_melted,
258 ice_created_by_refreeze = 0.0;
259
260 if (refreeze_ice_melt) {
261 ice_created_by_refreeze = melt * ddf.refreeze_fraction;
262 } else {
263 // Should this only be snow melted?
264 ice_created_by_refreeze = (firn_melted + snow_melted) * ddf.refreeze_fraction;
265 }
266
267 const double runoff = melt - ice_created_by_refreeze;
268
269 snow_depth = std::max(snow_depth - snow_melted, 0.0);
270 firn_depth = std::max(firn_depth - firn_melted, 0.0);
271
272 // FIXME: need to add snow that hasn't melted, is this correct?
273 // firn_depth += (snow_depth - snow_melted);
274 // Turn firn into ice at X times accumulation
275 // firn_depth -= accumulation * m_config->get_number("surface.pdd.firn_compaction_to_accumulation_ratio");
276
277 const double smb = accumulation - runoff;
278
279 Changes result;
280 // Ensure that we never generate negative ice thicknesses. As far as I can tell the code
281 // above guarantees that thickness + smb >= *in exact arithmetic*. The check below
282 // should make sure that we don't get bitten by rounding errors.
283 result.smb = thickness + smb >= 0 ? smb : -thickness;
284 result.firn_depth = firn_depth - old_firn_depth;
285 result.snow_depth = snow_depth - old_snow_depth;
286 result.melt = melt;
287 result.runoff = runoff;
288
289 assert(thickness + result.smb >= 0);
290
291 return result;
292 }
293
294
295 /*!
296 Initializes the random number generator (RNG). The RNG is GSL's recommended default,
297 which seems to be "mt19937" and is DIEHARD (whatever that means ...). Seed with
298 wall clock time in seconds in non-repeatable case, and with 0 in repeatable case.
299 */
300 PDDrandMassBalance::PDDrandMassBalance(Config::ConstPtr config, units::System::Ptr system,
301 Kind kind)
302 : PDDMassBalance(config, system) {
303 pddRandGen = gsl_rng_alloc(gsl_rng_default); // so pddRandGen != NULL now
304 gsl_rng_set(pddRandGen, kind == REPEATABLE ? 0 : time(0));
305
306 m_method = (kind == NOT_REPEATABLE
307 ? "simulation of a random process"
308 : "repeatable simulation of a random process");
309 }
310
311
312 PDDrandMassBalance::~PDDrandMassBalance() {
313 if (pddRandGen != NULL) {
314 gsl_rng_free(pddRandGen);
315 pddRandGen = NULL;
316 }
317 }
318
319
320 /*! We need to compute simulated random temperature each actual \e
321 day, or at least as close as we can reasonably get. Output `N` is
322 number of days or number of days plus one.
323
324 Thus this method ignores
325 `config.get_number("surface.pdd.max_evals_per_year")`, which is
326 used in the base class PDDMassBalance.
327
328 Implementation of get_PDDs() requires returned N >= 2, so we
329 guarantee that.
330 */
331 unsigned int PDDrandMassBalance::get_timeseries_length(double dt) {
332 return std::max(static_cast<size_t>(ceil(dt / m_seconds_per_day)), (size_t)2);
333 }
334
335 /**
336 * Computes
337 * \f[
338 * \text{PDD} = \sum_{i=0}^{N-1} h_{\text{days}} \cdot \text{max}(T_i-T_{\text{threshold}}, 0).
339 * \f]
340 *
341 * @param S \f$\sigma\f$ (standard deviation for daily temperature excursions)
342 * @param dt_series time-series step, in seconds
343 * @param T air temperature
344 * @param N number of points in the temperature time-series, each corresponds to a sub-interval
345 * @param PDDs pointer to a pre-allocated array of length N
346 */
347 void PDDrandMassBalance::get_PDDs(double dt_series,
348 const std::vector<double> &S,
349 const std::vector<double> &T,
350 std::vector<double> &PDDs) {
351 assert(S.size() == T.size() and T.size() == PDDs.size());
352 assert(dt_series > 0.0);
353
354 const double h_days = dt_series / m_seconds_per_day;
355 const size_t N = S.size();
356
357 for (unsigned int k = 0; k < N; ++k) {
358 // average temperature in k-th interval
359 double T_k = T[k] + gsl_ran_gaussian(pddRandGen, S[k]); // add random: N(0,sigma)
360
361 if (T_k > pdd_threshold_temp) {
362 PDDs[k] = h_days * (T_k - pdd_threshold_temp);
363 }
364 }
365 }
366
367
368 FaustoGrevePDDObject::FaustoGrevePDDObject(IceGrid::ConstPtr g)
369 : m_grid(g), m_config(g->ctx()->config()) {
370
371 m_beta_ice_w = m_config->get_number("surface.pdd.fausto.beta_ice_w");
372 m_beta_snow_w = m_config->get_number("surface.pdd.fausto.beta_snow_w");
373
374 m_T_c = m_config->get_number("surface.pdd.fausto.T_c");
375 m_T_w = m_config->get_number("surface.pdd.fausto.T_w");
376 m_beta_ice_c = m_config->get_number("surface.pdd.fausto.beta_ice_c");
377 m_beta_snow_c = m_config->get_number("surface.pdd.fausto.beta_snow_c");
378
379 m_fresh_water_density = m_config->get_number("constants.fresh_water.density");
380 m_ice_density = m_config->get_number("constants.ice.density");
381 m_pdd_fausto_latitude_beta_w = m_config->get_number("surface.pdd.fausto.latitude_beta_w");
382 m_refreeze_fraction = m_config->get_number("surface.pdd.refreeze");
383
384
385 m_temp_mj.create(m_grid, "temp_mj_faustogreve", WITHOUT_GHOSTS);
386 m_temp_mj.set_attrs("internal",
387 "mean July air temp from Fausto et al (2009) parameterization",
388 "K", "K", "", 0);
389 }
390
391 FaustoGrevePDDObject::~FaustoGrevePDDObject() {
392 // empty
393 }
394
395 LocalMassBalance::DegreeDayFactors FaustoGrevePDDObject::degree_day_factors(int i, int j,
396 double latitude) {
397
398 LocalMassBalance::DegreeDayFactors ddf;
399 ddf.refreeze_fraction = m_refreeze_fraction;
400
401 IceModelVec::AccessList list(m_temp_mj);
402 const double T_mj = m_temp_mj(i,j);
403
404 if (latitude < m_pdd_fausto_latitude_beta_w) { // case latitude < 72 deg N
405 ddf.ice = m_beta_ice_w;
406 ddf.snow = m_beta_snow_w;
407 } else { // case > 72 deg N
408 if (T_mj >= m_T_w) {
409 ddf.ice = m_beta_ice_w;
410 ddf.snow = m_beta_snow_w;
411 } else if (T_mj <= m_T_c) {
412 ddf.ice = m_beta_ice_c;
413 ddf.snow = m_beta_snow_c;
414 } else { // middle case T_c < T_mj < T_w
415 const double
416 lam_i = pow((m_T_w - T_mj) / (m_T_w - m_T_c) , 3.0),
417 lam_s = (T_mj - m_T_c) / (m_T_w - m_T_c);
418 ddf.ice = m_beta_ice_w + (m_beta_ice_c - m_beta_ice_w) * lam_i;
419 ddf.snow = m_beta_snow_w + (m_beta_snow_c - m_beta_snow_w) * lam_s;
420 }
421 }
422
423 // degree-day factors in \ref Faustoetal2009 are water-equivalent
424 // thickness per degree day; ice-equivalent thickness melted per degree
425 // day is slightly larger; for example, iwfactor = 1000/910
426 const double iwfactor = m_fresh_water_density / m_ice_density;
427 ddf.snow *= iwfactor;
428 ddf.ice *= iwfactor;
429
430 return ddf;
431 }
432
433
434 //! Updates mean July near-surface air temperature.
435 /*!
436 Unfortunately this duplicates code in SeaRISEGreenland::update();
437 */
438 void FaustoGrevePDDObject::update_temp_mj(const IceModelVec2S &surfelev,
439 const IceModelVec2S &lat,
440 const IceModelVec2S &lon) {
441 const double
442 d_mj = m_config->get_number("atmosphere.fausto_air_temp.d_mj"), // K
443 gamma_mj = m_config->get_number("atmosphere.fausto_air_temp.gamma_mj"), // K m-1
444 c_mj = m_config->get_number("atmosphere.fausto_air_temp.c_mj"), // K (degN)-1
445 kappa_mj = m_config->get_number("atmosphere.fausto_air_temp.kappa_mj"); // K (degW)-1
446
447 const IceModelVec2S
448 &h = surfelev,
449 &lat_degN = lat,
450 &lon_degE = lon;
451
452 IceModelVec::AccessList list{&h, &lat_degN, &lon_degE, &m_temp_mj};
453
454 for (Points p(*m_grid); p; p.next()) {
455 const int i = p.i(), j = p.j();
456 m_temp_mj(i,j) = d_mj + gamma_mj * h(i,j) + c_mj * lat_degN(i,j) + kappa_mj * (-lon_degE(i,j));
457 }
458 }
459
460 } // end of namespace surface
461 } // end of namespace pism