Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/Glaciers.cc

1/***********************************************************************************/
2/* Copyright 2009-2015 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
3/***********************************************************************************/
4/* This file is part of Alpine3D.
5 Alpine3D is free software: you can redistribute it and/or modify
6 it under the terms of the GNU Lesser General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 Alpine3D is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public License
16 along with Alpine3D. If not, see <http://www.gnu.org/licenses/>.
17*/
18#include <alpine3d/Glaciers.h>
19#include <alpine3d/MPIControl.h>
20
21using namespace std;
22using namespace mio;
23
24Glaciers::Glaciers(const mio::Config& cfg)
25 : dem(), isGlacier(), flowpath(), src_altitude(), src_distance(), KBL(IOUtils::nodata), K(0.), scale(1.)
26{
27init(cfg);
28}
29
30Glaciers::Glaciers(const mio::Config& cfg, const mio::DEMObject& in_dem)
31 : dem(in_dem), isGlacier(), flowpath(), src_altitude(), src_distance(), KBL(IOUtils::nodata), K(0.), scale(1.)
32{
33init(cfg);
34}
35
36void Glaciers::init(const mio::Config& cfg)
37{
38cfg.getValue("KATABATIC_LAYER_HEIGHT", "Snowpack", KBL, IOUtils::dothrow); //Katabatic Boundary Layer (KBL) thickness, 17m in GB
39cfg.getValue("KATABATIC_SCALING", "Snowpack", scale, IOUtils::nothrow); //a value of .25 seems to work quite well
40cfg.getValue("KATABATIC_K_COEFFICIENT", "Snowpack", K, IOUtils::nothrow); //K=0 for original GB model
41
42if (MPIControl::instance().master())
43std::cout << "[i] Enabled katabatic flows TA corrections\n";
44}
45
46void Glaciers::setDEM(const mio::DEMObject& in_dem)
47{
48dem = in_dem;
49if (!isGlacier.empty()) {
50//apply dem mask
51for (size_t ii=0; ii<isGlacier.size(); ii++) {
52if (dem(ii)==IOUtils::nodata)
53isGlacier(ii) = IOUtils::nodata;
54}
55hillslope_flow(isGlacier); //compute the flow paths
56}
57}
58
59/**
60 * @brief Set the mask of where the glaciers are.
61 * @param[in] glacierMask grid containing IOUtils::nodata at glaciated pixels
62 */
63void Glaciers::setGlacierMap(const Grid2DObject& glacierMask)
64{
65//convert the mask that remove the glaciers to a mask that keeps the glaciers
66isGlacier.set(glacierMask, 0.); //non glacier pixels tagged with 0
67for (size_t idx=0; idx<glacierMask.size(); idx++) { //tag glacier pixels with 1
68if (glacierMask(idx)==IOUtils::nodata)
69isGlacier(idx) = 1.;
70}
71
72if (!dem.empty()) {
73//apply dem mask
74for (size_t ii=0; ii<isGlacier.size(); ii++) {
75if (dem(ii)==IOUtils::nodata)
76isGlacier(ii) = IOUtils::nodata;
77}
78hillslope_flow(isGlacier); //compute the flow paths
79}
80}
81
82void Glaciers::getGrids(Grid2DObject &alt, Grid2DObject &dist) const
83{
84alt = src_altitude;
85dist = src_distance;
86}
87
88/**
89 * @brief Based on the air and surface temperatures as well as the fraction of the domain that is snow covered,
90 * enable or disable the katabatic flows (the fractions are currently set at 20%).
91 * @param[in] hs snow height grid
92 * @param[in] tss surface temperature grid
93 * @param[in] ta air temperature grid
94 * @param[in] isGlacier map of glaciated pixels
95 * @return should katabatic flows be computed?
96 */
97bool Glaciers::enableKatabatikFlows(const mio::Grid2DObject& hs, const mio::Grid2DObject& tss, const mio::Grid2DObject& ta, const mio::Grid2DObject& isGlacier)
98{
99if (!ta.isSameGeolocalization(hs))
100throw IOException("The temperature grid and the provided snow height grids don't have the same geolocalization!", AT);
101if (!ta.isSameGeolocalization(tss))
102throw IOException("The temperature grid and the provided dem don't have the same geolocalization!", AT);
103
104size_t glacier_pixels = 0;
105size_t non_glacier_pixels = 0;
106size_t n_snowfree = 0; //count cells that should experience katabatic winds
107size_t n_glacier_katabatic = 0; //count glacier cells that should experience katabatic winds
108
109for (size_t idx=0; idx<isGlacier.size(); idx++) {
110if (isGlacier(idx)==IOUtils::nodata) continue; //not in the simulation
111
112if (isGlacier(idx)==0) { //out of the glacier
113non_glacier_pixels++;
114if (hs(idx)<=0.) n_snowfree++;
115} else { //on the glacier
116glacier_pixels++;
117if (ta(idx)>tss(idx)) n_glacier_katabatic++;
118}
119}
120
121if (n_snowfree>(non_glacier_pixels/5) && n_glacier_katabatic>(glacier_pixels/5))
122return true;
123
124return false;
125}
126
127const mio::Grid2DObject Glaciers::correctTemperatures(const mio::Grid2DObject& hs, const mio::Grid2DObject& tss, const mio::Grid2DObject& ta) const
128{
129Grid2DObject new_ta(ta);
130correctTemperatures(hs, tss, new_ta);
131return new_ta;
132}
133
134//TODO see KLAM_21 for cold air flows
135
136/**
137 * @brief Greuell and Bohm katabatic flow temperature correction.
138 * This corrects the air temperature over the glaciated areas according to Greuell, W. and Bohm, R.,
139 * <i>"2 m temperatures along melting mid-latitude glaciers, and implications for the sensitivity of the mass
140 * balance to variations in temperature"</i>, J. Glaciol., 44(146), 1998.
141 * The implementation here follows the overview given by Petersen, L., Pellicciotti, F., Juszak, I., Carenzo, M.
142 * and Brock, B. <i>"Suitability of a constant air temperature lapse rate over an Alpine glacier:
143 * testing the Greuell and Bohm model as an alternative"</i>, Annals of Glaciology, 54(63), 2013.
144 *
145 * The improvements suggested by Ayala, A., Pellicciotti, F., Shea, J. M., <i>"Modeling 2m air temperatures over
146 * mountain glaciers: Exploring the influence of katabatic cooling and external warming"</i>, J. Geophys. Res. Atmos., <b>120</b>,
147 * 2015, doi:10.1002/2015JDO23137.
148 *
149 * The katabatic flows will only be considered active (ie computed and corrected for) if more than 20% of the non-glaciated
150 * pixels are snow-free and if the air temperature is higher than the surface temperature at glaciated pixels.
151 *
152 * @param hs grid containing the snow heights (used to decide if the katabatic flows are active)
153 * @param tss grid containing the surface temperatures (used to decide if the katabatic flows are active)
154 * @param ta grid containing the air temperatures that will be corrected and returned.
155 */
156void Glaciers::correctTemperatures(const mio::Grid2DObject& hs, const mio::Grid2DObject& tss, mio::Grid2DObject& ta) const
157{
158if (!dem.isSameGeolocalization(ta))
159throw IOException("The temperature grid and the provided dem don't have the same geolocalization!", AT);
160
161if (src_distance.empty())
162throw IOException("Please provide both dem and glacier map before computing temperature corrections!", AT);
163
164if (enableKatabatikFlows(hs, tss, ta, isGlacier))
165correctTemperatures(ta);
166else
167return;
168
169if (MPIControl::instance().master())
170std::cout << "[i] TA field corrected for glacier katabatic flow\n";
171}
172
173/**
174 * @brief Greuell and Bohm katabatic flow temperature correction.
175 * This method forces the air temperature correction, without checking if it is appropriate or not. This is only
176 * useful for calibration.
177 * @param ta grid containing the air temperatures that will be corrected and returned.
178 */
179void Glaciers::correctTemperatures(mio::Grid2DObject& ta) const
180{
181static const double CH = 0.002; //bulk heat transfer coefficient
182static const double Gamma = -mio::Cst::dry_adiabatique_lapse_rate;
183
184std::vector<double> X, Y;
185for (size_t idx=0; idx<isGlacier.size(); idx++) {
186if (isGlacier(idx)==1) {
187X.push_back( dem(idx) );
188Y.push_back( ta(idx) );
189}
190}
191
192const mio::Fit1D temp_fit(mio::Fit1D::NOISY_LINEAR, X, Y);
193
194for (size_t idx=0; idx<dem.size(); idx++) {
195if (dem(idx)==IOUtils::nodata) continue;
196
197const double slope = dem.slope(idx);
198const double downhill_dist = src_distance(idx)*scale;
199const double z0 = (src_altitude(idx)-dem(idx))*scale + dem(idx); //altitude where the air parcel enters the KBL
200//const double z0 = src_altitude(idx);
201
202const double T0 = temp_fit(z0) - mio::Cst::t_water_freezing_pt; //convert to celsius
203
204if (downhill_dist==IOUtils::nodata || T0==IOUtils::nodata || slope==IOUtils::nodata)
205continue;
206
207if (downhill_dist<=0.) continue;
208const double tan_slope = (z0-dem(idx)) / downhill_dist;
209const double cos_slope = Optim::invSqrt(1.+tan_slope*tan_slope); //1+tan^2 = 1/cos^2 and this is 30% faster than computing cos()
210const double Lr = KBL / CH * cos_slope; //KBL is H in the original paper
211const double b = Gamma * tan_slope;
212const double Teq = b * Lr;
213
214if (Lr>0.) {
215const double corr = ((T0 - Teq) * exp(-downhill_dist/Lr) + Teq + K*(downhill_dist/Lr) + mio::Cst::t_water_freezing_pt) - ta(idx); //in the paper: -(x-x0)/Lr
216ta(idx) = ta(idx) + 1.*corr;
217}
218}
219}
220
221
222/**
223 * @brief This distributes the flow from one cell to its neighbours
224 * @param[in] dem DEM of the domain
225 * @param[in] A upsload area of the current cell
226 * @param[in] ii first coordinate of the current cell
227 * @param[in] jj second coordinate of the current cell
228 * @param[out] flow grid that will be updated with the flow from the current cell
229 * @param[out] src_altitude altitude of the source of the flow reaching each cell
230 * @param[out] src_distance distance to the source of the flow reaching each cell
231 * @return true if it could be distributed, false if some surrounding cells were higher
232 */
233bool Glaciers::hillslope_distribute_cell(const Grid2DObject& dem, const Grid2DObject& glacier_mask, const double& A, const size_t ii, const size_t jj, Grid2DObject &flow, Grid2DObject &src_altitude, Grid2DObject &src_distance)
234{
235//set box around the current cell (ii,jj)
236const size_t jjmin = (jj>0)? jj-1 : 0;
237const size_t jjmax = (jj<dem.getNy()-1)? jj+1 : dem.getNy()-1;
238const size_t iimin = (ii>0)? ii-1 : 0;
239const size_t iimax = (ii<dem.getNx()-1)? ii+1 : dem.getNx()-1;
240
241//check that this is a top cell and get normalization factor
242double sum = 0.;
243for (size_t ll=jjmin; ll<=jjmax; ll++) {
244for (size_t mm=iimin; mm<=iimax; mm++) {
245if (glacier_mask(mm,ll)==1 && dem(mm,ll) > dem(ii,jj)) //some cells are higher and must be processed before
246return false;
247
248if (dem(ii,jj) > dem(mm,ll) && dem(mm,ll)!=IOUtils::nodata) { //(ii,jj) cell is naturally excluded
249sum += (dem(ii,jj) - dem(mm,ll)); //tan(beta) = delta_alt / L so tan(beta)*L = delta_alt
250}
251}
252}
253
254//initializa cell altitude and distance if necessary
255if (src_altitude(ii,jj) == IOUtils::nodata) {
256//ie cell is a top cell
257src_altitude(ii,jj) = dem(ii,jj);
258src_distance(ii,jj) = 0.;
259}
260
261//flow into the lower cells
262if (sum==0.) return true; //current cell is flat terrain
263const double C = A / sum;
264for (size_t ll=jjmin; ll<=jjmax; ll++) {
265for (size_t mm=iimin; mm<=iimax; mm++) {
266if (glacier_mask(mm,ll)!=1) continue; //not a glacier cell
267
268if (dem(ii,jj) > dem(mm,ll)) {
269const double flow_contrib = C * (dem(ii,jj) - dem(mm,ll));
270const double weight = flow_contrib/(flow(mm,ll)+flow_contrib);
271
272flow(mm,ll) += flow_contrib;
273
274//set altitude of the flow source
275if (src_altitude(mm,ll)==IOUtils::nodata)
276src_altitude(mm,ll) = src_altitude(ii,jj);
277else
278src_altitude(mm,ll) = weight * src_altitude(ii,jj) + (1.-weight)*src_altitude(mm,ll);
279
280//set distance of the flow source
281if (src_distance(mm,ll)==IOUtils::nodata)
282src_distance(mm,ll) = src_distance(ii,jj)+1.;
283else
284src_distance(mm,ll) = weight * (src_distance(ii,jj)+1.) + (1.-weight)*src_distance(mm,ll);
285}
286}
287}
288
289return true;
290}
291
292/**
293 * @brief This is an implementation of the multiple flow direction algorithm.
294 * This associates to each cell an upslope area that is based on the number of cells
295 * that flow into it. This follows <i>"The prediction of hillslope flow paths for distributed
296 * hydrological modelling using digital terrain models"</i>, Quinn P., Chevallier P., Planchon O.,
297 * hydrological processes, <b>5</b>, 1991, pp 59-79.
298 *
299 * In this implementation, the cells that are higher than their neighbours are computed first,
300 * and their elevation set to nodata. This is performed on the whole grid until all cells have
301 * nodata as their elevation.
302 * @param[in] glacier_mask pixels that are not glaciated are marked as IOUtils::nodata
303 */
304void Glaciers::hillslope_flow(Grid2DObject glacier_mask)
305{
306if (!dem.isSameGeolocalization(glacier_mask))
307throw IOException("The glacier mask and the provided dem don't have the same geolocalization!", AT);
308
309flowpath.set(dem, 1.);
310src_altitude.set(dem, IOUtils::nodata);
311src_distance.set(dem, IOUtils::nodata);
312
313unsigned int nr_cells_left;
314const size_t ncols = dem.getNx(), nrows = dem.getNy();
315do {
316nr_cells_left = 0;
317for (size_t jj=0; jj<nrows; jj++) {
318for (size_t ii=0; ii<ncols; ii++) {
319if (glacier_mask(ii,jj)!=1) continue; //this also skips out of domain cells
320
321const double A = flowpath(ii,jj);
322const bool distributed = hillslope_distribute_cell(dem, glacier_mask, A, ii, jj, flowpath, src_altitude, src_distance);
323if (distributed)
324glacier_mask(ii,jj) = IOUtils::nodata; //mark the cell as done
325else
326nr_cells_left++;
327}
328}
329} while (nr_cells_left>0);
330
331src_distance *= dem.cellsize;
332}

Archive Download this file

Revision: HEAD