Rev 490 | Rev 523 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
468 | werner | 2 | #include "snag.h" |
3 | #include "tree.h" |
||
4 | #include "species.h" |
||
5 | #include "globalsettings.h" |
||
6 | #include "expression.h" |
||
490 | werner | 7 | // for calculation of climate decomposition |
8 | #include "resourceunit.h" |
||
9 | #include "watercycle.h" |
||
10 | #include "climate.h" |
||
468 | werner | 11 | |
12 | /** @class Snag |
||
13 | Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools. |
||
490 | werner | 14 | Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags |
468 | werner | 15 | is subsequently forwarded to the soil sub model. |
522 | werner | 16 | Carbon is stored in three classes (depending on the size) |
468 | werner | 17 | |
18 | */ |
||
19 | // static variables |
||
522 | werner | 20 | double Snag::mDBHLower = 0.; |
21 | double Snag::mDBHHigher = 0.; |
||
22 | double Snag::mCarbonThreshold[] = {0., 0., 0.}; |
||
23 | |||
468 | werner | 24 | double CNPool::biomassCFraction = biomassCFraction; // get global from globalsettings.h |
25 | |||
522 | werner | 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 | |||
42 | |||
468 | werner | 43 | // species specific soil paramters |
44 | struct SoilParameters |
||
45 | { |
||
522 | werner | 46 | SoilParameters(): kyl(0.15), kyr(0.0807), ksw(0.015), cnFoliage(75.), cnFineroot(40.), cnWood(300.), halfLife(15.) {} |
468 | werner | 47 | double kyl; // litter decomposition rate |
48 | double kyr; // downed woody debris (dwd) decomposition rate |
||
49 | double ksw; // standing woody debris (swd) decomposition rate |
||
50 | double cnFoliage; // C/N foliage litter |
||
51 | double cnFineroot; // C/N ratio fine root |
||
52 | double cnWood; // C/N Wood: used for brances, stem and coarse root |
||
522 | werner | 53 | double halfLife; // half-life period of the given species (for calculation of SWD->DWD transition) |
468 | werner | 54 | } soilparams; |
55 | |||
56 | |||
57 | Snag::Snag() |
||
58 | { |
||
490 | werner | 59 | mRU = 0; |
476 | werner | 60 | CNPool::setCFraction(biomassCFraction); |
522 | werner | 61 | Snag::setupThresholds(mDBHLower, mDBHHigher); |
468 | werner | 62 | } |
63 | |||
490 | werner | 64 | void Snag::setup( const ResourceUnit *ru) |
468 | werner | 65 | { |
490 | werner | 66 | mRU = ru; |
67 | mClimateFactor = 0.; |
||
468 | werner | 68 | // branches |
69 | mBranchCounter=0; |
||
70 | for (int i=0;i<3;i++) { |
||
71 | mTimeSinceDeath[i] = 0.; |
||
72 | mNumberOfSnags[i] = 0.; |
||
522 | werner | 73 | mAvgDbh[i] = 0.; |
74 | mAvgHeight[i] = 0.; |
||
75 | mAvgVolume[i] = 0.; |
||
76 | mKSW[i] = 0.; |
||
77 | mCurrentKSW[i] = 0.; |
||
468 | werner | 78 | } |
475 | werner | 79 | mTotalSnagCarbon = 0.; |
468 | werner | 80 | } |
81 | |||
475 | werner | 82 | // debug outputs |
83 | QList<QVariant> Snag::debugList() |
||
84 | { |
||
85 | // list columns |
||
86 | // for three pools |
||
87 | QList<QVariant> list; |
||
88 | |||
89 | list << mTotalSnagCarbon; |
||
477 | werner | 90 | // fluxes to labile soil pool and to refractory soil pool |
91 | list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N << mSWDtoSoil.C << mSWDtoSoil.N; |
||
475 | werner | 92 | |
93 | for (int i=0;i<3;i++) { |
||
94 | // pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n" |
||
95 | list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N; |
||
96 | } |
||
97 | |||
98 | // branch pools (5 yrs) |
||
99 | list << mBranches[mBranchCounter].C << mBranches[mBranchCounter].N |
||
100 | << mBranches[(mBranchCounter+1)%5].C << mBranches[(mBranchCounter+1)%5].N |
||
101 | << mBranches[(mBranchCounter+2)%5].C << mBranches[(mBranchCounter+2)%5].N |
||
102 | << mBranches[(mBranchCounter+3)%5].C << mBranches[(mBranchCounter+3)%5].N |
||
103 | << mBranches[(mBranchCounter+4)%5].C << mBranches[(mBranchCounter+4)%5].N; |
||
104 | return list; |
||
105 | } |
||
106 | |||
468 | werner | 107 | void Snag::newYear() |
108 | { |
||
109 | for (int i=0;i<3;i++) { |
||
110 | mToSWD[i].clear(); // clear transfer pools to standing-woody-debris |
||
522 | werner | 111 | mCurrentKSW[i] = 0.; |
468 | werner | 112 | } |
113 | mLabileFlux.clear(); |
||
114 | mRefractoryFlux.clear(); |
||
476 | werner | 115 | mTotalToAtm.clear(); |
116 | mTotalToExtern.clear(); |
||
117 | mTotalIn.clear(); |
||
477 | werner | 118 | mSWDtoSoil.clear(); |
468 | werner | 119 | } |
120 | |||
490 | werner | 121 | /// calculate the dynamic climate modifier for decomposition 're' |
522 | werner | 122 | /// calculation is done on the level of ResourceUnit because the water content per day is needed. |
490 | werner | 123 | double Snag::calculateClimateFactors() |
124 | { |
||
125 | double rel_wc; |
||
126 | double ft, fw; |
||
127 | double f_sum = 0.; |
||
128 | for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day) |
||
129 | { |
||
130 | rel_wc = mRU->waterCycle()->relContent(day->day); // relative water content (0..1) of the day |
||
131 | ft = exp(308.56*(1./56.02-1./((273.+day->temperature)-227.13))); // empirical variable Q10 model of Lloyd and Taylor (1994), see also Adair et al. (2008) |
||
132 | fw = pow(1.-exp(-0.2*rel_wc),5.); // # see Standcarb for the 'stable soil' pool |
||
133 | f_sum += ft*fw; |
||
134 | } |
||
135 | // the climate factor is defined as the arithmentic annual mean value |
||
136 | mClimateFactor = f_sum / double(mRU->climate()->daysOfYear()); |
||
137 | return mClimateFactor; |
||
138 | } |
||
139 | |||
522 | werner | 140 | /// do the yearly calculation |
141 | /// see http://iland.boku.ac.at/snag+dynamics |
||
468 | werner | 142 | void Snag::processYear() |
143 | { |
||
522 | werner | 144 | mSWDtoSoil.clear(); |
477 | werner | 145 | if (isEmpty()) // nothing to do |
475 | werner | 146 | return; |
147 | |||
468 | werner | 148 | // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool |
149 | mRefractoryFlux+=mBranches[mBranchCounter]; |
||
150 | mBranches[mBranchCounter].clear(); |
||
151 | mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0. |
||
152 | |||
153 | // process standing snags. |
||
154 | // the input of the current year is in the mToSWD-Pools |
||
522 | werner | 155 | |
490 | werner | 156 | const double climate_factor_re = calculateClimateFactors(); |
468 | werner | 157 | for (int i=0;i<3;i++) { |
158 | |||
522 | werner | 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 | } |
||
475 | werner | 166 | |
522 | werner | 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 |
||
468 | werner | 171 | |
522 | werner | 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 |
||
181 | // calculate the transition probability of SWD to downed dead wood |
||
468 | werner | 182 | |
522 | werner | 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 | |||
468 | werner | 192 | // calculate flow to soil pool... |
522 | werner | 193 | mSWDtoSoil += mSWD[i] * transfer; |
194 | mRefractoryFlux += mSWD[i] * transfer; |
||
195 | mSWD[i] *= (1.-transfer); // reduce pool |
||
468 | werner | 196 | // calculate the stem number of remaining snags |
522 | werner | 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 |
||
468 | werner | 203 | mRefractoryFlux += mSWD[i]; |
522 | werner | 204 | mSWDtoSoil += mSWD[i]; |
468 | werner | 205 | mSWD[i].clear(); |
522 | werner | 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.; |
||
468 | werner | 213 | } |
522 | werner | 214 | |
468 | werner | 215 | } |
522 | werner | 216 | |
468 | werner | 217 | } |
522 | werner | 218 | // total carbon in the snag-container on the RU *after* processing is the content of the |
475 | werner | 219 | // standing woody debris pools + the branches |
220 | mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C + |
||
221 | mBranches[0].C + mBranches[1].C + mBranches[2].C + mBranches[3].C + mBranches[4].C; |
||
468 | werner | 222 | } |
223 | |||
224 | /// foliage and fineroot litter is transferred during tree growth. |
||
225 | void Snag::addTurnoverLitter(const Tree *tree, const double litter_foliage, const double litter_fineroot) |
||
226 | { |
||
227 | const SoilParameters &soil_params = soilparams; // to be updated |
||
228 | mLabileFlux.addBiomass(litter_foliage, soil_params.cnFoliage); |
||
229 | mLabileFlux.addBiomass(litter_fineroot, soil_params.cnFineroot); |
||
230 | } |
||
231 | |||
232 | /// after the death of the tree the five biomass compartments are processed. |
||
233 | void Snag::addMortality(const Tree *tree) |
||
234 | { |
||
235 | const SoilParameters &soil_params = soilparams; // to be updated |
||
236 | |||
237 | // immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool |
||
238 | // 100% of coarse root biomass: that goes to the refractory pool |
||
239 | mLabileFlux.addBiomass(tree->biomassFoliage(), soil_params.cnFoliage); |
||
240 | mLabileFlux.addBiomass(tree->biomassFineRoot(), soil_params.cnFineroot); |
||
241 | mRefractoryFlux.addBiomass(tree->biomassCoarseRoot(), soil_params.cnWood); |
||
242 | |||
243 | // branches are equally distributed over five years: |
||
477 | werner | 244 | double biomass_branch = tree->biomassBranch() * 0.2; |
468 | werner | 245 | for (int i=0;i<5; i++) |
477 | werner | 246 | mBranches[i].addBiomass(biomass_branch, soil_params.cnWood); |
468 | werner | 247 | |
248 | // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool |
||
522 | werner | 249 | int pi = poolIndex(tree->dbh()); // get right transfer pool |
250 | CNPool &swd = mToSWD[pi]; |
||
468 | werner | 251 | swd.addBiomass(tree->biomassStem(), soil_params.cnWood); |
522 | werner | 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]++; |
||
468 | werner | 271 | } |
272 | |||
273 | /// add residual biomass of 'tree' after harvesting. |
||
522 | werner | 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) |
468 | werner | 275 | /// the harvested biomass is collected. |
276 | void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction ) |
||
277 | { |
||
278 | const SoilParameters &soil_params = soilparams; // to be updated |
||
279 | |||
280 | // immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool |
||
281 | // 100% of coarse root biomass: that goes to the refractory pool |
||
282 | mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), soil_params.cnFoliage); |
||
283 | mLabileFlux.addBiomass(tree->biomassFineRoot(), soil_params.cnFineroot); |
||
284 | mRefractoryFlux.addBiomass(tree->biomassCoarseRoot(), soil_params.cnWood); |
||
285 | |||
286 | // residual branches are equally distributed over five years: |
||
287 | for (int i=0;i<5; i++) |
||
288 | mBranches[i].addBiomass(tree->biomassBranch() * remove_branch_fraction * 0.2, soil_params.cnWood); |
||
289 | |||
290 | // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool |
||
522 | werner | 291 | // TODO: what to do with harvest and stems??? I think harvested stems (that are not removed) should |
292 | // go directly to the DWD?? |
||
468 | werner | 293 | CNPool &swd = mToSWD[poolIndex(tree->dbh())]; // get right transfer pool |
294 | swd.addBiomass(tree->biomassStem() * remove_stem_fraction, soil_params.cnWood); |
||
295 | if (remove_stem_fraction < 1.) |
||
296 | mNumberOfSnags[poolIndex(tree->dbh())]++; |
||
297 | } |
||
298 |