Rev 269 | Rev 281 | 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 |
19 | |||
269 | werner | 20 | mSoilDepth = Model::settings().waterholdingCapacity * 10; // convert from cm to mm |
240 | werner | 21 | mCanopy.setup(); |
266 | werner | 22 | mPsi_koeff_b = -3; |
23 | mPsi_ref = -0.35; // kPa |
||
24 | mRho_ref = 0.42; |
||
25 | mPermanentWiltingPoint = heightFromPsi(-1500); |
||
26 | mFieldCapacity = heightFromPsi(-15); |
||
27 | mContent = mFieldCapacity; // start with full water content (in the middle of winter) |
||
237 | werner | 28 | } |
29 | |||
266 | werner | 30 | /** functions to calculate |
31 | */ |
||
32 | inline double WaterCycle::psiFromHeight(const double mm) const |
||
33 | { |
||
34 | // psi_x = psi_ref * ( rho_x / rho_ref) ^ b |
||
35 | if (mm<0.001) |
||
36 | return -100000000; |
||
37 | double psi_x = mPsi_ref * pow((mm / mSoilDepth / mRho_ref),mPsi_koeff_b); |
||
38 | return psi_x; |
||
39 | } |
||
40 | |||
41 | inline double WaterCycle::heightFromPsi(const double psi_kpa) const |
||
42 | { |
||
43 | // rho_x = rho_ref * (psi_x / psi_ref)^(1/b) |
||
44 | double h = mSoilDepth * mRho_ref * pow(psi_kpa / mPsi_ref, 1./mPsi_koeff_b); |
||
45 | return h; |
||
46 | } |
||
47 | |||
237 | werner | 48 | void WaterCycle::getStandValues() |
49 | { |
||
246 | werner | 50 | mLAINeedle=mLAIBroadleaved=0.; |
237 | werner | 51 | mCanopyConductance=0.; |
52 | double lai; |
||
268 | werner | 53 | double psi_max = 0.; |
237 | werner | 54 | foreach(const ResourceUnitSpecies &rus, mRU->ruSpecies()) { |
55 | lai = rus.constStatistics().leafAreaIndex(); |
||
56 | if (rus.species()->isConiferous()) |
||
57 | mLAINeedle+=lai; |
||
58 | else |
||
59 | mLAIBroadleaved+=lai; |
||
268 | werner | 60 | mCanopyConductance += rus.species()->canopyConductance() * lai; // weigh with LAI |
61 | psi_max += rus.species()->psiMax() * lai; // weigh with LAI |
||
237 | werner | 62 | } |
63 | double total_lai = mLAIBroadleaved+mLAINeedle; |
||
64 | if (total_lai>0.) { |
||
65 | mCanopyConductance /= total_lai; |
||
268 | werner | 66 | mPermanentWiltingPoint = heightFromPsi(1000. * psi_max / total_lai); // (psi/total_lai) = weighted average Psi (MPa) -> convert to kPa |
67 | } else { |
||
68 | mCanopyConductance = 0.02; |
||
69 | mPermanentWiltingPoint = heightFromPsi(-1500);// -1.5MPa |
||
237 | werner | 70 | } |
240 | werner | 71 | if (total_lai < 3.) { |
237 | werner | 72 | // following Landsberg and Waring: when LAI is < 3, a linear "ramp" from 0 to 3 is assumed |
73 | mCanopyConductance *= total_lai / 3.; |
||
74 | } |
||
75 | qDebug() << "WaterCycle:getStandValues: LAI needle" << mLAINeedle << "LAI Broadl:"<< mLAIBroadleaved << "weighted avg. Conductance (m/2):" << mCanopyConductance; |
||
76 | } |
||
77 | |||
78 | /// Main Water Cycle function. This function triggers all water related tasks for |
||
79 | /// one simulation year. |
||
80 | void WaterCycle::run() |
||
81 | { |
||
82 | // preparations (once a year) |
||
83 | getStandValues(); // fetch canopy characteristics from iLand |
||
84 | mCanopy.setStandParameters(mLAINeedle, |
||
85 | mLAIBroadleaved, |
||
86 | mCanopyConductance); |
||
87 | |||
88 | |||
89 | // main loop over all days of the year |
||
239 | werner | 90 | double prec_mm, prec_after_interception, prec_to_soil, et, excess; |
237 | werner | 91 | const ClimateDay *day = mRU->climate()->begin(); |
92 | const ClimateDay *end = mRU->climate()->end(); |
||
93 | int doy=0; |
||
94 | double total_excess = 0.; |
||
95 | for (; day<end; ++day, ++doy) { |
||
96 | // (1) precipitation of the day |
||
97 | prec_mm = day->preciptitation; |
||
98 | // (2) interception by the crown |
||
99 | prec_after_interception = mCanopy.flow(prec_mm, day->temperature); |
||
100 | // (3) storage in the snow pack |
||
101 | prec_to_soil = mSnowPack.flow(prec_after_interception, day->temperature); |
||
102 | // (4) add rest to soil |
||
103 | mContent += prec_to_soil; |
||
266 | werner | 104 | |
105 | excess = 0.; |
||
106 | if (mContent>mFieldCapacity) { |
||
107 | // excess water runoff |
||
108 | excess = mContent - mFieldCapacity; |
||
109 | total_excess += excess; |
||
110 | mContent = mFieldCapacity; |
||
111 | } |
||
112 | |||
237 | werner | 113 | // calculate the relative water content |
255 | werner | 114 | mRelativeContent[doy] = currentRelContent(); |
266 | werner | 115 | mPsi[doy] = psiFromHeight(mContent); |
237 | werner | 116 | // (5) transpiration of the vegetation |
256 | werner | 117 | et = mCanopy.evapotranspiration3PG(day, mRU->climate()->daylength_h(doy)); |
238 | werner | 118 | |
241 | werner | 119 | mContent -= et; // reduce content (transpiration) |
120 | mContent += mCanopy.interception(); // add intercepted water that is *not* evaporated to the soil again |
||
121 | |||
266 | werner | 122 | if (mContent<mPermanentWiltingPoint) { |
271 | werner | 123 | et -= mPermanentWiltingPoint - mContent; // reduce et (for bookkeeping) |
266 | werner | 124 | mContent = mPermanentWiltingPoint; |
237 | werner | 125 | } |
266 | werner | 126 | |
271 | werner | 127 | //DBGMODE( |
239 | werner | 128 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dWaterCycle)) { |
240 | werner | 129 | DebugList &out = GlobalSettings::instance()->debugList(day->id(), GlobalSettings::dWaterCycle); |
239 | werner | 130 | // climatic variables |
240 | werner | 131 | out << day->id() << day->temperature << day->vpd << day->preciptitation << day->radiation; |
239 | werner | 132 | // fluxes |
271 | werner | 133 | out << prec_after_interception << prec_to_soil << et << mCanopy.evaporationCanopy() |
134 | << mRelativeContent[doy]*mSoilDepth << mPsi[doy] << excess; |
||
239 | werner | 135 | // other states |
136 | out << mSnowPack.snowPack(); |
||
137 | |||
138 | } |
||
271 | werner | 139 | //); // DBGMODE() |
239 | werner | 140 | |
237 | werner | 141 | } |
142 | } |
||
143 | |||
144 | |||
145 | namespace Water { |
||
146 | |||
147 | /** calculates the input/output of water that is stored in the snow pack. |
||
148 | The approach is similar to Picus 1.3 and ForestBGC (Running, 1988). |
||
149 | Returns the amount of water that exits the snowpack (precipitation, snow melt) */ |
||
150 | double SnowPack::flow(const double &preciptitation_mm, const double &temperature) |
||
151 | { |
||
152 | if (temperature>0.) { |
||
153 | if (mSnowPack==0.) |
||
154 | return preciptitation_mm; // no change |
||
155 | else { |
||
156 | // snow melts |
||
157 | const double melting_coefficient = 0.7; // mm/°C |
||
240 | werner | 158 | double melt = qMin(temperature * melting_coefficient, mSnowPack); |
159 | mSnowPack -=melt; |
||
237 | werner | 160 | return preciptitation_mm + melt; |
161 | } |
||
162 | } else { |
||
163 | // snow: |
||
164 | mSnowPack += preciptitation_mm; |
||
165 | return 0.; // no output. |
||
166 | } |
||
167 | |||
168 | } |
||
169 | |||
170 | |||
171 | /** Interception in crown canopy. |
||
172 | Calculates the amount of preciptitation that does not reach the ground and |
||
173 | is stored in the canopy. The approach is adopted from Picus 1.3. |
||
174 | Returns the amount of precipitation (mm) that surpasses the canopy layer. */ |
||
175 | double Canopy::flow(const double &preciptitation_mm, const double &temperature) |
||
176 | { |
||
177 | // sanity checks |
||
178 | mInterception = 0.; |
||
271 | werner | 179 | mEvaporation = 0.; |
237 | werner | 180 | if (!mLAI) |
181 | return preciptitation_mm; |
||
182 | if (!preciptitation_mm) |
||
183 | return 0.; |
||
184 | double max_interception_mm=0.; // maximum interception based on the current foliage |
||
185 | double max_storage_mm=0.; // maximum storage in canopy |
||
186 | |||
187 | if (mLAINeedle>0.) { |
||
188 | // (1) calculate maximum fraction of thru-flow the crown (based on precipitation) |
||
189 | double max_flow_needle = 0.9 * sqrt(1.03 - exp(-0.055*preciptitation_mm)); |
||
190 | max_interception_mm += preciptitation_mm * (1. - max_flow_needle * mLAINeedle/mLAI); |
||
191 | // (2) calculate maximum storage potential based on the current LAI |
||
192 | double max_storage_needle = 4. * (1. - exp(-0.55*mLAINeedle) ); |
||
193 | max_storage_mm += max_storage_needle; |
||
194 | } |
||
195 | |||
196 | if (mLAIBroadleaved>0.) { |
||
197 | // (1) calculate maximum fraction of thru-flow the crown (based on precipitation) |
||
198 | double max_flow_broad = 0.9 * sqrt(1.03 - exp(-0.055*preciptitation_mm)); |
||
199 | max_interception_mm += preciptitation_mm * (1. - max_flow_broad * mLAIBroadleaved/mLAI); |
||
200 | // (2) calculate maximum storage potential based on the current LAI |
||
201 | double max_storage_broad = 4. * (1. - exp(-0.55*mLAIBroadleaved) ); |
||
202 | max_storage_mm += max_storage_broad; |
||
203 | } |
||
204 | |||
205 | // (3) calculate actual interception and store for evaporation calculation |
||
206 | mInterception = qMin( max_storage_mm, max_interception_mm ); |
||
207 | |||
208 | // (4) reduce preciptitaion by the amount that is intercepted by the canopy |
||
241 | werner | 209 | DBG_IF(preciptitation_mm < mInterception,"water_cycle", "prec < interception"); |
237 | werner | 210 | return preciptitation_mm - mInterception; |
211 | |||
212 | } |
||
213 | |||
239 | werner | 214 | /// sets up the canopy. fetch some global parameter values... |
215 | void Canopy::setup() |
||
237 | werner | 216 | { |
240 | werner | 217 | mHeatCapacityAir = Model::settings().heatCapacityAir; // J/kg/°C |
218 | mAirDensity = Model::settings().airDensity; // kg / m3 |
||
219 | double airPressure = Model::settings().airPressure; // mbar |
||
220 | double heat_capacity_kJ = mHeatCapacityAir / 1000; // convert to: kJ/kg/°C |
||
221 | const double latent_heat_water = 2450;// kJ/kg |
||
222 | // calc psychrometric constant (in mbar/°C): see also http://en.wikipedia.org/wiki/Psychrometric_constant |
||
223 | mPsychrometricConstant = heat_capacity_kJ*airPressure/latent_heat_water/0.622; |
||
237 | werner | 224 | } |
225 | |||
226 | void Canopy::setStandParameters(const double LAIneedle, const double LAIbroadleave, const double maxCanopyConductance) |
||
227 | { |
||
228 | mLAINeedle = LAIneedle; |
||
229 | mLAIBroadleaved=LAIbroadleave; |
||
230 | mLAI=LAIneedle+LAIbroadleave; |
||
239 | werner | 231 | mAvgMaxCanopyConductance = maxCanopyConductance; |
237 | werner | 232 | } |
233 | |||
234 | /** calculate the daily evaporation/transpiration using the Penman-Monteith-Equation. |
||
235 | The application of the equation follows broadly Running (1988). |
||
236 | Returns the total sum of evaporation+transpiration in mm of the day. */ |
||
256 | werner | 237 | double Canopy::evapotranspirationBGC(const ClimateDay *climate, const double daylength_h) |
237 | werner | 238 | { |
238 | werner | 239 | double vpd_mbar = climate->vpd * 10.; // convert from kPa to mbar |
240 | double temperature = climate->temperature; // average temperature of the day (°C) |
||
241 | double daylength = daylength_h * 3600.; // daylength in seconds (convert from length in hours) |
||
242 | double rad = climate->radiation / daylength * 1000000; //convert from MJ/m2 (day sum) to average radiation flow W/m2 [MJ=MWs -> /s * 1,000,000 |
||
237 | werner | 243 | const double aerodynamic_resistance = 5.; // m/s: aerodynamic resistance of the canopy is considered being constant |
244 | const double latent_heat = 2257000.; // Latent heat of vaporization. Energy required per unit mass of water vaporized [J kg-1] |
||
245 | |||
246 | // (1) calculate some intermediaries |
||
247 | // current canopy conductance: is calculated following Landsberg & Waring |
||
240 | werner | 248 | // note: here we use vpd again in kPa. |
237 | werner | 249 | double current_canopy_conductance; |
240 | werner | 250 | current_canopy_conductance = mAvgMaxCanopyConductance * exp(-2.5 * climate->vpd); |
237 | werner | 251 | |
240 | werner | 252 | // saturation vapor pressure (Running 1988, Eq. 1) in mbar |
237 | werner | 253 | double svp = 6.1078 * exp((17.269 * temperature) / (237.3 + temperature) ); |
254 | // the slope of svp is, thanks to http://www.wolframalpha.com/input/?i=derive+y%3D6.1078+exp+((17.269x)/(237.3%2Bx)) |
||
255 | double svp_slope = svp * ( 17.269/(237.3+temperature) - 17.269*temperature/((237.3+temperature)*(237.3+temperature)) ); |
||
256 | |||
239 | werner | 257 | double et; // transpiration in mm (follows Eq.(8) of Running, 1988). |
258 | // note: RC (resistance of canopy) = 1/CC (conductance of canopy) |
||
240 | werner | 259 | double dim = svp_slope + mPsychrometricConstant*(1. + 1. / (current_canopy_conductance*aerodynamic_resistance)); |
246 | werner | 260 | double dayl = daylength*mLAI / latent_heat; |
238 | werner | 261 | double upper = svp_slope*rad + mHeatCapacityAir*mAirDensity * vpd_mbar/aerodynamic_resistance; |
262 | et = upper / dim * dayl; |
||
239 | werner | 263 | |
264 | // now calculate the evaporation from intercepted preciptitaion in the canopy: 1+rc/ra -> 1 |
||
265 | if (mInterception>0.) { |
||
271 | werner | 266 | double dim_evap = svp_slope + mPsychrometricConstant; |
267 | double pot_evap = upper / dim_evap * dayl; |
||
268 | double evap = qMin(pot_evap, mInterception); |
||
269 | mEvaporation = evap; |
||
239 | werner | 270 | } |
237 | werner | 271 | return et; |
272 | } |
||
273 | |||
274 | |||
256 | werner | 275 | /** calculate the daily evaporation/transpiration using the Penman-Monteith-Equation. |
276 | This version is based on 3PG. See the Visual Basic Code in 3PGjs.xls. |
||
277 | Returns the total sum of evaporation+transpiration in mm of the day. */ |
||
278 | double Canopy::evapotranspiration3PG(const ClimateDay *climate, const double daylength_h) |
||
279 | { |
||
280 | double vpd_mbar = climate->vpd * 10.; // convert from kPa to mbar |
||
281 | double temperature = climate->temperature; // average temperature of the day (°C) |
||
282 | double daylength = daylength_h * 3600.; // daylength in seconds (convert from length in hours) |
||
283 | double rad = climate->radiation / daylength * 1000000; //convert from MJ/m2 (day sum) to average radiation flow W/m2 [MJ=MWs -> /s * 1,000,000 |
||
284 | |||
285 | //: Landsberg original: const double e20 = 2.2; //rate of change of saturated VP with T at 20C |
||
286 | const double VPDconv = 0.000622; //convert VPD to saturation deficit = 18/29/1000 |
||
287 | const double latent_heat = 2460000.; // Latent heat of vaporization. Energy required per unit mass of water vaporized [J kg-1] |
||
288 | |||
289 | double gBL = 0.2; // boundary layer conductance |
||
290 | // canopy conductance scales linearly from 0 to LAI=3 and is constant for LAI > 3. |
||
291 | double gC = mAvgMaxCanopyConductance * exp(-2.5 * climate->vpd); |
||
292 | if (mLAI<3.) |
||
293 | gC *= mLAI / 3.; |
||
294 | |||
295 | double defTerm = mAirDensity * latent_heat * (vpd_mbar * VPDconv) * gBL; |
||
296 | // saturation vapor pressure (Running 1988, Eq. 1) in mbar |
||
297 | double svp = 6.1078 * exp((17.269 * temperature) / (237.3 + temperature) ); |
||
298 | // the slope of svp is, thanks to http://www.wolframalpha.com/input/?i=derive+y%3D6.1078+exp+((17.269x)/(237.3%2Bx)) |
||
299 | double svp_slope = svp * ( 17.269/(237.3+temperature) - 17.269*temperature/((237.3+temperature)*(237.3+temperature)) ); |
||
300 | |||
301 | double div = (1. + svp_slope + gBL / gC); |
||
302 | double Etransp = (svp_slope * rad + defTerm) / div; |
||
303 | double canopy_transiration = Etransp / latent_heat * daylength; |
||
304 | |||
305 | if (mInterception>0.) { |
||
306 | // we assume that for evaporation from leaf surface gBL/gC -> 0 |
||
307 | double div_evap = 1 + svp_slope; |
||
308 | double evap = (svp_slope*rad + defTerm) / div_evap / latent_heat * daylength; |
||
309 | evap = qMin(evap, mInterception); |
||
271 | werner | 310 | mEvaporation = evap; |
256 | werner | 311 | } |
312 | return canopy_transiration; |
||
313 | } |
||
314 | |||
237 | werner | 315 | } // end namespace |