Rev 372 | Rev 379 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
237 | werner | 2 | #include "global.h" |
3 | #include "watercycle.h" |
||
4 | #include "climate.h" |
||
5 | #include "resourceunit.h" |
||
6 | #include "species.h" |
||
239 | werner | 7 | #include "model.h" |
237 | werner | 8 | |
9 | WaterCycle::WaterCycle() |
||
10 | { |
||
266 | werner | 11 | mSoilDepth = 0; |
237 | werner | 12 | } |
13 | |||
239 | werner | 14 | void WaterCycle::setup(const ResourceUnit *ru) |
237 | werner | 15 | { |
16 | mRU = ru; |
||
17 | // get values... |
||
266 | werner | 18 | mFieldCapacity = 0.; // on top |
281 | werner | 19 | const XmlHelper &xml=GlobalSettings::instance()->settings(); |
20 | mSoilDepth = xml.valueDouble("model.site.soilDepth", 0.) * 10; // convert from cm to mm |
||
338 | werner | 21 | double pct_sand = xml.valueDouble("model.site.pctSand"); |
22 | double pct_silt = xml.valueDouble("model.site.pctSilt"); |
||
23 | double pct_clay = xml.valueDouble("model.site.pctClay"); |
||
24 | if (pct_sand + pct_silt + pct_clay != 100.) |
||
25 | 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)); |
||
26 | |||
27 | // calculate soil characteristics based on empirical functions (Schwalm & Ek, 2004) |
||
28 | // note: the variables are percentages [0..100] |
||
29 | mPsi_ref = -exp((1.54 - 0.0095*pct_sand + 0.0063*pct_silt) * log(10)) * 0.000098; // Eq. 83 |
||
30 | mPsi_koeff_b = -( 3.1 + 0.157*pct_clay - 0.003*pct_sand ); // Eq. 84 |
||
31 | mRho_ref = 0.01 * (50.5 - 0.142*pct_sand - 0.037*pct_clay); // Eq. 78 |
||
240 | werner | 32 | mCanopy.setup(); |
339 | werner | 33 | |
34 | mPermanentWiltingPoint = heightFromPsi(-4000); // maximum psi is set to a constant of -4MPa |
||
266 | werner | 35 | mFieldCapacity = heightFromPsi(-15); |
36 | mContent = mFieldCapacity; // start with full water content (in the middle of winter) |
||
339 | werner | 37 | qDebug() << "setup of water: Psi_ref (kPa)" << mPsi_ref << "Rho_ref" << mRho_ref << "coeff. b" << mPsi_koeff_b; |
237 | werner | 38 | } |
39 | |||
331 | werner | 40 | /** function to calculate the water pressure [saugspannung] for a given amount of water. |
338 | werner | 41 | returns water potential in kPa. |
42 | see http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance */ |
||
266 | werner | 43 | inline double WaterCycle::psiFromHeight(const double mm) const |
44 | { |
||
45 | // psi_x = psi_ref * ( rho_x / rho_ref) ^ b |
||
46 | if (mm<0.001) |
||
47 | return -100000000; |
||
48 | double psi_x = mPsi_ref * pow((mm / mSoilDepth / mRho_ref),mPsi_koeff_b); |
||
338 | werner | 49 | return psi_x; // pis |
266 | werner | 50 | } |
51 | |||
331 | werner | 52 | /// calculate the height of the water column for a given pressure |
338 | werner | 53 | /// return water amount in mm |
54 | /// see http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance |
||
266 | werner | 55 | inline double WaterCycle::heightFromPsi(const double psi_kpa) const |
56 | { |
||
57 | // rho_x = rho_ref * (psi_x / psi_ref)^(1/b) |
||
58 | double h = mSoilDepth * mRho_ref * pow(psi_kpa / mPsi_ref, 1./mPsi_koeff_b); |
||
59 | return h; |
||
60 | } |
||
61 | |||
331 | werner | 62 | /// get canopy characteristics of the resource unit. |
63 | /// It is important, that species-statistics are valid when this function is called (LAI)! |
||
237 | werner | 64 | void WaterCycle::getStandValues() |
65 | { |
||
246 | werner | 66 | mLAINeedle=mLAIBroadleaved=0.; |
237 | werner | 67 | mCanopyConductance=0.; |
68 | double lai; |
||
69 | foreach(const ResourceUnitSpecies &rus, mRU->ruSpecies()) { |
||
70 | lai = rus.constStatistics().leafAreaIndex(); |
||
71 | if (rus.species()->isConiferous()) |
||
72 | mLAINeedle+=lai; |
||
73 | else |
||
74 | mLAIBroadleaved+=lai; |
||
268 | werner | 75 | mCanopyConductance += rus.species()->canopyConductance() * lai; // weigh with LAI |
237 | werner | 76 | } |
77 | double total_lai = mLAIBroadleaved+mLAINeedle; |
||
78 | if (total_lai>0.) { |
||
79 | mCanopyConductance /= total_lai; |
||
268 | werner | 80 | } else { |
372 | werner | 81 | mCanopyConductance = Model::settings().boundaryLayerConductance; // defaults to 0.2 |
237 | werner | 82 | } |
368 | werner | 83 | if (total_lai < Model::settings().laiThresholdForClosedStands) { |
84 | // following Landsberg and Waring: when LAI is < 3 (default for laiThresholdForClosedStands), a linear "ramp" from 0 to 3 is assumed |
||
299 | werner | 85 | // http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance |
368 | werner | 86 | mCanopyConductance *= total_lai / Model::settings().laiThresholdForClosedStands; |
237 | werner | 87 | } |
88 | qDebug() << "WaterCycle:getStandValues: LAI needle" << mLAINeedle << "LAI Broadl:"<< mLAIBroadleaved << "weighted avg. Conductance (m/2):" << mCanopyConductance; |
||
89 | } |
||
90 | |||
367 | werner | 91 | /// calculate combined VPD and soilwaterresponse for all species |
92 | /// on the RU. This is used for the calc. of the transpiration. |
||
93 | inline double WaterCycle::calculateSoilAtmosphereResponse(const double psi_kpa, const double vpd_kpa) |
||
94 | { |
||
95 | const ResourceUnitSpecies *i; |
||
96 | const QVector<ResourceUnitSpecies>::const_iterator iend = mRU->ruSpecies().end(); |
||
97 | double min_response; |
||
98 | double total_response = 0; |
||
99 | for (i=mRU->ruSpecies().begin();i<iend;++i) { |
||
100 | if (i->LAIfactor()>0) { |
||
101 | i->speciesResponse()->soilAtmosphereResponses(psi_kpa, vpd_kpa, min_response); |
||
102 | total_response += min_response * i->LAIfactor(); |
||
103 | } |
||
104 | } |
||
377 | werner | 105 | // add an aging factor to the total response (averageAging: leaf area weighted mean aging value): |
106 | // conceptually: response = min(vpd_response, water_response)*aging |
||
107 | total_response *= mRU->averageAging(); |
||
367 | werner | 108 | return total_response; |
109 | } |
||
110 | |||
111 | |||
237 | werner | 112 | /// Main Water Cycle function. This function triggers all water related tasks for |
113 | /// one simulation year. |
||
299 | werner | 114 | /// @sa http://iland.boku.ac.at/water+cycle |
237 | werner | 115 | void WaterCycle::run() |
116 | { |
||
117 | // preparations (once a year) |
||
367 | werner | 118 | getStandValues(); // fetch canopy characteristics from iLand (including weighted average for mCanopyConductance) |
237 | werner | 119 | mCanopy.setStandParameters(mLAINeedle, |
120 | mLAIBroadleaved, |
||
121 | mCanopyConductance); |
||
122 | |||
123 | |||
124 | // main loop over all days of the year |
||
239 | werner | 125 | double prec_mm, prec_after_interception, prec_to_soil, et, excess; |
338 | werner | 126 | const Climate *climate = mRU->climate(); |
127 | const ClimateDay *day = climate->begin(); |
||
128 | const ClimateDay *end = climate->end(); |
||
237 | werner | 129 | int doy=0; |
130 | double total_excess = 0.; |
||
131 | for (; day<end; ++day, ++doy) { |
||
132 | // (1) precipitation of the day |
||
133 | prec_mm = day->preciptitation; |
||
134 | // (2) interception by the crown |
||
135 | prec_after_interception = mCanopy.flow(prec_mm, day->temperature); |
||
136 | // (3) storage in the snow pack |
||
137 | prec_to_soil = mSnowPack.flow(prec_after_interception, day->temperature); |
||
138 | // (4) add rest to soil |
||
139 | mContent += prec_to_soil; |
||
266 | werner | 140 | |
141 | excess = 0.; |
||
142 | if (mContent>mFieldCapacity) { |
||
143 | // excess water runoff |
||
144 | excess = mContent - mFieldCapacity; |
||
145 | total_excess += excess; |
||
146 | mContent = mFieldCapacity; |
||
147 | } |
||
148 | |||
237 | werner | 149 | // calculate the relative water content |
255 | werner | 150 | mRelativeContent[doy] = currentRelContent(); |
367 | werner | 151 | double current_psi = psiFromHeight(mContent); |
152 | mPsi[doy] = current_psi; |
||
335 | werner | 153 | // (5) transpiration of the vegetation (and of water intercepted in canopy) |
367 | werner | 154 | // calculate the LAI-weighted response values for soil water and vpd: |
155 | double combined_response = calculateSoilAtmosphereResponse( current_psi, day->vpd); |
||
156 | et = mCanopy.evapotranspiration3PG(day, climate->daylength_h(doy), combined_response); |
||
238 | werner | 157 | |
241 | werner | 158 | mContent -= et; // reduce content (transpiration) |
338 | werner | 159 | // add intercepted water (that is *not* evaporated) again to the soil (or add to snow if temp too low -> call to snowpack) |
160 | mContent += mSnowPack.add(mCanopy.interception(),day->temperature); |
||
241 | werner | 161 | |
338 | werner | 162 | // do not remove water below the PWP (fixed value) |
266 | werner | 163 | if (mContent<mPermanentWiltingPoint) { |
271 | werner | 164 | et -= mPermanentWiltingPoint - mContent; // reduce et (for bookkeeping) |
266 | werner | 165 | mContent = mPermanentWiltingPoint; |
237 | werner | 166 | } |
266 | werner | 167 | |
271 | werner | 168 | //DBGMODE( |
239 | werner | 169 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dWaterCycle)) { |
240 | werner | 170 | DebugList &out = GlobalSettings::instance()->debugList(day->id(), GlobalSettings::dWaterCycle); |
239 | werner | 171 | // climatic variables |
364 | werner | 172 | out << day->id() << mRU->index() << day->temperature << day->vpd << day->preciptitation << day->radiation; |
368 | werner | 173 | out << combined_response; // combined response of all species on RU (min(water, vpd)) |
239 | werner | 174 | // fluxes |
271 | werner | 175 | out << prec_after_interception << prec_to_soil << et << mCanopy.evaporationCanopy() |
176 | << mRelativeContent[doy]*mSoilDepth << mPsi[doy] << excess; |
||
239 | werner | 177 | // other states |
178 | out << mSnowPack.snowPack(); |
||
364 | werner | 179 | //special sanity check: |
180 | if (prec_to_soil>0. && mCanopy.interception()>0.) |
||
181 | if (mSnowPack.snowPack()==0. && day->preciptitation==0) |
||
182 | qDebug() << "watercontent increase without precipititaion"; |
||
239 | werner | 183 | |
184 | } |
||
271 | werner | 185 | //); // DBGMODE() |
239 | werner | 186 | |
237 | werner | 187 | } |
188 | } |
||
189 | |||
190 | |||
191 | namespace Water { |
||
192 | |||
193 | /** calculates the input/output of water that is stored in the snow pack. |
||
194 | The approach is similar to Picus 1.3 and ForestBGC (Running, 1988). |
||
195 | Returns the amount of water that exits the snowpack (precipitation, snow melt) */ |
||
196 | double SnowPack::flow(const double &preciptitation_mm, const double &temperature) |
||
197 | { |
||
198 | if (temperature>0.) { |
||
199 | if (mSnowPack==0.) |
||
200 | return preciptitation_mm; // no change |
||
201 | else { |
||
202 | // snow melts |
||
203 | const double melting_coefficient = 0.7; // mm/°C |
||
240 | werner | 204 | double melt = qMin(temperature * melting_coefficient, mSnowPack); |
205 | mSnowPack -=melt; |
||
237 | werner | 206 | return preciptitation_mm + melt; |
207 | } |
||
208 | } else { |
||
209 | // snow: |
||
210 | mSnowPack += preciptitation_mm; |
||
211 | return 0.; // no output. |
||
212 | } |
||
213 | |||
214 | } |
||
215 | |||
216 | |||
334 | werner | 217 | inline double SnowPack::add(const double &preciptitation_mm, const double &temperature) |
218 | { |
||
219 | // do nothing for temps > 0° |
||
220 | if (temperature>0.) |
||
221 | return preciptitation_mm; |
||
222 | |||
223 | // temps < 0°: add to snow |
||
224 | mSnowPack += preciptitation_mm; |
||
225 | return 0.; |
||
226 | } |
||
227 | |||
237 | werner | 228 | /** Interception in crown canopy. |
229 | Calculates the amount of preciptitation that does not reach the ground and |
||
230 | is stored in the canopy. The approach is adopted from Picus 1.3. |
||
299 | werner | 231 | Returns the amount of precipitation (mm) that surpasses the canopy layer. |
232 | @sa http://iland.boku.ac.at/water+cycle#precipitation_and_interception */ |
||
237 | werner | 233 | double Canopy::flow(const double &preciptitation_mm, const double &temperature) |
234 | { |
||
235 | // sanity checks |
||
236 | mInterception = 0.; |
||
271 | werner | 237 | mEvaporation = 0.; |
237 | werner | 238 | if (!mLAI) |
239 | return preciptitation_mm; |
||
240 | if (!preciptitation_mm) |
||
241 | return 0.; |
||
242 | double max_interception_mm=0.; // maximum interception based on the current foliage |
||
243 | double max_storage_mm=0.; // maximum storage in canopy |
||
244 | |||
245 | if (mLAINeedle>0.) { |
||
246 | // (1) calculate maximum fraction of thru-flow the crown (based on precipitation) |
||
247 | double max_flow_needle = 0.9 * sqrt(1.03 - exp(-0.055*preciptitation_mm)); |
||
248 | max_interception_mm += preciptitation_mm * (1. - max_flow_needle * mLAINeedle/mLAI); |
||
249 | // (2) calculate maximum storage potential based on the current LAI |
||
250 | double max_storage_needle = 4. * (1. - exp(-0.55*mLAINeedle) ); |
||
251 | max_storage_mm += max_storage_needle; |
||
252 | } |
||
253 | |||
254 | if (mLAIBroadleaved>0.) { |
||
255 | // (1) calculate maximum fraction of thru-flow the crown (based on precipitation) |
||
299 | werner | 256 | double max_flow_broad = 0.9 * pow(1.22 - exp(-0.055*preciptitation_mm), 0.35); |
237 | werner | 257 | max_interception_mm += preciptitation_mm * (1. - max_flow_broad * mLAIBroadleaved/mLAI); |
258 | // (2) calculate maximum storage potential based on the current LAI |
||
299 | werner | 259 | double max_storage_broad = 2. * (1. - exp(-0.5*mLAIBroadleaved) ); |
237 | werner | 260 | max_storage_mm += max_storage_broad; |
261 | } |
||
262 | |||
263 | // (3) calculate actual interception and store for evaporation calculation |
||
264 | mInterception = qMin( max_storage_mm, max_interception_mm ); |
||
265 | |||
335 | werner | 266 | // (4) limit interception with amount of precipitation |
267 | mInterception = qMin( mInterception, preciptitation_mm); |
||
268 | |||
269 | // (5) reduce preciptitaion by the amount that is intercepted by the canopy |
||
237 | werner | 270 | return preciptitation_mm - mInterception; |
271 | |||
272 | } |
||
273 | |||
239 | werner | 274 | /// sets up the canopy. fetch some global parameter values... |
275 | void Canopy::setup() |
||
237 | werner | 276 | { |
240 | werner | 277 | mAirDensity = Model::settings().airDensity; // kg / m3 |
237 | werner | 278 | } |
279 | |||
280 | void Canopy::setStandParameters(const double LAIneedle, const double LAIbroadleave, const double maxCanopyConductance) |
||
281 | { |
||
282 | mLAINeedle = LAIneedle; |
||
283 | mLAIBroadleaved=LAIbroadleave; |
||
284 | mLAI=LAIneedle+LAIbroadleave; |
||
239 | werner | 285 | mAvgMaxCanopyConductance = maxCanopyConductance; |
237 | werner | 286 | } |
287 | |||
288 | |||
289 | |||
256 | werner | 290 | /** calculate the daily evaporation/transpiration using the Penman-Monteith-Equation. |
291 | This version is based on 3PG. See the Visual Basic Code in 3PGjs.xls. |
||
292 | Returns the total sum of evaporation+transpiration in mm of the day. */ |
||
367 | werner | 293 | double Canopy::evapotranspiration3PG(const ClimateDay *climate, const double daylength_h, const double combined_response) |
256 | werner | 294 | { |
295 | double vpd_mbar = climate->vpd * 10.; // convert from kPa to mbar |
||
296 | double temperature = climate->temperature; // average temperature of the day (°C) |
||
297 | double daylength = daylength_h * 3600.; // daylength in seconds (convert from length in hours) |
||
298 | double rad = climate->radiation / daylength * 1000000; //convert from MJ/m2 (day sum) to average radiation flow W/m2 [MJ=MWs -> /s * 1,000,000 |
||
299 | |||
300 | //: Landsberg original: const double e20 = 2.2; //rate of change of saturated VP with T at 20C |
||
301 | const double VPDconv = 0.000622; //convert VPD to saturation deficit = 18/29/1000 |
||
302 | const double latent_heat = 2460000.; // Latent heat of vaporization. Energy required per unit mass of water vaporized [J kg-1] |
||
303 | |||
368 | werner | 304 | double gBL = Model::settings().boundaryLayerConductance; // boundary layer conductance |
305 | |||
306 | // canopy conductance. |
||
307 | // The species traits are weighted by LAI on the RU. |
||
308 | // maximum canopy conductance: see getStandValues() |
||
309 | // current response: see calculateSoilAtmosphereResponse(). This is basically a weighted average of min(water_response, vpd_response) for each species |
||
367 | werner | 310 | double gC = mAvgMaxCanopyConductance * combined_response; |
256 | werner | 311 | |
367 | werner | 312 | |
256 | werner | 313 | double defTerm = mAirDensity * latent_heat * (vpd_mbar * VPDconv) * gBL; |
314 | // saturation vapor pressure (Running 1988, Eq. 1) in mbar |
||
315 | double svp = 6.1078 * exp((17.269 * temperature) / (237.3 + temperature) ); |
||
316 | // the slope of svp is, thanks to http://www.wolframalpha.com/input/?i=derive+y%3D6.1078+exp+((17.269x)/(237.3%2Bx)) |
||
317 | double svp_slope = svp * ( 17.269/(237.3+temperature) - 17.269*temperature/((237.3+temperature)*(237.3+temperature)) ); |
||
318 | |||
319 | double div = (1. + svp_slope + gBL / gC); |
||
320 | double Etransp = (svp_slope * rad + defTerm) / div; |
||
335 | werner | 321 | double canopy_transpiration = Etransp / latent_heat * daylength; |
256 | werner | 322 | |
323 | if (mInterception>0.) { |
||
324 | // we assume that for evaporation from leaf surface gBL/gC -> 0 |
||
325 | double div_evap = 1 + svp_slope; |
||
326 | double evap = (svp_slope*rad + defTerm) / div_evap / latent_heat * daylength; |
||
327 | evap = qMin(evap, mInterception); |
||
335 | werner | 328 | mInterception -= evap; // reduce interception |
329 | mEvaporation = evap; // evaporation from intercepted water |
||
256 | werner | 330 | } |
335 | werner | 331 | return canopy_transpiration; |
256 | werner | 332 | } |
333 | |||
237 | werner | 334 | } // end namespace |