Rev 490 | Rev 523 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 490 | Rev 522 | ||
---|---|---|---|
Line 11... | Line 11... | ||
11 | 11 | ||
12 | /** @class Snag
|
12 | /** @class Snag
|
13 | Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.
|
13 | Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.
|
14 | Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags
|
14 | Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags
|
15 | is subsequently forwarded to the soil sub model.
|
15 | is subsequently forwarded to the soil sub model.
|
- | 16 | Carbon is stored in three classes (depending on the size)
|
|
16 | 17 | ||
17 | */
|
18 | */
|
18 | // static variables
|
19 | // static variables
|
19 | double Snag::mDBHLower = 10.; |
- | |
20 | double Snag::mDBHHigher = 30.; |
- | |
- | 20 | double Snag::mDBHLower = 0.; |
|
- | 21 | double Snag::mDBHHigher = 0.; |
|
- | 22 | double Snag::mCarbonThreshold[] = {0., 0., 0.}; |
|
- | 23 | ||
21 | double CNPool::biomassCFraction = biomassCFraction; // get global from globalsettings.h |
24 | double CNPool::biomassCFraction = biomassCFraction; // get global from globalsettings.h |
- | 25 | ||
- | 26 | void Snag::setupThresholds(const double lower, const double upper) |
|
- | 27 | {
|
|
- | 28 | if (mDBHLower == lower) |
|
- | 29 | return; |
|
- | 30 | mDBHLower = lower; |
|
- | 31 | mDBHHigher = upper; |
|
- | 32 | mCarbonThreshold[0] = lower / 2.; |
|
- | 33 | mCarbonThreshold[1] = lower + (upper - lower)/2.; |
|
- | 34 | mCarbonThreshold[2] = upper + (upper - lower)/2.; |
|
- | 35 | //# threshold levels for emptying out the dbh-snag-classes
|
|
- | 36 | //# derived from Psme woody allometry, converted to C, with a threshold level set to 10%
|
|
- | 37 | //# values in kg!
|
|
- | 38 | for (int i=0;i<3;i++) |
|
- | 39 | mCarbonThreshold[i] = 0.10568*pow(mCarbonThreshold[i],2.4247)*0.5*0.1; |
|
- | 40 | }
|
|
- | 41 | ||
22 | 42 | ||
23 | // species specific soil paramters
|
43 | // species specific soil paramters
|
24 | struct SoilParameters
|
44 | struct SoilParameters
|
25 | {
|
45 | {
|
26 | SoilParameters(): kyl(0.15), kyr(0.0807), ksw(0.015), cnFoliage(75.), cnFineroot(40.), cnWood(300.) {} |
- | |
- | 46 | SoilParameters(): kyl(0.15), kyr(0.0807), ksw(0.015), cnFoliage(75.), cnFineroot(40.), cnWood(300.), halfLife(15.) {} |
|
27 | double kyl; // litter decomposition rate |
47 | double kyl; // litter decomposition rate |
28 | double kyr; // downed woody debris (dwd) decomposition rate |
48 | double kyr; // downed woody debris (dwd) decomposition rate |
29 | double ksw; // standing woody debris (swd) decomposition rate |
49 | double ksw; // standing woody debris (swd) decomposition rate |
30 | double cnFoliage; // C/N foliage litter |
50 | double cnFoliage; // C/N foliage litter |
31 | double cnFineroot; // C/N ratio fine root |
51 | double cnFineroot; // C/N ratio fine root |
32 | double cnWood; // C/N Wood: used for brances, stem and coarse root |
52 | double cnWood; // C/N Wood: used for brances, stem and coarse root |
33 | Expression pDWDformula; // expression that calculates annual transition probability for standing to downed deadwood. variable: 'tsd' (time since death) |
- | |
- | 53 | double halfLife; // half-life period of the given species (for calculation of SWD->DWD transition) |
|
34 | } soilparams; |
54 | } soilparams; |
35 | 55 | ||
36 | 56 | ||
37 | Snag::Snag() |
57 | Snag::Snag() |
38 | {
|
58 | {
|
39 | mRU = 0; |
59 | mRU = 0; |
40 | soilparams.pDWDformula.setExpression("1-1/(1+exp(-6.78+0.262*tsd))"); |
- | |
41 | CNPool::setCFraction(biomassCFraction); |
60 | CNPool::setCFraction(biomassCFraction); |
- | 61 | Snag::setupThresholds(mDBHLower, mDBHHigher); |
|
42 | }
|
62 | }
|
43 | 63 | ||
44 | void Snag::setup( const ResourceUnit *ru) |
64 | void Snag::setup( const ResourceUnit *ru) |
45 | {
|
65 | {
|
46 | mRU = ru; |
66 | mRU = ru; |
Line 48... | Line 68... | ||
48 | // branches
|
68 | // branches
|
49 | mBranchCounter=0; |
69 | mBranchCounter=0; |
50 | for (int i=0;i<3;i++) { |
70 | for (int i=0;i<3;i++) { |
51 | mTimeSinceDeath[i] = 0.; |
71 | mTimeSinceDeath[i] = 0.; |
52 | mNumberOfSnags[i] = 0.; |
72 | mNumberOfSnags[i] = 0.; |
- | 73 | mAvgDbh[i] = 0.; |
|
- | 74 | mAvgHeight[i] = 0.; |
|
- | 75 | mAvgVolume[i] = 0.; |
|
- | 76 | mKSW[i] = 0.; |
|
- | 77 | mCurrentKSW[i] = 0.; |
|
53 | }
|
78 | }
|
54 | mTotalSnagCarbon = 0.; |
79 | mTotalSnagCarbon = 0.; |
55 | }
|
80 | }
|
56 | 81 | ||
57 | // debug outputs
|
82 | // debug outputs
|
Line 67... | Line 92... | ||
67 | 92 | ||
68 | for (int i=0;i<3;i++) { |
93 | for (int i=0;i<3;i++) { |
69 | // pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n"
|
94 | // pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n"
|
70 | list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N; |
95 | list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N; |
71 | }
|
96 | }
|
72 | - | ||
73 | 97 | ||
74 | // branch pools (5 yrs)
|
98 | // branch pools (5 yrs)
|
75 | list << mBranches[mBranchCounter].C << mBranches[mBranchCounter].N |
99 | list << mBranches[mBranchCounter].C << mBranches[mBranchCounter].N |
76 | << mBranches[(mBranchCounter+1)%5].C << mBranches[(mBranchCounter+1)%5].N |
100 | << mBranches[(mBranchCounter+1)%5].C << mBranches[(mBranchCounter+1)%5].N |
77 | << mBranches[(mBranchCounter+2)%5].C << mBranches[(mBranchCounter+2)%5].N |
101 | << mBranches[(mBranchCounter+2)%5].C << mBranches[(mBranchCounter+2)%5].N |
Line 82... | Line 106... | ||
82 | 106 | ||
83 | void Snag::newYear() |
107 | void Snag::newYear() |
84 | {
|
108 | {
|
85 | for (int i=0;i<3;i++) { |
109 | for (int i=0;i<3;i++) { |
86 | mToSWD[i].clear(); // clear transfer pools to standing-woody-debris |
110 | mToSWD[i].clear(); // clear transfer pools to standing-woody-debris |
- | 111 | mCurrentKSW[i] = 0.; |
|
87 | }
|
112 | }
|
88 | mLabileFlux.clear(); |
113 | mLabileFlux.clear(); |
89 | mRefractoryFlux.clear(); |
114 | mRefractoryFlux.clear(); |
90 | mTotalToAtm.clear(); |
115 | mTotalToAtm.clear(); |
91 | mTotalToExtern.clear(); |
116 | mTotalToExtern.clear(); |
92 | mTotalIn.clear(); |
117 | mTotalIn.clear(); |
93 | mSWDtoSoil.clear(); |
118 | mSWDtoSoil.clear(); |
94 | }
|
119 | }
|
95 | 120 | ||
96 | /// calculate the dynamic climate modifier for decomposition 're'
|
121 | /// calculate the dynamic climate modifier for decomposition 're'
|
- | 122 | /// calculation is done on the level of ResourceUnit because the water content per day is needed.
|
|
97 | double Snag::calculateClimateFactors() |
123 | double Snag::calculateClimateFactors() |
98 | {
|
124 | {
|
99 | double rel_wc; |
125 | double rel_wc; |
100 | double ft, fw; |
126 | double ft, fw; |
101 | double f_sum = 0.; |
127 | double f_sum = 0.; |
Line 109... | Line 135... | ||
109 | // the climate factor is defined as the arithmentic annual mean value
|
135 | // the climate factor is defined as the arithmentic annual mean value
|
110 | mClimateFactor = f_sum / double(mRU->climate()->daysOfYear()); |
136 | mClimateFactor = f_sum / double(mRU->climate()->daysOfYear()); |
111 | return mClimateFactor; |
137 | return mClimateFactor; |
112 | }
|
138 | }
|
113 | 139 | ||
114 | // do the yearly calculation
|
- | |
115 | // see http://iland.boku.ac.at/snag+dynamics
|
- | |
- | 140 | /// do the yearly calculation
|
|
- | 141 | /// see http://iland.boku.ac.at/snag+dynamics
|
|
116 | void Snag::processYear() |
142 | void Snag::processYear() |
117 | {
|
143 | {
|
- | 144 | mSWDtoSoil.clear(); |
|
118 | if (isEmpty()) // nothing to do |
145 | if (isEmpty()) // nothing to do |
119 | return; |
146 | return; |
120 | - | ||
121 | const SoilParameters &soil_params = soilparams; // to be updated |
- | |
122 | 147 | ||
123 | // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool
|
148 | // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool
|
124 | mRefractoryFlux+=mBranches[mBranchCounter]; |
149 | mRefractoryFlux+=mBranches[mBranchCounter]; |
125 | mBranches[mBranchCounter].clear(); |
150 | mBranches[mBranchCounter].clear(); |
126 | mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0. |
151 | mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0. |
127 | 152 | ||
128 | mSWDtoSoil.clear(); |
- | |
129 | // process standing snags.
|
153 | // process standing snags.
|
130 | // the input of the current year is in the mToSWD-Pools
|
154 | // the input of the current year is in the mToSWD-Pools
|
131 | double tsd; |
- | |
- | 155 | ||
132 | const double climate_factor_re = calculateClimateFactors(); |
156 | const double climate_factor_re = calculateClimateFactors(); |
133 | for (int i=0;i<3;i++) { |
157 | for (int i=0;i<3;i++) { |
134 | // calculate 'tsd', i.e. time-since-death (see SnagDecay.xls for calculation details)
|
- | |
135 | // time_since_death = tsd_last_year*state_before / (state_before+input) + 1
|
- | |
136 | if (mSWD[i].C>0.) |
- | |
137 | tsd = mTimeSinceDeath[i]*mSWD[i].C / (mSWD[i].C+mToSWD[i].C) + 1.; |
- | |
138 | else
|
- | |
139 | tsd = 0.; |
- | |
140 | 158 | ||
141 | // update the swd-pool: move content to the SWD pool
|
- | |
142 | mSWD[i] += mToSWD[i]; |
- | |
- | 159 | // update the swd-pool with this years' input
|
|
- | 160 | if (!mToSWD[i].isEmpty()) { |
|
- | 161 | // update decay rate (apply average yearly input to the state parameters)
|
|
- | 162 | mKSW[i] = mKSW[i]*(mSWD[i].C/(mSWD[i].C+mToSWD[i].C)) + mCurrentKSW[i]*(mToSWD[i].C/(mSWD[i].C+mToSWD[i].C)); |
|
- | 163 | //move content to the SWD pool
|
|
- | 164 | mSWD[i] += mToSWD[i]; |
|
- | 165 | }
|
|
143 | 166 | ||
144 | mTimeSinceDeath[i] = tsd; |
- | |
145 | - | ||
146 | if (mSWD[i].C>0.) { |
- | |
147 | // calculate decay of SWD.
|
- | |
148 | // Eq. (1): mass (C) is lost, N remains unchanged.
|
- | |
149 | double factor = exp(-soil_params.ksw * climate_factor_re * 1.); |
- | |
150 | mSWD[i].C *= factor; |
- | |
- | 167 | if (mSWD[i].C > 0) { |
|
- | 168 | // reduce the Carbon (note: the N stays, thus the CN ratio changes)
|
|
- | 169 | // use the decay rate that is derived as a weighted average of all standing woody debris
|
|
- | 170 | mSWD[i].C *= exp(-mKSW[i] *climate_factor_re * 1. ); // 1: timestep |
|
151 | 171 | ||
- | 172 | // transition to downed woody debris
|
|
- | 173 | // update: use negative exponential decay, species parameter: half-life
|
|
- | 174 | // modified for the climatic effect on decomposition, i.e. if decomp is slower, snags stand longer and vice versa
|
|
- | 175 | // this is loosely oriented on Standcarb2 (http://andrewsforest.oregonstate.edu/pubs/webdocs/models/standcarb2.htm),
|
|
- | 176 | // where lag times for cohort transitions are linearly modified with re although here individual good or bad years have
|
|
- | 177 | // an immediate effect, the average climatic influence should come through (and it is inherently transient)
|
|
- | 178 | // note that swd.hl is species-specific, and thus a weighted average over the species in the input (=mortality)
|
|
- | 179 | // needs to be calculated, followed by a weighted update of the previous swd.hl.
|
|
- | 180 | // As weights here we use stem number, as the processes here pertain individual snags
|
|
152 | // calculate the transition probability of SWD to downed dead wood
|
181 | // calculate the transition probability of SWD to downed dead wood
|
153 | double pDWD = soil_params.pDWDformula.calculate(tsd); |
- | |
154 | pDWD = limit( pDWD * climate_factor_re, 0., 1.); // modified transition rate with climate decomp factor |
- | |
- | 182 | ||
- | 183 | double half_life = mHalfLife[i] / climate_factor_re; |
|
- | 184 | double rate = -M_LN2 / half_life; // M_LN2: math. constant |
|
- | 185 | ||
- | 186 | // higher decay rate for the class with smallest diameters
|
|
- | 187 | if (i==0) |
|
- | 188 | rate*=2.; |
|
- | 189 | ||
- | 190 | double transfer = exp(rate); |
|
- | 191 | ||
155 | // calculate flow to soil pool...
|
192 | // calculate flow to soil pool...
|
156 | mSWDtoSoil += mSWD[i] * pDWD; |
- | |
157 | mRefractoryFlux += mSWD[i] * pDWD; |
- | |
158 | mSWD[i] *= (1.-pDWD); // reduce pool |
- | |
- | 193 | mSWDtoSoil += mSWD[i] * transfer; |
|
- | 194 | mRefractoryFlux += mSWD[i] * transfer; |
|
- | 195 | mSWD[i] *= (1.-transfer); // reduce pool |
|
159 | // calculate the stem number of remaining snags
|
196 | // calculate the stem number of remaining snags
|
160 | mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - pDWD); |
- | |
161 | if (mNumberOfSnags[i] < 0.5) { |
- | |
162 | // clear the pool: add the rest to the soil
|
- | |
- | 197 | mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer); |
|
- | 198 | // if stems<0.5, empty the whole cohort into DWD, i.e. release the last bit of C and N and clear the stats
|
|
- | 199 | // also, if the Carbon of an average snag is less than 10% of the original average tree
|
|
- | 200 | // (derived from allometries for the three diameter classes), the whole cohort is emptied out to DWD
|
|
- | 201 | if (mNumberOfSnags[i] < 0.5 || mSWD[i].C / mNumberOfSnags[i] < mCarbonThreshold[i]) { |
|
- | 202 | // clear the pool: add the rest to the soil, clear statistics of the pool
|
|
163 | mRefractoryFlux += mSWD[i]; |
203 | mRefractoryFlux += mSWD[i]; |
- | 204 | mSWDtoSoil += mSWD[i]; |
|
164 | mSWD[i].clear(); |
205 | mSWD[i].clear(); |
- | 206 | mAvgDbh[i] = 0.; |
|
- | 207 | mAvgHeight[i] = 0.; |
|
- | 208 | mAvgVolume[i] = 0.; |
|
- | 209 | mKSW[i] = 0.; |
|
- | 210 | mCurrentKSW[i] = 0.; |
|
- | 211 | mHalfLife[i] = 0.; |
|
- | 212 | mTimeSinceDeath[i] = 0.; |
|
165 | }
|
213 | }
|
- | 214 | ||
166 | }
|
215 | }
|
- | 216 | ||
167 | }
|
217 | }
|
168 | // total carbon in the snag *after* processing is the content of the
|
- | |
- | 218 | // total carbon in the snag-container on the RU *after* processing is the content of the
|
|
169 | // standing woody debris pools + the branches
|
219 | // standing woody debris pools + the branches
|
170 | mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C + |
220 | mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C + |
171 | mBranches[0].C + mBranches[1].C + mBranches[2].C + mBranches[3].C + mBranches[4].C; |
221 | mBranches[0].C + mBranches[1].C + mBranches[2].C + mBranches[3].C + mBranches[4].C; |
172 | }
|
222 | }
|
173 | 223 | ||
Line 194... | Line 244... | ||
194 | double biomass_branch = tree->biomassBranch() * 0.2; |
244 | double biomass_branch = tree->biomassBranch() * 0.2; |
195 | for (int i=0;i<5; i++) |
245 | for (int i=0;i<5; i++) |
196 | mBranches[i].addBiomass(biomass_branch, soil_params.cnWood); |
246 | mBranches[i].addBiomass(biomass_branch, soil_params.cnWood); |
197 | 247 | ||
198 | // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
|
248 | // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
|
199 | CNPool &swd = mToSWD[poolIndex(tree->dbh())]; // get right transfer pool |
- | |
- | 249 | int pi = poolIndex(tree->dbh()); // get right transfer pool |
|
- | 250 | CNPool &swd = mToSWD[pi]; |
|
200 | swd.addBiomass(tree->biomassStem(), soil_params.cnWood); |
251 | swd.addBiomass(tree->biomassStem(), soil_params.cnWood); |
201 | mNumberOfSnags[poolIndex(tree->dbh())]++; |
- | |
- | 252 | ||
- | 253 | // update statistics - stemnumber-weighted averages
|
|
- | 254 | // note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)
|
|
- | 255 | double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers) |
|
- | 256 | double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1). |
|
- | 257 | mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new; |
|
- | 258 | mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new; |
|
- | 259 | mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new; |
|
- | 260 | mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new; |
|
- | 261 | mHalfLife[pi] = mHalfLife[pi]*p_old + soil_params.halfLife * p_new; |
|
- | 262 | ||
- | 263 | // average the decay rate (ksw); this is done based on the carbon content
|
|
- | 264 | // aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)
|
|
- | 265 | if (tree->biomassStem()==0) |
|
- | 266 | throw IException("Snag::addMortality: tree without stem biomass!!"); |
|
- | 267 | p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction); |
|
- | 268 | p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction); |
|
- | 269 | mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + soil_params.ksw * p_new; |
|
- | 270 | mNumberOfSnags[pi]++; |
|
202 | }
|
271 | }
|
203 | 272 | ||
204 | /// add residual biomass of 'tree' after harvesting.
|
273 | /// add residual biomass of 'tree' after harvesting.
|
205 | /// remove_(stem, branch, foliage)_fraction: percentage of biomass compartment that is *removed* by the harvest operation.
|
- | |
- | 274 | /// remove_{stem, branch, foliage}_fraction: percentage of biomass compartment that is *removed* by the harvest operation (i.e.: not to stay in the system)
|
|
206 | /// the harvested biomass is collected.
|
275 | /// the harvested biomass is collected.
|
207 | void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction ) |
276 | void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction ) |
208 | {
|
277 | {
|
209 | const SoilParameters &soil_params = soilparams; // to be updated |
278 | const SoilParameters &soil_params = soilparams; // to be updated |
210 | 279 | ||
Line 217... | Line 286... | ||
217 | // residual branches are equally distributed over five years:
|
286 | // residual branches are equally distributed over five years:
|
218 | for (int i=0;i<5; i++) |
287 | for (int i=0;i<5; i++) |
219 | mBranches[i].addBiomass(tree->biomassBranch() * remove_branch_fraction * 0.2, soil_params.cnWood); |
288 | mBranches[i].addBiomass(tree->biomassBranch() * remove_branch_fraction * 0.2, soil_params.cnWood); |
220 | 289 | ||
221 | // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
|
290 | // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
|
- | 291 | // TODO: what to do with harvest and stems??? I think harvested stems (that are not removed) should
|
|
- | 292 | // go directly to the DWD??
|
|
222 | CNPool &swd = mToSWD[poolIndex(tree->dbh())]; // get right transfer pool |
293 | CNPool &swd = mToSWD[poolIndex(tree->dbh())]; // get right transfer pool |
223 | swd.addBiomass(tree->biomassStem() * remove_stem_fraction, soil_params.cnWood); |
294 | swd.addBiomass(tree->biomassStem() * remove_stem_fraction, soil_params.cnWood); |
224 | if (remove_stem_fraction < 1.) |
295 | if (remove_stem_fraction < 1.) |
225 | mNumberOfSnags[poolIndex(tree->dbh())]++; |
296 | mNumberOfSnags[poolIndex(tree->dbh())]++; |
226 | }
|
297 | }
|