Subversion Repositories public iLand

Rev

Rev 605 | Rev 609 | 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"
541 werner 11
#include "model.h"
468 werner 12
 
13
/** @class Snag
14
  Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.
490 werner 15
  Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags
468 werner 16
  is subsequently forwarded to the soil sub model.
522 werner 17
  Carbon is stored in three classes (depending on the size)
528 werner 18
  The Snag dynamics class uses the following species parameter:
19
  cnFoliage, cnFineroot, cnWood, snagHalflife, snagKSW
468 werner 20
 
21
  */
22
// static variables
528 werner 23
double Snag::mDBHLower = -1.;
522 werner 24
double Snag::mDBHHigher = 0.;
25
double Snag::mCarbonThreshold[] = {0., 0., 0.};
26
 
534 werner 27
double CNPair::biomassCFraction = biomassCFraction; // get global from globalsettings.h
468 werner 28
 
534 werner 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)
31
{
32
    if (biomass==0.)
33
        return;
34
    double new_c = biomass*biomassCFraction;
35
    double p_old = C / (new_c + C);
36
    mParameter = mParameter*p_old + parameter_value*(1.-p_old);
37
    CNPair::addBiomass(biomass, CNratio);
38
}
39
 
40
// increase pool (and weigh the value)
41
void CNPool::operator+=(const CNPool &s)
42
{
43
    if (s.C==0.)
44
        return;
45
    mParameter = parameter(s); // calculate weighted parameter
46
    C+=s.C;
47
    N+=s.N;
48
}
49
 
50
double CNPool::parameter(const CNPool &s) const
51
{
52
    if (s.C==0.)
53
        return parameter();
54
    double p_old = C / (s.C + C);
55
    double result =  mParameter*p_old + s.parameter()*(1.-p_old);
56
    return result;
57
}
58
 
59
 
522 werner 60
void Snag::setupThresholds(const double lower, const double upper)
61
{
62
    if (mDBHLower == lower)
63
        return;
64
    mDBHLower = lower;
65
    mDBHHigher = upper;
66
    mCarbonThreshold[0] = lower / 2.;
67
    mCarbonThreshold[1] = lower + (upper - lower)/2.;
68
    mCarbonThreshold[2] = upper + (upper - lower)/2.;
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%
71
    //# values in kg!
72
    for (int i=0;i<3;i++)
73
        mCarbonThreshold[i] = 0.10568*pow(mCarbonThreshold[i],2.4247)*0.5*0.1;
74
}
75
 
76
 
468 werner 77
Snag::Snag()
78
{
490 werner 79
    mRU = 0;
534 werner 80
    CNPair::setCFraction(biomassCFraction);
468 werner 81
}
82
 
490 werner 83
void Snag::setup( const ResourceUnit *ru)
468 werner 84
{
490 werner 85
    mRU = ru;
86
    mClimateFactor = 0.;
468 werner 87
    // branches
88
    mBranchCounter=0;
89
    for (int i=0;i<3;i++) {
90
        mTimeSinceDeath[i] = 0.;
91
        mNumberOfSnags[i] = 0.;
522 werner 92
        mAvgDbh[i] = 0.;
93
        mAvgHeight[i] = 0.;
94
        mAvgVolume[i] = 0.;
95
        mKSW[i] = 0.;
96
        mCurrentKSW[i] = 0.;
557 werner 97
        mHalfLife[i] = 0.;
468 werner 98
    }
475 werner 99
    mTotalSnagCarbon = 0.;
528 werner 100
    if (mDBHLower<=0)
101
        throw IException("Snag::setupThresholds() not called or called with invalid parameters.");
557 werner 102
 
103
    // Inital values from XML file
104
    XmlHelper xml=GlobalSettings::instance()->settings();
574 werner 105
    double kyr = xml.valueDouble("model.site.youngRefractoryDecompRate", -1);
557 werner 106
    // put carbon of snags to the middle size class
107
    xml.setCurrentNode("model.initialization.snags");
108
    mSWD[1].C = xml.valueDouble(".swdC");
109
    mSWD[1].N = mSWD[1].C / xml.valueDouble(".swdCN", 50.);
110
    mSWD[1].setParameter(kyr);
111
    mKSW[1] = xml.valueDouble(".swdDecompRate");
112
    mNumberOfSnags[1] = xml.valueDouble(".swdCount");
113
    mHalfLife[1] = xml.valueDouble(".swdHalfLife");
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 );
116
 
117
    mTotalSnagCarbon = other.C + mSWD[1].C;
118
 
119
    other *= 0.2;
120
    for (int i=0;i<5;i++)
121
        mOtherWood[i] = other;
468 werner 122
}
123
 
475 werner 124
// debug outputs
125
QList<QVariant> Snag::debugList()
126
{
127
    // list columns
128
    // for three pools
129
    QList<QVariant> list;
130
 
523 werner 131
    // totals
132
    list << mTotalSnagCarbon << mTotalIn.C << mTotalToAtm.C << mSWDtoSoil.C << mSWDtoSoil.N;
477 werner 133
    // fluxes to labile soil pool and to refractory soil pool
524 werner 134
    list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N;
475 werner 135
 
136
    for (int i=0;i<3;i++) {
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;
524 werner 139
        list << mAvgDbh[i] << mAvgHeight[i] << mAvgVolume[i];
475 werner 140
    }
141
 
540 werner 142
    // branch/coarse wood pools (5 yrs)
143
    for (int i=0;i<5;i++) {
144
        list << mOtherWood[i].C << mOtherWood[i].N;
145
    }
146
//    list << mOtherWood[mBranchCounter].C << mOtherWood[mBranchCounter].N
147
//            << mOtherWood[(mBranchCounter+1)%5].C << mOtherWood[(mBranchCounter+1)%5].N
148
//            << mOtherWood[(mBranchCounter+2)%5].C << mOtherWood[(mBranchCounter+2)%5].N
149
//            << mOtherWood[(mBranchCounter+3)%5].C << mOtherWood[(mBranchCounter+3)%5].N
150
//            << mOtherWood[(mBranchCounter+4)%5].C << mOtherWood[(mBranchCounter+4)%5].N;
475 werner 151
    return list;
152
}
153
 
468 werner 154
void Snag::newYear()
155
{
156
    for (int i=0;i<3;i++) {
157
        mToSWD[i].clear(); // clear transfer pools to standing-woody-debris
522 werner 158
        mCurrentKSW[i] = 0.;
468 werner 159
    }
160
    mLabileFlux.clear();
161
    mRefractoryFlux.clear();
476 werner 162
    mTotalToAtm.clear();
163
    mTotalToExtern.clear();
164
    mTotalIn.clear();
477 werner 165
    mSWDtoSoil.clear();
468 werner 166
}
167
 
490 werner 168
/// calculate the dynamic climate modifier for decomposition 're'
522 werner 169
/// calculation is done on the level of ResourceUnit because the water content per day is needed.
490 werner 170
double Snag::calculateClimateFactors()
171
{
172
    double ft, fw;
173
    double f_sum = 0.;
552 werner 174
    int iday=0;
553 werner 175
    // calculate the water-factor for each month (see Adair et al 2008)
176
    double fw_month[12];
177
    double ratio;
178
    for (int m=0;m<12;m++) {
562 werner 179
        if (mRU->waterCycle()->referenceEvapotranspiration()[m]>0.)
180
            ratio = mRU->climate()->precipitationMonth()[m] /  mRU->waterCycle()->referenceEvapotranspiration()[m];
553 werner 181
        else
182
            ratio = 0;
183
        fw_month[m] = 1. / (1. + 30.*exp(-8.5*ratio));
564 werner 184
        if (logLevelDebug()) qDebug() <<"month"<< m << "PET" << mRU->waterCycle()->referenceEvapotranspiration()[m] << "prec" <<mRU->climate()->precipitationMonth()[m];
553 werner 185
    }
186
 
552 werner 187
    for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day, ++iday)
490 werner 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)
553 werner 190
        fw = fw_month[day->month-1];
540 werner 191
 
490 werner 192
        f_sum += ft*fw;
193
    }
194
    // the climate factor is defined as the arithmentic annual mean value
195
    mClimateFactor = f_sum / double(mRU->climate()->daysOfYear());
196
    return mClimateFactor;
197
}
198
 
522 werner 199
/// do the yearly calculation
200
/// see http://iland.boku.ac.at/snag+dynamics
526 werner 201
void Snag::calculateYear()
468 werner 202
{
522 werner 203
    mSWDtoSoil.clear();
532 werner 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)
477 werner 205
    if (isEmpty()) // nothing to do
475 werner 206
        return;
207
 
468 werner 208
    // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool
540 werner 209
    mRefractoryFlux+=mOtherWood[mBranchCounter];
210
 
211
    mOtherWood[mBranchCounter].clear();
468 werner 212
    mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0.
540 werner 213
    // decay of branches/coarse roots
214
    for (int i=0;i<5;i++) {
215
        if (mOtherWood[i].C>0.) {
216
            double survive_rate = exp(- climate_factor_re * mOtherWood[i].parameter() ); // parameter: the "kyr" value...
217
            mOtherWood[i].C *= survive_rate;
218
        }
219
    }
468 werner 220
 
221
    // process standing snags.
222
    // the input of the current year is in the mToSWD-Pools
223
    for (int i=0;i<3;i++) {
224
 
522 werner 225
        // update the swd-pool with this years' input
226
        if (!mToSWD[i].isEmpty()) {
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));
229
            //move content to the SWD pool
230
            mSWD[i] += mToSWD[i];
231
        }
475 werner 232
 
522 werner 233
        if (mSWD[i].C > 0) {
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
523 werner 236
            double survive_rate = exp(-mKSW[i] *climate_factor_re * 1. ); // 1: timestep
237
            mTotalToAtm.C += mSWD[i].C * (1. - survive_rate);
238
            mSWD[i].C *= survive_rate;
468 werner 239
 
522 werner 240
            // transition to downed woody debris
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
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
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)
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
249
            // calculate the transition probability of SWD to downed dead wood
468 werner 250
 
522 werner 251
            double half_life = mHalfLife[i] / climate_factor_re;
252
            double rate = -M_LN2 / half_life; // M_LN2: math. constant
253
 
254
            // higher decay rate for the class with smallest diameters
255
            if (i==0)
256
                rate*=2.;
257
 
523 werner 258
            double transfer = 1. - exp(rate);
522 werner 259
 
468 werner 260
            // calculate flow to soil pool...
522 werner 261
            mSWDtoSoil += mSWD[i] * transfer;
262
            mRefractoryFlux += mSWD[i] * transfer;
263
            mSWD[i] *= (1.-transfer); // reduce pool
468 werner 264
            // calculate the stem number of remaining snags
522 werner 265
            mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer);
523 werner 266
 
267
            mTimeSinceDeath[i] += 1.;
522 werner 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
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]) {
272
                // clear the pool: add the rest to the soil, clear statistics of the pool
468 werner 273
                mRefractoryFlux += mSWD[i];
522 werner 274
                mSWDtoSoil += mSWD[i];
468 werner 275
                mSWD[i].clear();
522 werner 276
                mAvgDbh[i] = 0.;
277
                mAvgHeight[i] = 0.;
278
                mAvgVolume[i] = 0.;
279
                mKSW[i] = 0.;
280
                mCurrentKSW[i] = 0.;
281
                mHalfLife[i] = 0.;
282
                mTimeSinceDeath[i] = 0.;
468 werner 283
            }
522 werner 284
 
468 werner 285
        }
522 werner 286
 
468 werner 287
    }
522 werner 288
    // total carbon in the snag-container on the RU *after* processing is the content of the
475 werner 289
    // standing woody debris pools + the branches
290
    mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C +
540 werner 291
                       mOtherWood[0].C + mOtherWood[1].C + mOtherWood[2].C + mOtherWood[3].C + mOtherWood[4].C;
587 werner 292
    mTotalSWD = mSWD[0] + mSWD[1] + mSWD[2];
293
    mTotalOther = mOtherWood[0] + mOtherWood[1] + mOtherWood[2] + mOtherWood[3] + mOtherWood[4];
468 werner 294
}
295
 
296
/// foliage and fineroot litter is transferred during tree growth.
588 werner 297
void Snag::addTurnoverLitter(const Species *species, const double litter_foliage, const double litter_fineroot)
468 werner 298
{
588 werner 299
    mLabileFlux.addBiomass(litter_foliage, species->cnFoliage(), species->snagKyl());
300
    mLabileFlux.addBiomass(litter_fineroot, species->cnFineroot(), species->snagKyl());
468 werner 301
}
302
 
595 werner 303
void Snag::addTurnoverWood(const Species *species, const double woody_biomass)
304
{
305
    mRefractoryFlux.addBiomass(woody_biomass, species->cnWood(), species->snagKyr());
306
}
307
 
468 werner 308
/// after the death of the tree the five biomass compartments are processed.
309
void Snag::addMortality(const Tree *tree)
310
{
528 werner 311
    const Species *species = tree->species();
468 werner 312
 
313
    // immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool
534 werner 314
    mLabileFlux.addBiomass(tree->biomassFoliage(), species->cnFoliage(), tree->species()->snagKyl());
315
    mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());
468 werner 316
 
540 werner 317
    // branches and coarse roots are equally distributed over five years:
318
    double biomass_rest = (tree->biomassBranch()+tree->biomassCoarseRoot()) * 0.2;
468 werner 319
    for (int i=0;i<5; i++)
540 werner 320
        mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());
468 werner 321
 
540 werner 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());
468 werner 324
    // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
522 werner 325
    int pi = poolIndex(tree->dbh()); // get right transfer pool
326
 
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)
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).
331
    mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new;
332
    mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new;
333
    mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new;
334
    mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new;
528 werner 335
    mHalfLife[pi] = mHalfLife[pi]*p_old + species->snagHalflife()* p_new;
522 werner 336
 
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)
339
    if (tree->biomassStem()==0)
340
        throw IException("Snag::addMortality: tree without stem biomass!!");
341
    p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
342
    p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
534 werner 343
    mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + species->snagKsw() * p_new;
522 werner 344
    mNumberOfSnags[pi]++;
523 werner 345
 
346
    // finally add the biomass
534 werner 347
    CNPool &to_swd = mToSWD[pi];
348
    to_swd.addBiomass(tree->biomassStem(), species->cnWood(), tree->species()->snagKyr());
468 werner 349
}
350
 
351
/// add residual biomass of 'tree' after harvesting.
522 werner 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)
528 werner 353
/// records on harvested biomass is collected (mTotalToExtern-pool).
468 werner 354
void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction )
355
{
528 werner 356
    const Species *species = tree->species();
468 werner 357
 
358
    // immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool
534 werner 359
    mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), species->cnFoliage(), tree->species()->snagKyl());
360
    mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());
540 werner 361
 
528 werner 362
    // for branches, add all biomass that remains in the forest to the soil
534 werner 363
    mRefractoryFlux.addBiomass(tree->biomassBranch()*(1.-remove_branch_fraction), species->cnWood(), tree->species()->snagKyr());
528 werner 364
    // the same treatment for stem residuals
605 werner 365
    mRefractoryFlux.addBiomass(tree->biomassStem() * (1. - remove_stem_fraction), species->cnWood(), tree->species()->snagKyr());
468 werner 366
 
540 werner 367
    // split the corase wood biomass into parts (slower decay)
368
    double biomass_rest = (tree->biomassCoarseRoot()) * 0.2;
369
    for (int i=0;i<5; i++)
370
        mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());
371
 
372
 
528 werner 373
    // for book-keeping...
374
    mTotalToExtern.addBiomass(tree->biomassFoliage()*remove_foliage_fraction +
375
                              tree->biomassBranch()*remove_branch_fraction +
376
                              tree->biomassStem()*remove_stem_fraction, species->cnWood());
468 werner 377
}
378
 
588 werner 379
// add flow from regeneration layer (dead trees) to soil
595 werner 380
void Snag::addToSoil(const Species *species, const CNPair &woody_pool, const CNPair &litter_pool)
588 werner 381
{
382
    mLabileFlux.add(litter_pool, species->snagKyl());
383
    mRefractoryFlux.add(woody_pool, species->snagKyr());
384
}
534 werner 385
 
607 werner 386
/// disturbance function: remove the fraction of 'factor' of biomass from the SWD pools; 0: remove nothing, 1: remove all
387
/// biomass removed by this function goes to the atmosphere
388
void Snag::removeCarbon(const double factor)
389
{
390
    // reduce pools of currently standing dead wood and also of pools that are added
391
    // during (previous) management operations of the current year
392
    for (int i=0;i<3;i++) {
393
        mTotalToAtm += (mSWD[i] + mToSWD[i]) * factor;
394
        mSWD[i] *= (1. - factor);
395
        mToSWD[i] *= (1. - factor);
396
    }
534 werner 397
 
607 werner 398
    for (int i=0;i<5;i++) {
399
        mTotalToAtm += mOtherWood[i]*factor;
400
        mOtherWood[i] *= (1. - factor);
401
    }
402
}
403
 
404
 
405
/// cut down swd (and branches) and move to soil pools
406
/// @param factor 0: cut 0%, 1: cut and slash 100% of the wood
407
void Snag::management(const double factor)
408
{
409
    if (factor<0. || factor>1.)
410
        throw IException(QString("Invalid factor in Snag::management: '%1'").arg(factor));
411
    // swd pools
412
    for (int i=0;i<3;i++) {
413
        mSWDtoSoil += mSWD[i] * factor;
414
        mSWD[i] *= (1. - factor);
415
        mSWDtoSoil += mToSWD[i] * factor;
416
        mToSWD[i] *= (1. - factor);
417
    }
418
    // what to do with the branches: now move also all wood to soil (note: this is note
419
    // very good w.r.t the coarse roots...
420
    for (int i=0;i<5;i++) {
421
        mRefractoryFlux+=mOtherWood[i]*factor;
422
        mOtherWood[i]*=(1. - factor);
423
    }
424
 
425
}
426
 
427