Rev 808 | Rev 946 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
671 | werner | 2 | /******************************************************************************************** |
3 | ** iLand - an individual based forest landscape and disturbance model |
||
4 | ** http://iland.boku.ac.at |
||
5 | ** Copyright (C) 2009- Werner Rammer, Rupert Seidl |
||
6 | ** |
||
7 | ** This program is free software: you can redistribute it and/or modify |
||
8 | ** it under the terms of the GNU General Public License as published by |
||
9 | ** the Free Software Foundation, either version 3 of the License, or |
||
10 | ** (at your option) any later version. |
||
11 | ** |
||
12 | ** This program is distributed in the hope that it will be useful, |
||
13 | ** but WITHOUT ANY WARRANTY; without even the implied warranty of |
||
14 | ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
||
15 | ** GNU General Public License for more details. |
||
16 | ** |
||
17 | ** You should have received a copy of the GNU General Public License |
||
18 | ** along with this program. If not, see <http://www.gnu.org/licenses/>. |
||
19 | ********************************************************************************************/ |
||
20 | |||
237 | werner | 21 | #include "global.h" |
22 | #include "watercycle.h" |
||
23 | #include "climate.h" |
||
24 | #include "resourceunit.h" |
||
25 | #include "species.h" |
||
239 | werner | 26 | #include "model.h" |
808 | werner | 27 | #include "debugtimer.h" |
646 | werner | 28 | #include "modules.h" |
237 | werner | 29 | |
697 | werner | 30 | /** @class WaterCycle |
31 | @ingroup core |
||
32 | simulates the water cycle on a ResourceUnit. |
||
33 | The WaterCycle is simulated with a daily time step on the spatial level of a ResourceUnit. Related are |
||
34 | the snow module (SnowPack), and Canopy module that simulates the interception (and evaporation) of precipitation and the |
||
35 | transpiration from the canopy. |
||
36 | The WaterCycle covers the "soil water bucket". Main entry function is run(). |
||
37 | |||
38 | See http://iland.boku.ac.at/water+cycle |
||
39 | */ |
||
40 | |||
237 | werner | 41 | WaterCycle::WaterCycle() |
42 | { |
||
266 | werner | 43 | mSoilDepth = 0; |
496 | werner | 44 | mLastYear = -1; |
237 | werner | 45 | } |
46 | |||
239 | werner | 47 | void WaterCycle::setup(const ResourceUnit *ru) |
237 | werner | 48 | { |
49 | mRU = ru; |
||
50 | // get values... |
||
266 | werner | 51 | mFieldCapacity = 0.; // on top |
281 | werner | 52 | const XmlHelper &xml=GlobalSettings::instance()->settings(); |
53 | mSoilDepth = xml.valueDouble("model.site.soilDepth", 0.) * 10; // convert from cm to mm |
||
338 | werner | 54 | double pct_sand = xml.valueDouble("model.site.pctSand"); |
55 | double pct_silt = xml.valueDouble("model.site.pctSilt"); |
||
56 | double pct_clay = xml.valueDouble("model.site.pctClay"); |
||
572 | werner | 57 | if (fabs(100. - (pct_sand + pct_silt + pct_clay)) > 0.01) |
338 | werner | 58 | throw IException(QString("Setup Watercycle: soil composition percentages do not sum up to 100. Sand: %1, Silt: %2 Clay: %3").arg(pct_sand).arg(pct_silt).arg(pct_clay)); |
59 | |||
60 | // calculate soil characteristics based on empirical functions (Schwalm & Ek, 2004) |
||
61 | // note: the variables are percentages [0..100] |
||
802 | werner | 62 | mPsi_sat = -exp((1.54 - 0.0095*pct_sand + 0.0063*pct_silt) * log(10)) * 0.000098; // Eq. 83 |
338 | werner | 63 | mPsi_koeff_b = -( 3.1 + 0.157*pct_clay - 0.003*pct_sand ); // Eq. 84 |
802 | werner | 64 | mTheta_sat = 0.01 * (50.5 - 0.142*pct_sand - 0.037*pct_clay); // Eq. 78 |
240 | werner | 65 | mCanopy.setup(); |
339 | werner | 66 | |
67 | mPermanentWiltingPoint = heightFromPsi(-4000); // maximum psi is set to a constant of -4MPa |
||
379 | werner | 68 | if (xml.valueBool("model.settings.waterUseSoilSaturation",false)==false) { |
69 | mFieldCapacity = heightFromPsi(-15); |
||
70 | } else { |
||
71 | // =-EXP((1.54-0.0095* pctSand +0.0063* pctSilt)*LN(10))*0.000098 |
||
72 | double psi_sat = -exp((1.54-0.0095 * pct_sand + 0.0063*pct_silt)*log(10.))*0.000098; |
||
73 | mFieldCapacity = heightFromPsi(psi_sat); |
||
431 | werner | 74 | if (logLevelDebug()) qDebug() << "psi: saturation " << psi_sat << "field capacity:" << mFieldCapacity; |
379 | werner | 75 | } |
76 | |||
266 | werner | 77 | mContent = mFieldCapacity; // start with full water content (in the middle of winter) |
802 | werner | 78 | if (logLevelDebug()) qDebug() << "setup of water: Psi_sat (kPa)" << mPsi_sat << "Theta_sat" << mTheta_sat << "coeff. b" << mPsi_koeff_b; |
566 | werner | 79 | mCanopyConductance = 0.; |
496 | werner | 80 | mLastYear = -1; |
621 | werner | 81 | |
82 | // canopy settings |
||
83 | mCanopy.mNeedleFactor = xml.valueDouble("model.settings.interceptionStorageNeedle", 4.); |
||
84 | mCanopy.mDecidousFactor = xml.valueDouble("model.settings.interceptionStorageBroadleaf", 2.); |
||
85 | mSnowPack.mSnowTemperature = xml.valueDouble("model.settings.snowMeltTemperature", 0.); |
||
237 | werner | 86 | } |
87 | |||
331 | werner | 88 | /** function to calculate the water pressure [saugspannung] for a given amount of water. |
338 | werner | 89 | returns water potential in kPa. |
90 | see http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance */ |
||
266 | werner | 91 | inline double WaterCycle::psiFromHeight(const double mm) const |
92 | { |
||
93 | // psi_x = psi_ref * ( rho_x / rho_ref) ^ b |
||
94 | if (mm<0.001) |
||
95 | return -100000000; |
||
802 | werner | 96 | double psi_x = mPsi_sat * pow((mm / mSoilDepth / mTheta_sat),mPsi_koeff_b); |
338 | werner | 97 | return psi_x; // pis |
266 | werner | 98 | } |
99 | |||
331 | werner | 100 | /// calculate the height of the water column for a given pressure |
338 | werner | 101 | /// return water amount in mm |
102 | /// see http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance |
||
266 | werner | 103 | inline double WaterCycle::heightFromPsi(const double psi_kpa) const |
104 | { |
||
105 | // rho_x = rho_ref * (psi_x / psi_ref)^(1/b) |
||
802 | werner | 106 | double h = mSoilDepth * mTheta_sat * pow(psi_kpa / mPsi_sat, 1./mPsi_koeff_b); |
266 | werner | 107 | return h; |
108 | } |
||
109 | |||
331 | werner | 110 | /// get canopy characteristics of the resource unit. |
111 | /// It is important, that species-statistics are valid when this function is called (LAI)! |
||
237 | werner | 112 | void WaterCycle::getStandValues() |
113 | { |
||
246 | werner | 114 | mLAINeedle=mLAIBroadleaved=0.; |
237 | werner | 115 | mCanopyConductance=0.; |
553 | werner | 116 | const double ground_vegetationCC = 0.02; |
237 | werner | 117 | double lai; |
720 | werner | 118 | foreach(const ResourceUnitSpecies *rus, mRU->ruSpecies()) { |
455 | werner | 119 | lai = rus->constStatistics().leafAreaIndex(); |
120 | if (rus->species()->isConiferous()) |
||
237 | werner | 121 | mLAINeedle+=lai; |
122 | else |
||
123 | mLAIBroadleaved+=lai; |
||
455 | werner | 124 | mCanopyConductance += rus->species()->canopyConductance() * lai; // weigh with LAI |
237 | werner | 125 | } |
126 | double total_lai = mLAIBroadleaved+mLAINeedle; |
||
502 | werner | 127 | |
128 | // handle cases with LAI < 1 (use generic "ground cover characteristics" instead) |
||
723 | werner | 129 | /* The LAI used here is derived from the "stockable" area (and not the stocked area). |
130 | If the stand has gaps, the available trees are "thinned" across the whole area. Otherwise (when stocked area is used) |
||
131 | the LAI would overestimate the transpiring canopy. However, the current solution overestimates e.g. the interception. |
||
132 | If the "thinned out" LAI is below one, the rest (i.e. the gaps) are thought to be covered by ground vegetation. |
||
133 | */ |
||
502 | werner | 134 | if (total_lai<1.) { |
553 | werner | 135 | mCanopyConductance+=(ground_vegetationCC)*(1. - total_lai); |
502 | werner | 136 | total_lai = 1.; |
237 | werner | 137 | } |
502 | werner | 138 | mCanopyConductance /= total_lai; |
139 | |||
368 | werner | 140 | if (total_lai < Model::settings().laiThresholdForClosedStands) { |
141 | // following Landsberg and Waring: when LAI is < 3 (default for laiThresholdForClosedStands), a linear "ramp" from 0 to 3 is assumed |
||
299 | werner | 142 | // http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance |
368 | werner | 143 | mCanopyConductance *= total_lai / Model::settings().laiThresholdForClosedStands; |
237 | werner | 144 | } |
431 | werner | 145 | if (logLevelInfo()) qDebug() << "WaterCycle:getStandValues: LAI needle" << mLAINeedle << "LAI Broadl:"<< mLAIBroadleaved << "weighted avg. Conductance (m/2):" << mCanopyConductance; |
237 | werner | 146 | } |
147 | |||
502 | werner | 148 | /// calculate responses for ground vegetation, i.e. for "unstocked" areas. |
149 | /// this duplicates calculations done in Species. |
||
150 | /// @return Minimum of vpd and soilwater response for default |
||
151 | inline double WaterCycle::calculateBaseSoilAtmosphereResponse(const double psi_kpa, const double vpd_kpa) |
||
152 | { |
||
153 | // constant parameters used for ground vegetation: |
||
503 | werner | 154 | const double mPsiMin = 1.5; // MPa |
502 | werner | 155 | const double mRespVpdExponent = -0.6; |
156 | // see SpeciesResponse::soilAtmosphereResponses() |
||
157 | double water_resp; |
||
158 | // see Species::soilwaterResponse: |
||
159 | const double psi_mpa = psi_kpa / 1000.; // convert to MPa |
||
160 | water_resp = limit( 1. - psi_mpa / mPsiMin, 0., 1.); |
||
161 | // see species::vpdResponse |
||
162 | |||
163 | double vpd_resp; |
||
164 | vpd_resp = exp(mRespVpdExponent * vpd_kpa); |
||
165 | return qMin(water_resp, vpd_resp); |
||
166 | } |
||
167 | |||
367 | werner | 168 | /// calculate combined VPD and soilwaterresponse for all species |
169 | /// on the RU. This is used for the calc. of the transpiration. |
||
170 | inline double WaterCycle::calculateSoilAtmosphereResponse(const double psi_kpa, const double vpd_kpa) |
||
171 | { |
||
172 | double min_response; |
||
502 | werner | 173 | double total_response = 0; // LAI weighted minimum response for all speices on the RU |
174 | double total_lai_factor = 0.; |
||
720 | werner | 175 | foreach(const ResourceUnitSpecies *rus, mRU->ruSpecies()) { |
176 | if (rus->LAIfactor()>0.) { |
||
502 | werner | 177 | // retrieve the minimum of VPD / soil water response for that species |
455 | werner | 178 | rus->speciesResponse()->soilAtmosphereResponses(psi_kpa, vpd_kpa, min_response); |
179 | total_response += min_response * rus->LAIfactor(); |
||
502 | werner | 180 | total_lai_factor += rus->LAIfactor(); |
367 | werner | 181 | } |
182 | } |
||
455 | werner | 183 | |
502 | werner | 184 | if (total_lai_factor<1.) { |
185 | // the LAI is below 1: the rest is considered as "ground vegetation" |
||
186 | total_response += calculateBaseSoilAtmosphereResponse(psi_kpa, vpd_kpa) * (1. - total_lai_factor); |
||
187 | } |
||
188 | |||
377 | werner | 189 | // add an aging factor to the total response (averageAging: leaf area weighted mean aging value): |
190 | // conceptually: response = min(vpd_response, water_response)*aging |
||
503 | werner | 191 | if (total_lai_factor==1.) |
192 | total_response *= mRU->averageAging(); // no ground cover: use aging value for all LA |
||
193 | else if (total_lai_factor>0. && mRU->averageAging()>0.) |
||
194 | total_response *= (1.-total_lai_factor)*1. + (total_lai_factor * mRU->averageAging()); // between 0..1: a part of the LAI is "ground cover" (aging=1) |
||
502 | werner | 195 | |
720 | werner | 196 | DBGMODE( |
197 | if (mRU->averageAging()>1. || mRU->averageAging()<0. || total_response<0 || total_response>1.) |
||
198 | qDebug() << "water cycle: average aging invalid. aging:" << mRU->averageAging() << "total response" << total_response << "total lai factor:" << total_lai_factor; |
||
484 | werner | 199 | ); |
482 | werner | 200 | |
201 | //DBG_IF(mRU->averageAging()>1. || mRU->averageAging()<0.,"water cycle", "average aging invalid!" ); |
||
367 | werner | 202 | return total_response; |
203 | } |
||
204 | |||
205 | |||
237 | werner | 206 | /// Main Water Cycle function. This function triggers all water related tasks for |
207 | /// one simulation year. |
||
299 | werner | 208 | /// @sa http://iland.boku.ac.at/water+cycle |
237 | werner | 209 | void WaterCycle::run() |
210 | { |
||
496 | werner | 211 | // necessary? |
212 | if (GlobalSettings::instance()->currentYear() == mLastYear) |
||
213 | return; |
||
626 | werner | 214 | DebugTimer tw("water:run"); |
646 | werner | 215 | WaterCycleData add_data; |
626 | werner | 216 | |
237 | werner | 217 | // preparations (once a year) |
367 | werner | 218 | getStandValues(); // fetch canopy characteristics from iLand (including weighted average for mCanopyConductance) |
237 | werner | 219 | mCanopy.setStandParameters(mLAINeedle, |
220 | mLAIBroadleaved, |
||
221 | mCanopyConductance); |
||
222 | |||
223 | // main loop over all days of the year |
||
239 | werner | 224 | double prec_mm, prec_after_interception, prec_to_soil, et, excess; |
338 | werner | 225 | const Climate *climate = mRU->climate(); |
226 | const ClimateDay *day = climate->begin(); |
||
227 | const ClimateDay *end = climate->end(); |
||
237 | werner | 228 | int doy=0; |
229 | double total_excess = 0.; |
||
230 | for (; day<end; ++day, ++doy) { |
||
231 | // (1) precipitation of the day |
||
232 | prec_mm = day->preciptitation; |
||
233 | // (2) interception by the crown |
||
777 | werner | 234 | prec_after_interception = mCanopy.flow(prec_mm); |
237 | werner | 235 | // (3) storage in the snow pack |
236 | prec_to_soil = mSnowPack.flow(prec_after_interception, day->temperature); |
||
646 | werner | 237 | // save extra data (used by e.g. fire module) |
238 | add_data.water_to_ground[doy] = prec_to_soil; |
||
239 | add_data.snow_cover[doy] = mSnowPack.snowPack(); |
||
237 | werner | 240 | // (4) add rest to soil |
241 | mContent += prec_to_soil; |
||
266 | werner | 242 | |
243 | excess = 0.; |
||
244 | if (mContent>mFieldCapacity) { |
||
245 | // excess water runoff |
||
246 | excess = mContent - mFieldCapacity; |
||
247 | total_excess += excess; |
||
248 | mContent = mFieldCapacity; |
||
249 | } |
||
250 | |||
367 | werner | 251 | double current_psi = psiFromHeight(mContent); |
252 | mPsi[doy] = current_psi; |
||
553 | werner | 253 | |
335 | werner | 254 | // (5) transpiration of the vegetation (and of water intercepted in canopy) |
367 | werner | 255 | // calculate the LAI-weighted response values for soil water and vpd: |
680 | werner | 256 | double interception_before_transpiration = mCanopy.interception(); |
367 | werner | 257 | double combined_response = calculateSoilAtmosphereResponse( current_psi, day->vpd); |
258 | et = mCanopy.evapotranspiration3PG(day, climate->daylength_h(doy), combined_response); |
||
680 | werner | 259 | // if there is some flow from intercepted water to the ground -> add to "water_to_the_ground" |
260 | if (mCanopy.interception() < interception_before_transpiration) |
||
261 | add_data.water_to_ground[doy]+= interception_before_transpiration - mCanopy.interception(); |
||
238 | werner | 262 | |
241 | werner | 263 | mContent -= et; // reduce content (transpiration) |
338 | werner | 264 | // add intercepted water (that is *not* evaporated) again to the soil (or add to snow if temp too low -> call to snowpack) |
265 | mContent += mSnowPack.add(mCanopy.interception(),day->temperature); |
||
241 | werner | 266 | |
338 | werner | 267 | // do not remove water below the PWP (fixed value) |
266 | werner | 268 | if (mContent<mPermanentWiltingPoint) { |
271 | werner | 269 | et -= mPermanentWiltingPoint - mContent; // reduce et (for bookkeeping) |
266 | werner | 270 | mContent = mPermanentWiltingPoint; |
237 | werner | 271 | } |
266 | werner | 272 | |
546 | werner | 273 | |
271 | werner | 274 | //DBGMODE( |
239 | werner | 275 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dWaterCycle)) { |
240 | werner | 276 | DebugList &out = GlobalSettings::instance()->debugList(day->id(), GlobalSettings::dWaterCycle); |
239 | werner | 277 | // climatic variables |
605 | werner | 278 | out << day->id() << mRU->index() << mRU->id() << day->temperature << day->vpd << day->preciptitation << day->radiation; |
368 | werner | 279 | out << combined_response; // combined response of all species on RU (min(water, vpd)) |
239 | werner | 280 | // fluxes |
271 | werner | 281 | out << prec_after_interception << prec_to_soil << et << mCanopy.evaporationCanopy() |
540 | werner | 282 | << mContent << mPsi[doy] << excess; |
239 | werner | 283 | // other states |
284 | out << mSnowPack.snowPack(); |
||
364 | werner | 285 | //special sanity check: |
286 | if (prec_to_soil>0. && mCanopy.interception()>0.) |
||
287 | if (mSnowPack.snowPack()==0. && day->preciptitation==0) |
||
288 | qDebug() << "watercontent increase without precipititaion"; |
||
239 | werner | 289 | |
290 | } |
||
271 | werner | 291 | //); // DBGMODE() |
239 | werner | 292 | |
237 | werner | 293 | } |
646 | werner | 294 | // call external modules |
295 | GlobalSettings::instance()->model()->modules()->calculateWater(mRU, &add_data); |
||
496 | werner | 296 | mLastYear = GlobalSettings::instance()->currentYear(); |
297 | |||
237 | werner | 298 | } |
299 | |||
300 | |||
301 | namespace Water { |
||
302 | |||
303 | /** calculates the input/output of water that is stored in the snow pack. |
||
304 | The approach is similar to Picus 1.3 and ForestBGC (Running, 1988). |
||
305 | Returns the amount of water that exits the snowpack (precipitation, snow melt) */ |
||
306 | double SnowPack::flow(const double &preciptitation_mm, const double &temperature) |
||
307 | { |
||
621 | werner | 308 | if (temperature>mSnowTemperature) { |
237 | werner | 309 | if (mSnowPack==0.) |
310 | return preciptitation_mm; // no change |
||
311 | else { |
||
312 | // snow melts |
||
945 | werner | 313 | const double melting_coefficient = 0.7; // mm/C |
621 | werner | 314 | double melt = qMin( (temperature-mSnowTemperature) * melting_coefficient, mSnowPack); |
240 | werner | 315 | mSnowPack -=melt; |
237 | werner | 316 | return preciptitation_mm + melt; |
317 | } |
||
318 | } else { |
||
319 | // snow: |
||
320 | mSnowPack += preciptitation_mm; |
||
321 | return 0.; // no output. |
||
322 | } |
||
323 | |||
324 | } |
||
325 | |||
326 | |||
334 | werner | 327 | inline double SnowPack::add(const double &preciptitation_mm, const double &temperature) |
328 | { |
||
945 | werner | 329 | // do nothing for temps > 0 C |
621 | werner | 330 | if (temperature>mSnowTemperature) |
334 | werner | 331 | return preciptitation_mm; |
332 | |||
945 | werner | 333 | // temps < 0 C: add to snow |
334 | werner | 334 | mSnowPack += preciptitation_mm; |
335 | return 0.; |
||
336 | } |
||
337 | |||
237 | werner | 338 | /** Interception in crown canopy. |
339 | Calculates the amount of preciptitation that does not reach the ground and |
||
340 | is stored in the canopy. The approach is adopted from Picus 1.3. |
||
299 | werner | 341 | Returns the amount of precipitation (mm) that surpasses the canopy layer. |
342 | @sa http://iland.boku.ac.at/water+cycle#precipitation_and_interception */ |
||
777 | werner | 343 | double Canopy::flow(const double &preciptitation_mm) |
237 | werner | 344 | { |
345 | // sanity checks |
||
346 | mInterception = 0.; |
||
271 | werner | 347 | mEvaporation = 0.; |
237 | werner | 348 | if (!mLAI) |
349 | return preciptitation_mm; |
||
350 | if (!preciptitation_mm) |
||
351 | return 0.; |
||
352 | double max_interception_mm=0.; // maximum interception based on the current foliage |
||
945 | werner | 353 | double max_storage_mm=0.; // maximum storage in canopy (current LAI) |
354 | double max_storage_potentital = 0.; // storage capacity at very high LAI |
||
237 | werner | 355 | |
356 | if (mLAINeedle>0.) { |
||
357 | // (1) calculate maximum fraction of thru-flow the crown (based on precipitation) |
||
358 | double max_flow_needle = 0.9 * sqrt(1.03 - exp(-0.055*preciptitation_mm)); |
||
359 | max_interception_mm += preciptitation_mm * (1. - max_flow_needle * mLAINeedle/mLAI); |
||
360 | // (2) calculate maximum storage potential based on the current LAI |
||
945 | werner | 361 | // by weighing the needle/decidious storage capacity |
362 | max_storage_potentital += mNeedleFactor * mLAINeedle/mLAI; |
||
237 | werner | 363 | } |
364 | |||
365 | if (mLAIBroadleaved>0.) { |
||
366 | // (1) calculate maximum fraction of thru-flow the crown (based on precipitation) |
||
299 | werner | 367 | double max_flow_broad = 0.9 * pow(1.22 - exp(-0.055*preciptitation_mm), 0.35); |
237 | werner | 368 | max_interception_mm += preciptitation_mm * (1. - max_flow_broad * mLAIBroadleaved/mLAI); |
369 | // (2) calculate maximum storage potential based on the current LAI |
||
945 | werner | 370 | max_storage_potentital += mDecidousFactor * mLAIBroadleaved/mLAI; |
237 | werner | 371 | } |
372 | |||
945 | werner | 373 | // the extent to which the maximum stoarge capacity is exploited, depends on LAI: |
374 | max_storage_mm = max_storage_potentital * (1. - exp(-0.5 * mLAI)); |
||
375 | |||
237 | werner | 376 | // (3) calculate actual interception and store for evaporation calculation |
377 | mInterception = qMin( max_storage_mm, max_interception_mm ); |
||
378 | |||
335 | werner | 379 | // (4) limit interception with amount of precipitation |
945 | werner | 380 | mInterception = qMin( mInterception, preciptitation_mm ); |
335 | werner | 381 | |
382 | // (5) reduce preciptitaion by the amount that is intercepted by the canopy |
||
237 | werner | 383 | return preciptitation_mm - mInterception; |
384 | |||
385 | } |
||
386 | |||
239 | werner | 387 | /// sets up the canopy. fetch some global parameter values... |
388 | void Canopy::setup() |
||
237 | werner | 389 | { |
240 | werner | 390 | mAirDensity = Model::settings().airDensity; // kg / m3 |
237 | werner | 391 | } |
392 | |||
393 | void Canopy::setStandParameters(const double LAIneedle, const double LAIbroadleave, const double maxCanopyConductance) |
||
394 | { |
||
395 | mLAINeedle = LAIneedle; |
||
396 | mLAIBroadleaved=LAIbroadleave; |
||
397 | mLAI=LAIneedle+LAIbroadleave; |
||
239 | werner | 398 | mAvgMaxCanopyConductance = maxCanopyConductance; |
553 | werner | 399 | |
400 | // clear aggregation containers |
||
562 | werner | 401 | for (int i=0;i<12;++i) mET0[i]=0.; |
553 | werner | 402 | |
237 | werner | 403 | } |
404 | |||
405 | |||
406 | |||
256 | werner | 407 | /** calculate the daily evaporation/transpiration using the Penman-Monteith-Equation. |
408 | This version is based on 3PG. See the Visual Basic Code in 3PGjs.xls. |
||
409 | Returns the total sum of evaporation+transpiration in mm of the day. */ |
||
367 | werner | 410 | double Canopy::evapotranspiration3PG(const ClimateDay *climate, const double daylength_h, const double combined_response) |
256 | werner | 411 | { |
412 | double vpd_mbar = climate->vpd * 10.; // convert from kPa to mbar |
||
945 | werner | 413 | double temperature = climate->temperature; // average temperature of the day (degree C) |
256 | werner | 414 | double daylength = daylength_h * 3600.; // daylength in seconds (convert from length in hours) |
415 | double rad = climate->radiation / daylength * 1000000; //convert from MJ/m2 (day sum) to average radiation flow W/m2 [MJ=MWs -> /s * 1,000,000 |
||
416 | |||
561 | werner | 417 | // the radiation: based on linear empirical function |
418 | const double qa = -90.; |
||
419 | const double qb = 0.8; |
||
420 | double net_rad = qa + qb*rad; |
||
421 | |||
256 | werner | 422 | //: Landsberg original: const double e20 = 2.2; //rate of change of saturated VP with T at 20C |
423 | const double VPDconv = 0.000622; //convert VPD to saturation deficit = 18/29/1000 |
||
424 | const double latent_heat = 2460000.; // Latent heat of vaporization. Energy required per unit mass of water vaporized [J kg-1] |
||
425 | |||
368 | werner | 426 | double gBL = Model::settings().boundaryLayerConductance; // boundary layer conductance |
427 | |||
428 | // canopy conductance. |
||
429 | // The species traits are weighted by LAI on the RU. |
||
430 | // maximum canopy conductance: see getStandValues() |
||
431 | // current response: see calculateSoilAtmosphereResponse(). This is basically a weighted average of min(water_response, vpd_response) for each species |
||
367 | werner | 432 | double gC = mAvgMaxCanopyConductance * combined_response; |
256 | werner | 433 | |
367 | werner | 434 | |
256 | werner | 435 | double defTerm = mAirDensity * latent_heat * (vpd_mbar * VPDconv) * gBL; |
436 | |||
561 | werner | 437 | // with temperature-dependent slope of vapor pressure saturation curve |
438 | // (following Allen et al. (1998), http://www.fao.org/docrep/x0490e/x0490e07.htm#atmospheric%20parameters) |
||
439 | // svp_slope in mbar. |
||
621 | werner | 440 | //double svp_slope = 4098. * (6.1078 * exp(17.269 * temperature / (temperature + 237.3))) / ((237.3+temperature)*(237.3+temperature)); |
561 | werner | 441 | |
621 | werner | 442 | // alternatively: very simple variant (following here the original 3PG code). This |
443 | // keeps yields +- same results for summer, but slightly lower values in winter (2011/03/16) |
||
444 | double svp_slope = 2.2; |
||
445 | |||
256 | werner | 446 | double div = (1. + svp_slope + gBL / gC); |
561 | werner | 447 | double Etransp = (svp_slope * net_rad + defTerm) / div; |
335 | werner | 448 | double canopy_transpiration = Etransp / latent_heat * daylength; |
256 | werner | 449 | |
562 | werner | 450 | // calculate reference evapotranspiration |
451 | // see Adair et al 2008 |
||
452 | const double psychrometric_const = 0.0672718682328237; // kPa/degC |
||
453 | const double windspeed = 2.; // m/s |
||
454 | double net_rad_mj_day = net_rad*daylength/1000000.; // convert W/m2 again to MJ/m2*day |
||
455 | double et0_day = 0.408*svp_slope*net_rad_mj_day + psychrometric_const*900./(temperature+273.)*windspeed*climate->vpd; |
||
456 | double et0_div = svp_slope+psychrometric_const*(1.+0.34*windspeed); |
||
457 | et0_day = et0_day / et0_div; |
||
458 | mET0[climate->month-1] += et0_day; |
||
553 | werner | 459 | |
256 | werner | 460 | if (mInterception>0.) { |
461 | // we assume that for evaporation from leaf surface gBL/gC -> 0 |
||
562 | werner | 462 | double div_evap = 1. + svp_slope; |
620 | werner | 463 | double evap_canopy_potential = (svp_slope*net_rad + defTerm) / div_evap / latent_heat * daylength; |
464 | // reduce the amount of transpiration on a wet day based on the approach of |
||
465 | // Wigmosta et al (1994). see http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance |
||
466 | |||
467 | double ratio_T_E = canopy_transpiration / evap_canopy_potential; |
||
468 | double evap_canopy = qMin(evap_canopy_potential, mInterception); |
||
469 | |||
470 | // for interception -> 0, the canopy transpiration is unchanged |
||
471 | canopy_transpiration = (evap_canopy_potential - evap_canopy) * ratio_T_E; |
||
472 | |||
562 | werner | 473 | mInterception -= evap_canopy; // reduce interception |
474 | mEvaporation = evap_canopy; // evaporation from intercepted water |
||
620 | werner | 475 | |
256 | werner | 476 | } |
335 | werner | 477 | return canopy_transpiration; |
256 | werner | 478 | } |
479 | |||
237 | werner | 480 | } // end namespace |