Rev 595 | Rev 607 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 595 | Rev 605 | ||
---|---|---|---|
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/snag.cpp': |
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/snag.cpp': |
2 | #include "snag.h"
|
2 | #include "snag.h"
|
3 | #include "tree.h"
|
3 | #include "tree.h"
|
4 | #include "species.h"
|
4 | #include "species.h"
|
5 | #include "globalsettings.h"
|
5 | #include "globalsettings.h"
|
6 | #include "expression.h"
|
6 | #include "expression.h"
|
7 | // for calculation of climate decomposition
|
7 | // for calculation of climate decomposition
|
8 | #include "resourceunit.h"
|
8 | #include "resourceunit.h"
|
9 | #include "watercycle.h"
|
9 | #include "watercycle.h"
|
10 | #include "climate.h"
|
10 | #include "climate.h"
|
11 | #include "model.h"
|
11 | #include "model.h"
|
12 | 12 | ||
13 | /** @class Snag
|
13 | /** @class Snag
|
14 | Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.
|
14 | Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.
|
15 | Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags
|
15 | Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags
|
16 | is subsequently forwarded to the soil sub model.
|
16 | is subsequently forwarded to the soil sub model.
|
17 | Carbon is stored in three classes (depending on the size)
|
17 | Carbon is stored in three classes (depending on the size)
|
18 | The Snag dynamics class uses the following species parameter:
|
18 | The Snag dynamics class uses the following species parameter:
|
19 | cnFoliage, cnFineroot, cnWood, snagHalflife, snagKSW
|
19 | cnFoliage, cnFineroot, cnWood, snagHalflife, snagKSW
|
20 | 20 | ||
21 | */
|
21 | */
|
22 | // static variables
|
22 | // static variables
|
23 | double Snag::mDBHLower = -1.; |
23 | double Snag::mDBHLower = -1.; |
24 | double Snag::mDBHHigher = 0.; |
24 | double Snag::mDBHHigher = 0.; |
25 | double Snag::mCarbonThreshold[] = {0., 0., 0.}; |
25 | double Snag::mCarbonThreshold[] = {0., 0., 0.}; |
26 | 26 | ||
27 | double CNPair::biomassCFraction = biomassCFraction; // get global from globalsettings.h |
27 | double CNPair::biomassCFraction = biomassCFraction; // get global from globalsettings.h |
28 | 28 | ||
29 | /// add biomass and weigh the parameter_value with the current C-content of the pool
|
29 | /// add biomass and weigh the parameter_value with the current C-content of the pool
|
30 | void CNPool::addBiomass(const double biomass, const double CNratio, const double parameter_value) |
30 | void CNPool::addBiomass(const double biomass, const double CNratio, const double parameter_value) |
31 | {
|
31 | {
|
32 | if (biomass==0.) |
32 | if (biomass==0.) |
33 | return; |
33 | return; |
34 | double new_c = biomass*biomassCFraction; |
34 | double new_c = biomass*biomassCFraction; |
35 | double p_old = C / (new_c + C); |
35 | double p_old = C / (new_c + C); |
36 | mParameter = mParameter*p_old + parameter_value*(1.-p_old); |
36 | mParameter = mParameter*p_old + parameter_value*(1.-p_old); |
37 | CNPair::addBiomass(biomass, CNratio); |
37 | CNPair::addBiomass(biomass, CNratio); |
38 | }
|
38 | }
|
39 | 39 | ||
40 | // increase pool (and weigh the value)
|
40 | // increase pool (and weigh the value)
|
41 | void CNPool::operator+=(const CNPool &s) |
41 | void CNPool::operator+=(const CNPool &s) |
42 | {
|
42 | {
|
43 | if (s.C==0.) |
43 | if (s.C==0.) |
44 | return; |
44 | return; |
45 | mParameter = parameter(s); // calculate weighted parameter |
45 | mParameter = parameter(s); // calculate weighted parameter |
46 | C+=s.C; |
46 | C+=s.C; |
47 | N+=s.N; |
47 | N+=s.N; |
48 | }
|
48 | }
|
49 | 49 | ||
50 | double CNPool::parameter(const CNPool &s) const |
50 | double CNPool::parameter(const CNPool &s) const |
51 | {
|
51 | {
|
52 | if (s.C==0.) |
52 | if (s.C==0.) |
53 | return parameter(); |
53 | return parameter(); |
54 | double p_old = C / (s.C + C); |
54 | double p_old = C / (s.C + C); |
55 | double result = mParameter*p_old + s.parameter()*(1.-p_old); |
55 | double result = mParameter*p_old + s.parameter()*(1.-p_old); |
56 | return result; |
56 | return result; |
57 | }
|
57 | }
|
58 | 58 | ||
59 | 59 | ||
60 | void Snag::setupThresholds(const double lower, const double upper) |
60 | void Snag::setupThresholds(const double lower, const double upper) |
61 | {
|
61 | {
|
62 | if (mDBHLower == lower) |
62 | if (mDBHLower == lower) |
63 | return; |
63 | return; |
64 | mDBHLower = lower; |
64 | mDBHLower = lower; |
65 | mDBHHigher = upper; |
65 | mDBHHigher = upper; |
66 | mCarbonThreshold[0] = lower / 2.; |
66 | mCarbonThreshold[0] = lower / 2.; |
67 | mCarbonThreshold[1] = lower + (upper - lower)/2.; |
67 | mCarbonThreshold[1] = lower + (upper - lower)/2.; |
68 | mCarbonThreshold[2] = upper + (upper - lower)/2.; |
68 | mCarbonThreshold[2] = upper + (upper - lower)/2.; |
69 | //# threshold levels for emptying out the dbh-snag-classes
|
69 | //# threshold levels for emptying out the dbh-snag-classes
|
70 | //# derived from Psme woody allometry, converted to C, with a threshold level set to 10%
|
70 | //# derived from Psme woody allometry, converted to C, with a threshold level set to 10%
|
71 | //# values in kg!
|
71 | //# values in kg!
|
72 | for (int i=0;i<3;i++) |
72 | for (int i=0;i<3;i++) |
73 | mCarbonThreshold[i] = 0.10568*pow(mCarbonThreshold[i],2.4247)*0.5*0.1; |
73 | mCarbonThreshold[i] = 0.10568*pow(mCarbonThreshold[i],2.4247)*0.5*0.1; |
74 | }
|
74 | }
|
75 | 75 | ||
76 | 76 | ||
77 | Snag::Snag() |
77 | Snag::Snag() |
78 | {
|
78 | {
|
79 | mRU = 0; |
79 | mRU = 0; |
80 | CNPair::setCFraction(biomassCFraction); |
80 | CNPair::setCFraction(biomassCFraction); |
81 | }
|
81 | }
|
82 | 82 | ||
83 | void Snag::setup( const ResourceUnit *ru) |
83 | void Snag::setup( const ResourceUnit *ru) |
84 | {
|
84 | {
|
85 | mRU = ru; |
85 | mRU = ru; |
86 | mClimateFactor = 0.; |
86 | mClimateFactor = 0.; |
87 | // branches
|
87 | // branches
|
88 | mBranchCounter=0; |
88 | mBranchCounter=0; |
89 | for (int i=0;i<3;i++) { |
89 | for (int i=0;i<3;i++) { |
90 | mTimeSinceDeath[i] = 0.; |
90 | mTimeSinceDeath[i] = 0.; |
91 | mNumberOfSnags[i] = 0.; |
91 | mNumberOfSnags[i] = 0.; |
92 | mAvgDbh[i] = 0.; |
92 | mAvgDbh[i] = 0.; |
93 | mAvgHeight[i] = 0.; |
93 | mAvgHeight[i] = 0.; |
94 | mAvgVolume[i] = 0.; |
94 | mAvgVolume[i] = 0.; |
95 | mKSW[i] = 0.; |
95 | mKSW[i] = 0.; |
96 | mCurrentKSW[i] = 0.; |
96 | mCurrentKSW[i] = 0.; |
97 | mHalfLife[i] = 0.; |
97 | mHalfLife[i] = 0.; |
98 | }
|
98 | }
|
99 | mTotalSnagCarbon = 0.; |
99 | mTotalSnagCarbon = 0.; |
100 | if (mDBHLower<=0) |
100 | if (mDBHLower<=0) |
101 | throw IException("Snag::setupThresholds() not called or called with invalid parameters."); |
101 | throw IException("Snag::setupThresholds() not called or called with invalid parameters."); |
102 | 102 | ||
103 | // Inital values from XML file
|
103 | // Inital values from XML file
|
104 | XmlHelper xml=GlobalSettings::instance()->settings(); |
104 | XmlHelper xml=GlobalSettings::instance()->settings(); |
105 | double kyr = xml.valueDouble("model.site.youngRefractoryDecompRate", -1); |
105 | double kyr = xml.valueDouble("model.site.youngRefractoryDecompRate", -1); |
106 | // put carbon of snags to the middle size class
|
106 | // put carbon of snags to the middle size class
|
107 | xml.setCurrentNode("model.initialization.snags"); |
107 | xml.setCurrentNode("model.initialization.snags"); |
108 | mSWD[1].C = xml.valueDouble(".swdC"); |
108 | mSWD[1].C = xml.valueDouble(".swdC"); |
109 | mSWD[1].N = mSWD[1].C / xml.valueDouble(".swdCN", 50.); |
109 | mSWD[1].N = mSWD[1].C / xml.valueDouble(".swdCN", 50.); |
110 | mSWD[1].setParameter(kyr); |
110 | mSWD[1].setParameter(kyr); |
111 | mKSW[1] = xml.valueDouble(".swdDecompRate"); |
111 | mKSW[1] = xml.valueDouble(".swdDecompRate"); |
112 | mNumberOfSnags[1] = xml.valueDouble(".swdCount"); |
112 | mNumberOfSnags[1] = xml.valueDouble(".swdCount"); |
113 | mHalfLife[1] = xml.valueDouble(".swdHalfLife"); |
113 | mHalfLife[1] = xml.valueDouble(".swdHalfLife"); |
114 | // and for the Branch/coarse root pools: split the init value into five chunks
|
114 | // and for the Branch/coarse root pools: split the init value into five chunks
|
115 | CNPool other(xml.valueDouble(".otherC"), xml.valueDouble(".otherC")/xml.valueDouble(".otherCN", 50.), kyr ); |
115 | CNPool other(xml.valueDouble(".otherC"), xml.valueDouble(".otherC")/xml.valueDouble(".otherCN", 50.), kyr ); |
116 | 116 | ||
117 | mTotalSnagCarbon = other.C + mSWD[1].C; |
117 | mTotalSnagCarbon = other.C + mSWD[1].C; |
118 | 118 | ||
119 | other *= 0.2; |
119 | other *= 0.2; |
120 | for (int i=0;i<5;i++) |
120 | for (int i=0;i<5;i++) |
121 | mOtherWood[i] = other; |
121 | mOtherWood[i] = other; |
122 | }
|
122 | }
|
123 | 123 | ||
124 | // debug outputs
|
124 | // debug outputs
|
125 | QList<QVariant> Snag::debugList() |
125 | QList<QVariant> Snag::debugList() |
126 | {
|
126 | {
|
127 | // list columns
|
127 | // list columns
|
128 | // for three pools
|
128 | // for three pools
|
129 | QList<QVariant> list; |
129 | QList<QVariant> list; |
130 | 130 | ||
131 | // totals
|
131 | // totals
|
132 | list << mTotalSnagCarbon << mTotalIn.C << mTotalToAtm.C << mSWDtoSoil.C << mSWDtoSoil.N; |
132 | list << mTotalSnagCarbon << mTotalIn.C << mTotalToAtm.C << mSWDtoSoil.C << mSWDtoSoil.N; |
133 | // fluxes to labile soil pool and to refractory soil pool
|
133 | // fluxes to labile soil pool and to refractory soil pool
|
134 | list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N; |
134 | list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N; |
135 | 135 | ||
136 | for (int i=0;i<3;i++) { |
136 | for (int i=0;i<3;i++) { |
137 | // pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n"
|
137 | // pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n"
|
138 | list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N; |
138 | list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N; |
139 | list << mAvgDbh[i] << mAvgHeight[i] << mAvgVolume[i]; |
139 | list << mAvgDbh[i] << mAvgHeight[i] << mAvgVolume[i]; |
140 | }
|
140 | }
|
141 | 141 | ||
142 | // branch/coarse wood pools (5 yrs)
|
142 | // branch/coarse wood pools (5 yrs)
|
143 | for (int i=0;i<5;i++) { |
143 | for (int i=0;i<5;i++) { |
144 | list << mOtherWood[i].C << mOtherWood[i].N; |
144 | list << mOtherWood[i].C << mOtherWood[i].N; |
145 | }
|
145 | }
|
146 | // list << mOtherWood[mBranchCounter].C << mOtherWood[mBranchCounter].N
|
146 | // list << mOtherWood[mBranchCounter].C << mOtherWood[mBranchCounter].N
|
147 | // << mOtherWood[(mBranchCounter+1)%5].C << mOtherWood[(mBranchCounter+1)%5].N
|
147 | // << mOtherWood[(mBranchCounter+1)%5].C << mOtherWood[(mBranchCounter+1)%5].N
|
148 | // << mOtherWood[(mBranchCounter+2)%5].C << mOtherWood[(mBranchCounter+2)%5].N
|
148 | // << mOtherWood[(mBranchCounter+2)%5].C << mOtherWood[(mBranchCounter+2)%5].N
|
149 | // << mOtherWood[(mBranchCounter+3)%5].C << mOtherWood[(mBranchCounter+3)%5].N
|
149 | // << mOtherWood[(mBranchCounter+3)%5].C << mOtherWood[(mBranchCounter+3)%5].N
|
150 | // << mOtherWood[(mBranchCounter+4)%5].C << mOtherWood[(mBranchCounter+4)%5].N;
|
150 | // << mOtherWood[(mBranchCounter+4)%5].C << mOtherWood[(mBranchCounter+4)%5].N;
|
151 | return list; |
151 | return list; |
152 | }
|
152 | }
|
153 | 153 | ||
154 | void Snag::newYear() |
154 | void Snag::newYear() |
155 | {
|
155 | {
|
156 | for (int i=0;i<3;i++) { |
156 | for (int i=0;i<3;i++) { |
157 | mToSWD[i].clear(); // clear transfer pools to standing-woody-debris |
157 | mToSWD[i].clear(); // clear transfer pools to standing-woody-debris |
158 | mCurrentKSW[i] = 0.; |
158 | mCurrentKSW[i] = 0.; |
159 | }
|
159 | }
|
160 | mLabileFlux.clear(); |
160 | mLabileFlux.clear(); |
161 | mRefractoryFlux.clear(); |
161 | mRefractoryFlux.clear(); |
162 | mTotalToAtm.clear(); |
162 | mTotalToAtm.clear(); |
163 | mTotalToExtern.clear(); |
163 | mTotalToExtern.clear(); |
164 | mTotalIn.clear(); |
164 | mTotalIn.clear(); |
165 | mSWDtoSoil.clear(); |
165 | mSWDtoSoil.clear(); |
166 | }
|
166 | }
|
167 | 167 | ||
168 | /// calculate the dynamic climate modifier for decomposition 're'
|
168 | /// calculate the dynamic climate modifier for decomposition 're'
|
169 | /// calculation is done on the level of ResourceUnit because the water content per day is needed.
|
169 | /// calculation is done on the level of ResourceUnit because the water content per day is needed.
|
170 | double Snag::calculateClimateFactors() |
170 | double Snag::calculateClimateFactors() |
171 | {
|
171 | {
|
172 | double ft, fw; |
172 | double ft, fw; |
173 | double f_sum = 0.; |
173 | double f_sum = 0.; |
174 | int iday=0; |
174 | int iday=0; |
175 | // calculate the water-factor for each month (see Adair et al 2008)
|
175 | // calculate the water-factor for each month (see Adair et al 2008)
|
176 | double fw_month[12]; |
176 | double fw_month[12]; |
177 | double ratio; |
177 | double ratio; |
178 | for (int m=0;m<12;m++) { |
178 | for (int m=0;m<12;m++) { |
179 | if (mRU->waterCycle()->referenceEvapotranspiration()[m]>0.) |
179 | if (mRU->waterCycle()->referenceEvapotranspiration()[m]>0.) |
180 | ratio = mRU->climate()->precipitationMonth()[m] / mRU->waterCycle()->referenceEvapotranspiration()[m]; |
180 | ratio = mRU->climate()->precipitationMonth()[m] / mRU->waterCycle()->referenceEvapotranspiration()[m]; |
181 | else
|
181 | else
|
182 | ratio = 0; |
182 | ratio = 0; |
183 | fw_month[m] = 1. / (1. + 30.*exp(-8.5*ratio)); |
183 | fw_month[m] = 1. / (1. + 30.*exp(-8.5*ratio)); |
184 | if (logLevelDebug()) qDebug() <<"month"<< m << "PET" << mRU->waterCycle()->referenceEvapotranspiration()[m] << "prec" <<mRU->climate()->precipitationMonth()[m]; |
184 | if (logLevelDebug()) qDebug() <<"month"<< m << "PET" << mRU->waterCycle()->referenceEvapotranspiration()[m] << "prec" <<mRU->climate()->precipitationMonth()[m]; |
185 | }
|
185 | }
|
186 | 186 | ||
187 | for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day, ++iday) |
187 | for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day, ++iday) |
188 | {
|
188 | {
|
189 | 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) |
189 | 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) |
190 | fw = fw_month[day->month-1]; |
190 | fw = fw_month[day->month-1]; |
191 | 191 | ||
192 | f_sum += ft*fw; |
192 | f_sum += ft*fw; |
193 | }
|
193 | }
|
194 | // the climate factor is defined as the arithmentic annual mean value
|
194 | // the climate factor is defined as the arithmentic annual mean value
|
195 | mClimateFactor = f_sum / double(mRU->climate()->daysOfYear()); |
195 | mClimateFactor = f_sum / double(mRU->climate()->daysOfYear()); |
196 | return mClimateFactor; |
196 | return mClimateFactor; |
197 | }
|
197 | }
|
198 | 198 | ||
199 | /// do the yearly calculation
|
199 | /// do the yearly calculation
|
200 | /// see http://iland.boku.ac.at/snag+dynamics
|
200 | /// see http://iland.boku.ac.at/snag+dynamics
|
201 | void Snag::calculateYear() |
201 | void Snag::calculateYear() |
202 | {
|
202 | {
|
203 | mSWDtoSoil.clear(); |
203 | mSWDtoSoil.clear(); |
204 | const double climate_factor_re = calculateClimateFactors(); // calculate anyway, because also the soil module needs it (and currently one can have Snag and Soil only as a couple) |
204 | const double climate_factor_re = calculateClimateFactors(); // calculate anyway, because also the soil module needs it (and currently one can have Snag and Soil only as a couple) |
205 | if (isEmpty()) // nothing to do |
205 | if (isEmpty()) // nothing to do |
206 | return; |
206 | return; |
207 | 207 | ||
208 | // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool
|
208 | // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool
|
209 | mRefractoryFlux+=mOtherWood[mBranchCounter]; |
209 | mRefractoryFlux+=mOtherWood[mBranchCounter]; |
210 | 210 | ||
211 | mOtherWood[mBranchCounter].clear(); |
211 | mOtherWood[mBranchCounter].clear(); |
212 | mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0. |
212 | mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0. |
213 | // decay of branches/coarse roots
|
213 | // decay of branches/coarse roots
|
214 | for (int i=0;i<5;i++) { |
214 | for (int i=0;i<5;i++) { |
215 | if (mOtherWood[i].C>0.) { |
215 | if (mOtherWood[i].C>0.) { |
216 | double survive_rate = exp(- climate_factor_re * mOtherWood[i].parameter() ); // parameter: the "kyr" value... |
216 | double survive_rate = exp(- climate_factor_re * mOtherWood[i].parameter() ); // parameter: the "kyr" value... |
217 | mOtherWood[i].C *= survive_rate; |
217 | mOtherWood[i].C *= survive_rate; |
218 | }
|
218 | }
|
219 | }
|
219 | }
|
220 | 220 | ||
221 | // process standing snags.
|
221 | // process standing snags.
|
222 | // the input of the current year is in the mToSWD-Pools
|
222 | // the input of the current year is in the mToSWD-Pools
|
223 | for (int i=0;i<3;i++) { |
223 | for (int i=0;i<3;i++) { |
224 | 224 | ||
225 | // update the swd-pool with this years' input
|
225 | // update the swd-pool with this years' input
|
226 | if (!mToSWD[i].isEmpty()) { |
226 | if (!mToSWD[i].isEmpty()) { |
227 | // update decay rate (apply average yearly input to the state parameters)
|
227 | // update decay rate (apply average yearly input to the state parameters)
|
228 | mKSW[i] = mKSW[i]*(mSWD[i].C/(mSWD[i].C+mToSWD[i].C)) + mCurrentKSW[i]*(mToSWD[i].C/(mSWD[i].C+mToSWD[i].C)); |
228 | mKSW[i] = mKSW[i]*(mSWD[i].C/(mSWD[i].C+mToSWD[i].C)) + mCurrentKSW[i]*(mToSWD[i].C/(mSWD[i].C+mToSWD[i].C)); |
229 | //move content to the SWD pool
|
229 | //move content to the SWD pool
|
230 | mSWD[i] += mToSWD[i]; |
230 | mSWD[i] += mToSWD[i]; |
231 | }
|
231 | }
|
232 | 232 | ||
233 | if (mSWD[i].C > 0) { |
233 | if (mSWD[i].C > 0) { |
234 | // reduce the Carbon (note: the N stays, thus the CN ratio changes)
|
234 | // reduce the Carbon (note: the N stays, thus the CN ratio changes)
|
235 | // use the decay rate that is derived as a weighted average of all standing woody debris
|
235 | // use the decay rate that is derived as a weighted average of all standing woody debris
|
236 | double survive_rate = exp(-mKSW[i] *climate_factor_re * 1. ); // 1: timestep |
236 | double survive_rate = exp(-mKSW[i] *climate_factor_re * 1. ); // 1: timestep |
237 | mTotalToAtm.C += mSWD[i].C * (1. - survive_rate); |
237 | mTotalToAtm.C += mSWD[i].C * (1. - survive_rate); |
238 | mSWD[i].C *= survive_rate; |
238 | mSWD[i].C *= survive_rate; |
239 | 239 | ||
240 | // transition to downed woody debris
|
240 | // transition to downed woody debris
|
241 | // update: use negative exponential decay, species parameter: half-life
|
241 | // update: use negative exponential decay, species parameter: half-life
|
242 | // modified for the climatic effect on decomposition, i.e. if decomp is slower, snags stand longer and vice versa
|
242 | // modified for the climatic effect on decomposition, i.e. if decomp is slower, snags stand longer and vice versa
|
243 | // this is loosely oriented on Standcarb2 (http://andrewsforest.oregonstate.edu/pubs/webdocs/models/standcarb2.htm),
|
243 | // this is loosely oriented on Standcarb2 (http://andrewsforest.oregonstate.edu/pubs/webdocs/models/standcarb2.htm),
|
244 | // where lag times for cohort transitions are linearly modified with re although here individual good or bad years have
|
244 | // where lag times for cohort transitions are linearly modified with re although here individual good or bad years have
|
245 | // an immediate effect, the average climatic influence should come through (and it is inherently transient)
|
245 | // an immediate effect, the average climatic influence should come through (and it is inherently transient)
|
246 | // note that swd.hl is species-specific, and thus a weighted average over the species in the input (=mortality)
|
246 | // note that swd.hl is species-specific, and thus a weighted average over the species in the input (=mortality)
|
247 | // needs to be calculated, followed by a weighted update of the previous swd.hl.
|
247 | // needs to be calculated, followed by a weighted update of the previous swd.hl.
|
248 | // As weights here we use stem number, as the processes here pertain individual snags
|
248 | // As weights here we use stem number, as the processes here pertain individual snags
|
249 | // calculate the transition probability of SWD to downed dead wood
|
249 | // calculate the transition probability of SWD to downed dead wood
|
250 | 250 | ||
251 | double half_life = mHalfLife[i] / climate_factor_re; |
251 | double half_life = mHalfLife[i] / climate_factor_re; |
252 | double rate = -M_LN2 / half_life; // M_LN2: math. constant |
252 | double rate = -M_LN2 / half_life; // M_LN2: math. constant |
253 | 253 | ||
254 | // higher decay rate for the class with smallest diameters
|
254 | // higher decay rate for the class with smallest diameters
|
255 | if (i==0) |
255 | if (i==0) |
256 | rate*=2.; |
256 | rate*=2.; |
257 | 257 | ||
258 | double transfer = 1. - exp(rate); |
258 | double transfer = 1. - exp(rate); |
259 | 259 | ||
260 | // calculate flow to soil pool...
|
260 | // calculate flow to soil pool...
|
261 | mSWDtoSoil += mSWD[i] * transfer; |
261 | mSWDtoSoil += mSWD[i] * transfer; |
262 | mRefractoryFlux += mSWD[i] * transfer; |
262 | mRefractoryFlux += mSWD[i] * transfer; |
263 | mSWD[i] *= (1.-transfer); // reduce pool |
263 | mSWD[i] *= (1.-transfer); // reduce pool |
264 | // calculate the stem number of remaining snags
|
264 | // calculate the stem number of remaining snags
|
265 | mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer); |
265 | mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer); |
266 | 266 | ||
267 | mTimeSinceDeath[i] += 1.; |
267 | mTimeSinceDeath[i] += 1.; |
268 | // if stems<0.5, empty the whole cohort into DWD, i.e. release the last bit of C and N and clear the stats
|
268 | // if stems<0.5, empty the whole cohort into DWD, i.e. release the last bit of C and N and clear the stats
|
269 | // also, if the Carbon of an average snag is less than 10% of the original average tree
|
269 | // also, if the Carbon of an average snag is less than 10% of the original average tree
|
270 | // (derived from allometries for the three diameter classes), the whole cohort is emptied out to DWD
|
270 | // (derived from allometries for the three diameter classes), the whole cohort is emptied out to DWD
|
271 | if (mNumberOfSnags[i] < 0.5 || mSWD[i].C / mNumberOfSnags[i] < mCarbonThreshold[i]) { |
271 | if (mNumberOfSnags[i] < 0.5 || mSWD[i].C / mNumberOfSnags[i] < mCarbonThreshold[i]) { |
272 | // clear the pool: add the rest to the soil, clear statistics of the pool
|
272 | // clear the pool: add the rest to the soil, clear statistics of the pool
|
273 | mRefractoryFlux += mSWD[i]; |
273 | mRefractoryFlux += mSWD[i]; |
274 | mSWDtoSoil += mSWD[i]; |
274 | mSWDtoSoil += mSWD[i]; |
275 | mSWD[i].clear(); |
275 | mSWD[i].clear(); |
276 | mAvgDbh[i] = 0.; |
276 | mAvgDbh[i] = 0.; |
277 | mAvgHeight[i] = 0.; |
277 | mAvgHeight[i] = 0.; |
278 | mAvgVolume[i] = 0.; |
278 | mAvgVolume[i] = 0.; |
279 | mKSW[i] = 0.; |
279 | mKSW[i] = 0.; |
280 | mCurrentKSW[i] = 0.; |
280 | mCurrentKSW[i] = 0.; |
281 | mHalfLife[i] = 0.; |
281 | mHalfLife[i] = 0.; |
282 | mTimeSinceDeath[i] = 0.; |
282 | mTimeSinceDeath[i] = 0.; |
283 | }
|
283 | }
|
284 | 284 | ||
285 | }
|
285 | }
|
286 | 286 | ||
287 | }
|
287 | }
|
288 | // total carbon in the snag-container on the RU *after* processing is the content of the
|
288 | // total carbon in the snag-container on the RU *after* processing is the content of the
|
289 | // standing woody debris pools + the branches
|
289 | // standing woody debris pools + the branches
|
290 | mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C + |
290 | mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C + |
291 | mOtherWood[0].C + mOtherWood[1].C + mOtherWood[2].C + mOtherWood[3].C + mOtherWood[4].C; |
291 | mOtherWood[0].C + mOtherWood[1].C + mOtherWood[2].C + mOtherWood[3].C + mOtherWood[4].C; |
292 | mTotalSWD = mSWD[0] + mSWD[1] + mSWD[2]; |
292 | mTotalSWD = mSWD[0] + mSWD[1] + mSWD[2]; |
293 | mTotalOther = mOtherWood[0] + mOtherWood[1] + mOtherWood[2] + mOtherWood[3] + mOtherWood[4]; |
293 | mTotalOther = mOtherWood[0] + mOtherWood[1] + mOtherWood[2] + mOtherWood[3] + mOtherWood[4]; |
294 | }
|
294 | }
|
295 | 295 | ||
296 | /// foliage and fineroot litter is transferred during tree growth.
|
296 | /// foliage and fineroot litter is transferred during tree growth.
|
297 | void Snag::addTurnoverLitter(const Species *species, const double litter_foliage, const double litter_fineroot) |
297 | void Snag::addTurnoverLitter(const Species *species, const double litter_foliage, const double litter_fineroot) |
298 | {
|
298 | {
|
299 | mLabileFlux.addBiomass(litter_foliage, species->cnFoliage(), species->snagKyl()); |
299 | mLabileFlux.addBiomass(litter_foliage, species->cnFoliage(), species->snagKyl()); |
300 | mLabileFlux.addBiomass(litter_fineroot, species->cnFineroot(), species->snagKyl()); |
300 | mLabileFlux.addBiomass(litter_fineroot, species->cnFineroot(), species->snagKyl()); |
301 | }
|
301 | }
|
302 | 302 | ||
303 | void Snag::addTurnoverWood(const Species *species, const double woody_biomass) |
303 | void Snag::addTurnoverWood(const Species *species, const double woody_biomass) |
304 | {
|
304 | {
|
305 | mRefractoryFlux.addBiomass(woody_biomass, species->cnWood(), species->snagKyr()); |
305 | mRefractoryFlux.addBiomass(woody_biomass, species->cnWood(), species->snagKyr()); |
306 | }
|
306 | }
|
307 | 307 | ||
308 | /// after the death of the tree the five biomass compartments are processed.
|
308 | /// after the death of the tree the five biomass compartments are processed.
|
309 | void Snag::addMortality(const Tree *tree) |
309 | void Snag::addMortality(const Tree *tree) |
310 | {
|
310 | {
|
311 | const Species *species = tree->species(); |
311 | const Species *species = tree->species(); |
312 | 312 | ||
313 | // immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool
|
313 | // immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool
|
314 | mLabileFlux.addBiomass(tree->biomassFoliage(), species->cnFoliage(), tree->species()->snagKyl()); |
314 | mLabileFlux.addBiomass(tree->biomassFoliage(), species->cnFoliage(), tree->species()->snagKyl()); |
315 | mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl()); |
315 | mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl()); |
316 | 316 | ||
317 | // branches and coarse roots are equally distributed over five years:
|
317 | // branches and coarse roots are equally distributed over five years:
|
318 | double biomass_rest = (tree->biomassBranch()+tree->biomassCoarseRoot()) * 0.2; |
318 | double biomass_rest = (tree->biomassBranch()+tree->biomassCoarseRoot()) * 0.2; |
319 | for (int i=0;i<5; i++) |
319 | for (int i=0;i<5; i++) |
320 | mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr()); |
320 | mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr()); |
321 | 321 | ||
322 | // just for book-keeping: keep track of all inputs into branches / roots / swd
|
322 | // just for book-keeping: keep track of all inputs into branches / roots / swd
|
323 | mTotalIn.addBiomass(tree->biomassBranch() + tree->biomassCoarseRoot() + tree->biomassStem(), species->cnWood()); |
323 | mTotalIn.addBiomass(tree->biomassBranch() + tree->biomassCoarseRoot() + tree->biomassStem(), species->cnWood()); |
324 | // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
|
324 | // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
|
325 | int pi = poolIndex(tree->dbh()); // get right transfer pool |
325 | int pi = poolIndex(tree->dbh()); // get right transfer pool |
326 | 326 | ||
327 | // update statistics - stemnumber-weighted averages
|
327 | // update statistics - stemnumber-weighted averages
|
328 | // note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)
|
328 | // note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)
|
329 | double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers) |
329 | double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers) |
330 | double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1). |
330 | double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1). |
331 | mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new; |
331 | mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new; |
332 | mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new; |
332 | mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new; |
333 | mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new; |
333 | mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new; |
334 | mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new; |
334 | mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new; |
335 | mHalfLife[pi] = mHalfLife[pi]*p_old + species->snagHalflife()* p_new; |
335 | mHalfLife[pi] = mHalfLife[pi]*p_old + species->snagHalflife()* p_new; |
336 | 336 | ||
337 | // average the decay rate (ksw); this is done based on the carbon content
|
337 | // average the decay rate (ksw); this is done based on the carbon content
|
338 | // aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)
|
338 | // aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)
|
339 | if (tree->biomassStem()==0) |
339 | if (tree->biomassStem()==0) |
340 | throw IException("Snag::addMortality: tree without stem biomass!!"); |
340 | throw IException("Snag::addMortality: tree without stem biomass!!"); |
341 | p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction); |
341 | p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction); |
342 | p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction); |
342 | p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction); |
343 | mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + species->snagKsw() * p_new; |
343 | mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + species->snagKsw() * p_new; |
344 | mNumberOfSnags[pi]++; |
344 | mNumberOfSnags[pi]++; |
345 | 345 | ||
346 | // finally add the biomass
|
346 | // finally add the biomass
|
347 | CNPool &to_swd = mToSWD[pi]; |
347 | CNPool &to_swd = mToSWD[pi]; |
348 | to_swd.addBiomass(tree->biomassStem(), species->cnWood(), tree->species()->snagKyr()); |
348 | to_swd.addBiomass(tree->biomassStem(), species->cnWood(), tree->species()->snagKyr()); |
349 | }
|
349 | }
|
350 | 350 | ||
351 | /// add residual biomass of 'tree' after harvesting.
|
351 | /// add residual biomass of 'tree' after harvesting.
|
352 | /// remove_{stem, branch, foliage}_fraction: percentage of biomass compartment that is *removed* by the harvest operation (i.e.: not to stay in the system)
|
352 | /// remove_{stem, branch, foliage}_fraction: percentage of biomass compartment that is *removed* by the harvest operation (i.e.: not to stay in the system)
|
353 | /// records on harvested biomass is collected (mTotalToExtern-pool).
|
353 | /// records on harvested biomass is collected (mTotalToExtern-pool).
|
354 | void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction ) |
354 | void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction ) |
355 | {
|
355 | {
|
356 | const Species *species = tree->species(); |
356 | const Species *species = tree->species(); |
357 | 357 | ||
358 | // immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool
|
358 | // immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool
|
359 | mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), species->cnFoliage(), tree->species()->snagKyl()); |
359 | mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), species->cnFoliage(), tree->species()->snagKyl()); |
360 | mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl()); |
360 | mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl()); |
361 | 361 | ||
362 | // for branches, add all biomass that remains in the forest to the soil
|
362 | // for branches, add all biomass that remains in the forest to the soil
|
363 | mRefractoryFlux.addBiomass(tree->biomassBranch()*(1.-remove_branch_fraction), species->cnWood(), tree->species()->snagKyr()); |
363 | mRefractoryFlux.addBiomass(tree->biomassBranch()*(1.-remove_branch_fraction), species->cnWood(), tree->species()->snagKyr()); |
364 | // the same treatment for stem residuals
|
364 | // the same treatment for stem residuals
|
365 | mRefractoryFlux.addBiomass(tree->biomassStem() * remove_stem_fraction, species->cnWood(), tree->species()->snagKyr()); |
- | |
- | 365 | mRefractoryFlux.addBiomass(tree->biomassStem() * (1. - remove_stem_fraction), species->cnWood(), tree->species()->snagKyr()); |
|
366 | 366 | ||
367 | // split the corase wood biomass into parts (slower decay)
|
367 | // split the corase wood biomass into parts (slower decay)
|
368 | double biomass_rest = (tree->biomassCoarseRoot()) * 0.2; |
368 | double biomass_rest = (tree->biomassCoarseRoot()) * 0.2; |
369 | for (int i=0;i<5; i++) |
369 | for (int i=0;i<5; i++) |
370 | mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr()); |
370 | mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr()); |
371 | 371 | ||
372 | 372 | ||
373 | // for book-keeping...
|
373 | // for book-keeping...
|
374 | mTotalToExtern.addBiomass(tree->biomassFoliage()*remove_foliage_fraction + |
374 | mTotalToExtern.addBiomass(tree->biomassFoliage()*remove_foliage_fraction + |
375 | tree->biomassBranch()*remove_branch_fraction + |
375 | tree->biomassBranch()*remove_branch_fraction + |
376 | tree->biomassStem()*remove_stem_fraction, species->cnWood()); |
376 | tree->biomassStem()*remove_stem_fraction, species->cnWood()); |
377 | }
|
377 | }
|
378 | 378 | ||
379 | // add flow from regeneration layer (dead trees) to soil
|
379 | // add flow from regeneration layer (dead trees) to soil
|
380 | void Snag::addToSoil(const Species *species, const CNPair &woody_pool, const CNPair &litter_pool) |
380 | void Snag::addToSoil(const Species *species, const CNPair &woody_pool, const CNPair &litter_pool) |
381 | {
|
381 | {
|
382 | mLabileFlux.add(litter_pool, species->snagKyl()); |
382 | mLabileFlux.add(litter_pool, species->snagKyl()); |
383 | mRefractoryFlux.add(woody_pool, species->snagKyr()); |
383 | mRefractoryFlux.add(woody_pool, species->snagKyr()); |
384 | }
|
384 | }
|
385 | 385 | ||
386 | 386 | ||
387 | 387 |