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