Predicting constant decay rates of coarse woody debris—A meta-analysis approach with a mixed model

Ecological Modelling 220 (2009) 904–912

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Predicting constant decay rates of coarse woody debris—A meta-analysis approach with a mixed model Jürgen Zell ∗ , Gerald Kändler, Marc Hanewinkel Forest Research Institute of Baden-Württemberg, Wonnhaldestr. 4, 79100 Freiburg, Germany

a r t i c l e

i n f o

Article history: Received 29 July 2008 Received in revised form 21 January 2009 Accepted 25 January 2009 Available online 25 February 2009 Keywords: Coarse woody debris Decomposition Deadwood Mixed model Meta-analysis

a b s t r a c t The aim of the study is the estimation of decay rates for coarse woody debris in large forest regions. These rates, together with estimations of the amount of deadwood, can be used to calculate the release of carbon from that pool into the atmosphere. The model can be used for predictions of decomposition rate constants in a wide range of forest areas (e.g. in process based ecological models, reporting of GHG-emissions), as only easily available predictor variables were used in the regression. Based on an intensive literature research a meta-analysis on influencing factors controlling the constant decay rate of coarse woody debris was set up. The included studies differed significantly in the survey methods as well as in the geographical origin. 39 studies were collected, 30 appeared in North America and nine in Europe. Based on these studies 291 observations of the remaining fraction of coarse woody debris were collected. To quantify the effects that influence the decomposition rates a nonlinear mixed effects model was constructed. Only physiologically interpretable variables were included. With this approach it was possible to determine influencing effects from mean temperature in July, annual rainfall (as quadratic term), diameter of woody material and grouping into hardwoods or conifers and mass- or density loss were significant variables. The mixed effects model also allowed an estimation of the species-specific effects on the decomposition process. These random effects are given for 42 tree species. The degrees of freedom were used efficiently. The model explains 79.6% of the variance and is superior to a comparable multiple regression model. © 2009 Elsevier B.V. All rights reserved.

1. Introduction and background Deadwood in forests is considered to be an important component of the ecosystem. It serves as habitat for many plant-, animaland fungi-species. In Germany, 25% of approximately 5200 beetle species rely on deadwood (Albrecht, 1991). Especially natural forests are distinguished by deadwood as an important element, thus representing an important storage of carbon and nutrients. With regard to climate change and the reporting needs established by the United Nations Framework Convention on Climate Change (UNFCCC), the amount of stored carbon warrants inventory. Instead of measuring the output, the pool can be multiplied by a constant decay rate (Good-Practice-Guidance, IPCC, 2003). Within Europe there are only a few studies with measured decay rates of coarse woody debris. Obviously there is a lack of information between the needs for the reporting and given studies,

∗ Corresponding author at: Forest Research Institute of Baden-Württemberg, Biometrics and Informatic, Wonnhaldestr. 4, 79100 Freiburg, Germany. Tel.: +49 761 4018 187; fax: +49 761 4018 355. E-mail address: [email protected] (J. Zell). 0304-3800/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2009.01.020

which report on decay rates for very specific situations. Measurements of decay rates are time consuming and the estimated decay rates are only valid for a given situation. Usually the studies come up with some influencing factors, which have a more or less significant influence on the decomposition process, but as there are usually only few observations, other maybe stronger relationships, will not be found. We tried to bring these factors together into unified mathematical background and estimate the effects based on observations from literature. The timely dynamic of the deadwood pools is determined by income and losses (Olson, 1963). Income is driven on the one hand through mortality of trees, either natural or by human intervention, and on the other hand by the drop-off of dead branches. Losses are driven by the decay of deadwood. The decay is a complex process, mainly driven by heterotrophic respiration of decomposers (Mackensen et al., 2003) leading to a direct release of CO2 . A systematic assessment of deadwood is usually assured by National Forest Inventories. However, only deadwood above a certain threshold is usually being assessed. E.g. the calculated mean volume of deadwood (>20 cm diameter) in the forests of Germany is 11.5 m3 /ha (BMELV, 2007). This rather small stock of deadwood

J. Zell et al. / Ecological Modelling 220 (2009) 904–912

is the result of a forest management that regularly harvests trees before they die. In order to describe the development of deadwood over time, harvesting interventions in the ecosystem have to be quantified and knowledge about the decay of deadwood is indispensable. In Germany, only two studies dealing with the decay of deadwood are available (Kahl, 2003; Müller-Using, 2005) whilst a lot of studies have been published in North America (see the reference lists in Harmon et al., 1986; Yin, 1999). Therefore, a comprehensive meta-analysis concerning this topic was compiled in the present investigation using the accessible studies reverted to literature to develop a decay model that can be applied to many tree species in a variety of climatic conditions. This statistical model could found the basis of climate- and species-specific constant decay rates that can serve to calculate steady states for deadwood and for predictions of deadwood pools in ecological models. Further it can be used for an advanced quantification of constant decay rates, as they are needed for the accounting of CO2 release within the reporting needs given by the IPCC (IPCC, 2003). Although it is known, that the process of decay is not constant over time, a constant decay rate was modeled. There are several

905

approaches describing the process of decomposition. Most often the single exponential model (Olson, 1963) is used. If only a short time span is observed, a linear decomposition model can also be estimated. Laiho and Prescott (1999) used a sigmoidal function to fit mass loss. Yatskov et al. (2003) divided the decay-process into phases, which itself can be described by different simple exponential models. Further it is possible to model the process based on physiological assumptions. E.g. Yin (1999) used a process described in Agren and Bosatta (1998), where the decay process is driven by a changing quality of the substrate. If a functional representation of the process is not needed, matrix models (Kruys et al., 2002) can also be applied. In order to refer as many studies as possible, mass loss as well as density loss based measurements where taken into account. Harmon et al. (2000) concluded that decomposition rates, based on density loss are lower, than those based on mass loss. This can be expected, since fragmentation is not incorporated in density losses. Furthermore studies deriving the remaining deadwood fraction using population monitoring have been included as well as studies relying on chronosequences. Population monitoring assesses

Table 1 Overview of the study sites and climatic conditions. Location

n

Elevation (m)

Mean temperature (◦ C) Ann.

Alberta 1 Alberta 2 Alberta 3

Precipitation (mm)

January

July

Ann.

−12.8 −11 −10.3

11.5 13.3 14

Air hum. (%)

P/EP (%)

July

July

60 60 60

71 70 69

90 81 78

January

July

660 660 660

35 35 35

Author

Taylor et al. (1991)

2 4 1

1830 1500 1410

0.2 2 2.7

California Germany 1 Germany 2 England 1 England 2 Indiana Lousiana Minnesota 1 Minnesota 2 Nevada New Brunswick

4 8 55 7 4 13 2 8 3 1 6

2250 495 765 80 45 200 9 400 360 2020 10

10 6.4 6.6 9 12 12.1 19.5 3.1 2.4 6.2 2

−1 0 −1.39 1 4.8 −1.6 10.4 −17 −17.7 −3.2 −13

18.8 16.3 15.12 14 19 24.4 27.6 19.6 18.9 17.2 18

1130 1035 1944 670 1144 1000 1420 639 620 232 900

271 57 284 80 258 76 131 21 23 43 84

8 73 64 20 72 106 132 99 98 7 83

47 n.a. n.a. n.a. n.a. 74 78 74 75 50 78

9 n.a. n.a. n.a. n.a. 88 99 98 99 9 89

New Hampshire 1 New Hampshire 2 New Hampshire 3 New Mexico Netherlands North Carolina

13 24 8 1 1 18

530 700 1235 3170 15 850

4.3 3.4 0 0 9 12.6

−10.2 −10 −12 −10 1 2.6

17.4 16.4 12 10 17 21.5

970 1340 1500 700 770 1820

72 112 104 32 70 178

80 118 143 135 70 129

74 76 79 77 n.a. 74

92 142 220 247 n.a. 126

Norway Oregon 1 Oregon 2 Oregon 3 Oregon 4 Oregon 5 Quebec Russia

5 5 2 14 5 12 2 12

370 1670 1070 200 460 850 522 240

5.7 4.4 8.5 10.1 6.5 7.8 −4.8 2.6

−6.8 −4.7 2.3 4.6 0.3 1 −22 −7.5

15.8 14.6 20.6 15.7 14.6 15.5 12 12.5

755 280 2200 3420 2200 2000 785 775

49 47 330 560 330 300 48 53

84 9 20 46 20 30 99 89

n.a. 57 54 76 72 68 79 n.a.

n.a. 12 19 56 27 38 139 n.a.

South Carolina Tennessee Utah Wales Washington 1 Washington 2

8 4 1 9 12 8

250 180 2400 335 210 450

15.9 15.1 2.6 8 9.4 9.7

5.5 2.8 −6 2 2.7 2

25.8 26.3 15 12 16.7 12

1380 1230 1140 2800 1440 1020

133 114 100 278 201 150

114 97 54 167 38 20

74 73 64 n.a. 74 75

91 75 71 n.a. 44 30

Washington 3 Washington 4 Washington 5 Washington 6 Washington 7

8 7 8 7 8

1130 430 60 750 111

4.7 9.2 10 8.3 9.7

−2.7 2.9 4.7 −1.9 4.6

11.8 17.2 15.5 19.6 14.7

2300 1040 2500 600 3300

340 150 400 104 501

54 20 39 10 60

74 68 77 56 78

81 23 47 10 77

Erickson et al. (1985)

Washington 8 Wyoming

8 5

1150 2800

5.4 2.8

−3.2 −8.5

14.4 12

2730 600

431 150

46 50

69 67

60 79

Edmonds (1987) Fahey (1983)

n.a.: not available.

Müller-Using (2005) Kahl (2003) Boddy and Swift (1984) Swift et al. (1976) Rice et al. (1997) Miller (1983) Yin (1999) Gosz et al. (1973) and De Vries and Kuyper (1988) Fahey Abbott and Crossley (1982) et and Mattson et al. (1987) al. Naesset (1999) (1988) Fogel and Cromack (1976)

Moore (1984) Shorohova and Shorohov (2001) Barber and van Lear (1984) Gosz (1980) Fahey et al. (1991) Edmonds (1987) Edmonds et al. (1986)

906

J. Zell et al. / Ecological Modelling 220 (2009) 904–912

different decay processes at the same point in time. Such an investigation bears the risk of not assessing the already disappeared pieces of deadwood that showed the highest decay rates. The decay rate may therefore be overestimated (Yin, 1999). 2. Database Altogether 39 studies have been analyzed (see Table 1), in order to develop a decay model based on published data. The NorthAmerican studies have already been used by Yin (1999) to develop a decay model specific to North America. The theoretical background of his model is based on a description of microbial decay processes by Agren and Bosatta (1998). Yin (1999) developed a model, which has a time-dependent decomposition process. This model is only applicable in North-America, as it is depending on the latitude, which is only a reasonable variable within this region. The dataset from Yin (1999) was adopted and the respective studies were analyzed, in order to check and enlarge the data by other variables. Additionally, the dataset was completed by nine European studies. Table 1 contains the climatic situation in the different study areas as well as the references. As the decay process is influenced by climatic conditions, the mean annual temperature as well as the mean temperature in July and January and the sum of precipitations were taken from these studies as potential explaining variables. Unlike in Yin (1999), in the present study a model was developed that uses the exponential decay as a basis. This results in a constant decay rate, for which climate and species-specific effects where analyzed. Mean remaining deadwood fraction of a tree species in a study area was used as dependent variable. If possible, the deadwood fraction was sub-divided into branch and stem-parts or into distinct different diameter classes. One study can contain several observations. The time to reach the remaining deadwood fraction was considered to be the independent variable. In addition, other co-variables at tree-level were derived, either directly from the studies (mean diameter of the pieces of deadwood), or using the existing literature (e.g. species-specific mass densities, Kollmann, 1982). In addition, an indicator variable was used, to separate coniferous from hardwood (ICH ). Table 5 contains an overview of the mean values of the species-specific variables. 3. Methods Usually studies on decomposition only report a mean constant decay rate, which is valid for the given observations under the specific conditions. Using a large number of observations it is sometimes possible to extract variables influencing this mean (e.g. Naesset, 1999; Tarasov and Birdsey, 2001). These statistical correlations are needed for process based ecological models as they are used for predictions in unobserved situations. For example, the process based decomposition model “Yasso” (Liski et al., 2005; Palosuo et al., 2005) uses a regression on mass loss data derived in the Leningrad Region in Russia (Tarasov and Birdsey, 2001). In this case decomposition rate constants were depending on summer drought and temperature. Even if this relation is physiological justified, it is applied (extrapolated) in the model “Yasso” outside Russia (Kaipainen et al., 2004; Thürig et al., 2005; Böttcher et al., 2008). It is therefore of high interest to bring different studies on decomposition together and perform an overall regression model, which is based on a broad data-basis. In this paper we tried to use as much as possible observations and extract physiological sound variables influencing the mean decomposition rate constant. These variables should be easily available, so that predictions on decomposition rate constants can be carried out straight forward. In the

following we give an overview of different equations describing the decay process. Studies concerning deadwood decay usually imply derived decay rates. If the period of decay is not too long, the decay can be described by a simple linear function. However, for longer periods distinct deviations from the linear development were found. The most common approach to model decay is that of a first order differential equation (Olson, 1963): dx = −kx(t) dt

(1)

This approach assures a decay rate k, constant over time. The quantity of decomposed wood is the product of this decay rate and the amount of wood at time t. The solution of this differential equation is the exponential function xt = x0 exp(−kt), with a constant x0 , that can be interpreted as the initial value. If there is more than one observation, the equation can be expressed in logarithmic form, so k and x0 can be calculated using least-square estimates. The decay rate k itself can depend on other co-variables which can be expressed as (k = X ˇ + ε) in a multiple linear regression. By taking the initial value x0 to the left side of the equation, we obtain xt /x0 = Rt as a measure of the relative fraction of deadwood at time t (Naesset, 1999). Rt is labeled remaining fraction, as mass and density losses where both used in the following equations: ln R(t) = ln ˇ0 − (X  ˇ)t + ε

(2)

with ln(ˇ0 ), a constant indicating the deviation from the relative starting value at time t = 0, which should be normally 1 or 100%. Estimating different decay rates for the phases of deadwood decomposition implies the problem of separating the different phases from each other. A separation into phases with the given data is impossible, as data from different studies were used. There are different definitions for the degree of decomposition of the deadwood on the one hand, and on the other hand, passing times for different phases of decomposition are rarely given. This problem can be solved by simultaneously estimating two (or more) decay rates, if the limits between the different phases are not known. For example, double exponential models are used to describe litter-decomposition of components with different durability (e.g. Berg and McClaugherty, 2003): R(t) = ˇ0 e(−k1 t) + (100% − ˇ0 )e(−k2 t)

(3)

This equation can lead to a successful nonlinear estimation of parameters, but is highly instable as parameter-estimations for k1 and k2 will considerably change if datasets are removed or covariables will be added to the exponent (further discussion with nonlinear-regression models see: Ratkowsky, 1990). With the given data, a model with three decay rates leads to difficulties in convergence. These problems in modeling were one reason for using a simple exponential decay in the present study. Further with a non-constant decay rate, one would have to simulate steady state of the deadwood pool as a time-dependent process. For many approaches such data are lacking or could only be estimated by simulation of standdevelopments. Only for reasons of comparison, such models where fit, to see if they are statistically better. One approach with a sigmoidal function (according to Laiho and Prescott, 1999) was tested. R(t) = 100% e−(ˇ0 ·t)

ˇ1

(4)

In the following the simple exponential model was further refined by finding co-variables (beside the time) which influence the constant decay rate. The constant decay rate should always be negative, for all possible predictions. Therefore a transformation of the rate −k was carried out (−k = −exp(X ˇ)). Further there are some observations with time close to zero, which have already lost some amount of mass or density. If we do not account for this, we would

J. Zell et al. / Ecological Modelling 220 (2009) 904–912

907

Fig. 1. Shows fitted values and corresponding residuals of the first two simple decay models. Models are expressed in Table 2.

get skewed residuals. This can be avoided by introducing a variable ˇ0 , which gives the mean amount of remaining at time = 0 (which is less than 100%). Another possibility to avoid this skewness is by introducing an extra error to the measurement of time, a “time lag” (lag). The latter is easier to interpret and can be seen as the amount of time the trees were already in the process of being decomposed, but were still standing. It can be regarded as the mean measurement error, according to the start of the decomposition process. This variable leads to better statistical results, compared to the variable ˇ0 . With increasing time of decomposition, the variance of the observed remaining fraction will obviously increase. So a weighting

according to the increasing variance was estimated based on the ˆ Finally the nonlinear multiple regression model showed data (ı). systematic deviations of the residuals concerning the different tree species. To account for these differences a nonlinear multiple regression model was further enlarged with a random effect (bi ) for the different species (S) (mixed model, see Pinheiro and Bates, 2000). Index i comprises the ith tree species and j the jth observation within a tree species. (ˇ0 +bi +X  ˇ)

Rij = 100% − 100% e(−kt) = 100% − 100% e[−e

·(t+lag)]

+ εij (5)

With an increasing variance, as decomposition proceeds: ˆ

var(εij ) =  2 Rˆ ij2ı And two random variables: bˆ i ∼ N(0, S2 )

εˆ ij ∼ N(0,  2 )

4. Results 4.1. Simple decay models

Fig. 2. Observations (n = 287) of remaining fraction of coarse woody debris, derived from 39 references, with four simple nonlinear-regression models (models see Table 2).

At first, an exponential model with an initial value of 100% of the remaining fraction of deadwood Rt and a fixed parameter k, representing the decay rate was developed based on nonlinear regression (model 1). As the initial value is forced to be at 100% of the remaining deadwood fraction, the model showed a tendency in the residuals over the fitted values. The second model avoids this by using a time-lag variable (R(t) = 100% e−k·(t+lag) , see Fig. 1). Thirdly a simultaneous estimation of two decay rates (see Eq. (3)) was carried out. And fourthly a time-dependent process with a sigmoidal function was fitted (see Eq. (4)).

Table 2 Parameters for four different simple regression models for coarse woody debris. No.

Model

Parameter

est

sep

t-Value

p-Value

sem

Nonl. coef. det. (%)

1

R(t) = 100% e−k·t

k

0.031

0.002

17.45

0.0000

20.04

13.9

2

R(t) = 100% e−k·(t+lag)

k lag

0.017 −11.36

0.002 1.937

10.54 5.88

0.0000 0.0000

17.00

38.2

ˇ0 k1 k2

79.51 0.015 1.072

0.002 0.414 2.240

8.20 2.60 35.49

0.0000 0.0101 0.0000

16.87

39.4

ˇ0 ˇ1

0.015 0.457

0.002 0.036

6.21 12.78

0.0000 0.0000

16.69

40.5

3

4

R(t) = ˇ0 e−k1 ·t + (100% − ˇ0 )e−k2 ·t R(t) = 100% e(ˇ0 ·t)ˇ1

R(t): remaining fraction at time t. est: estimated value; sep: standard error parameter; sem: standard error model.

908

J. Zell et al. / Ecological Modelling 220 (2009) 904–912

Table 3 Parameter estimates of the mixed nonlinear decay model. Parameter

Variable

Unit

Estimate

Standard error

t-Value

p-Value

ˇ0 ˇ1 ˇ2 ˇ3 ˇ4 ˇ5 ˇ6 ˇ6

intercept ICH d IMD tj py py lag

– – cm – ◦ C mm mm years

−4.07 −0.467 −1.07E−02 0.625 4.71E−02 −5.92E−04 −2.41E−07 −1.32

0.305 0.123 2.57E−03 9.04E−02 4.71E−02 5.92E−04 6.01E−08 0.359

−13.31 −3.79 −4.18 6.91 3.91 2.61 −4.02 −3.68

0.0000 0.0002 0.0000 0.0000 0.0001 0.0096 0.0001 0.0003

ICH : indicator-variable that takes the value 1 if the tree species is coniferous otherwise 0, d: mean diameter (cm), IMD : indicator-variable that takes the value 1 if the study has observed loss of mass, otherwise 0 (loss of density), tj : mean temperature in July (◦ C), py : sum of precipitation per year (mm), lag: starting time of decomposition process in years, before t = 0.

Table 2 presents summary statistics for the four models. It can be seen, that the incorporation of a time lag improves the fit of model 1 substantially and the residuals become more homogenous over the range of fitted values (see Fig. 1). On the other hand, the estimate of a time lag has a dramatic influence on the mean decay rate: it falls from 0.031 to 0.017, as it assumes, that the decay process has

already started 11.36 years before the given starting point in the studies. As with further improvements of the model the effect of the time lag decreases, it was decided to follow this approach due to its favorable statistics. Furthermore models 3 and 4 clearly indicate that the process of decomposition can be better approximated with timely not constant decay rates, either by dividing into artificial

Table 4 Overview of species-specific variables (mean values) and estimated best linear unbiased predictor (BLUP, random effects). n denotes number of observations derived from literature. Tree species

n

Diameter (cm)

Density (kg dm−3 )

Time (years)

Remaining fraction (%)

BLUP

Abies amabilis Abies balsamea Abies concolor Abies lasiocarpa Acer rubrum Acer saccarum Acer saccharinum Alnus rubra Betula alleghaniensis Betula papyrifera Betula pendula Carya spp. Cornus florida Corylus colurna Fagus grandifolia Fagus sylvatica Fraxinus exelsior Fraxinus profunda Juniperus communis Liriodendrum tulipifera Nyssa sylvatica Oxydendrum arboreum Picea abies Picea engelmannii Picea glauca Picea mariana Picea rubens Picea sitchensis Pinus banksiana Pinus contorta Pinus jeffreyi Pinus ponderosa Pinus resinosa Pinus rigida Pinus sylvestris Pinus taeda Populus tremoloides Populus tremula Prunus pensylvanica Pseudotsuga menziesii Quercus alba Quercus coccinea Quercus prinus Quercus robur Robinia pseudoacacia Thuja spec. Tsuga canadensis Tsuga heterophylla

14 14 4 2 2 10 1 13 4 1 6 4 4 1 7 20 1 2 1 1 1 1 11 2 3 2 20 13 4 12 1 6 2 1 2 8 8 5 1 26 5 1 4 2 1 4 2 31

5.59 13.25 20.00 1.38 7.75 12.60 5.00 4.08 8.25 0.50 6.33 28.50 4.08 4.00 12.86 13.79 4.00 1.50 7.50 6.40 11.00 13.00 9.65 1.38 8.83 0.60 20.43 20.53 7.50 14.65 1.80 5.75 14.50 6.90 8.75 3.13 7.77 7.80 0.50 28.77 39.40 15.00 4.30 5.75 9.70 37.00 14.80 25.03

0.48 0.41 0.50 0.35 0.62 0.74 0.51 0.47 0.75 0.66 0.61 0.76 1.00 0.58 0.77 0.68 0.65 0.59 0.55 0.43 0.55 0.59 0.43 0.35 0.52 0.48 0.46 0.42 0.51 0.49 0.47 0.42 0.51 0.54 0.49 0.62 0.43 0.43 0.56 0.47 0.71 0.71 0.79 0.65 0.73 0.35 0.43 0.50

3.9 27.9 17.1 4.0 4.0 6.9 3.5 1.9 4.8 2.0 16.0 9.5 5.3 2.0 7.3 5.0 2.0 2.5 8.0 6.0 6.0 6.0 20.1 4.0 11.7 2.0 19.9 14.7 8.0 17.0 2.0 11.0 15.5 6.0 8.0 2.8 6.8 14.3 2.0 22.2 12.2 6.0 3.8 5.0 6.0 34.3 43.5 13.7

90.5 59.4 64.3 82.4 64.6 63.9 73.2 78.9 59.8 80.0 46.5 54.7 59.4 51.8 57.3 62.3 69.7 45.3 76.7 55.8 59.4 86.3 62.5 84.8 47.1 69.5 56.5 76.0 66.0 69.9 87.1 85.9 45.5 75.0 51.8 85.1 65.0 48.8 89.0 76.8 65.0 67.9 73.1 49.4 92.3 76.5 59.1 82.9

−0.007 −0.329 0.081 −0.083 0.102 −0.158 −0.296 −0.009 0.101 −0.134 0.016 0.117 0.184 0.274 0.036 0.173 0.003 0.150 −0.253 0.142 0.109 −0.272 −0.082 −0.169 0.409 0.372 0.214 0.634 0.159 0.005 −0.030 −0.647 0.277 0.025 0.106 −0.083 0.027 0.185 −0.257 −0.264 −0.251 −0.001 0.257 0.123 −0.365 −0.517 −0.209 0.136

J. Zell et al. / Ecological Modelling 220 (2009) 904–912

909

Fig. 3. Predicted mass loss decay rates with the nonlinear mixed effects model for Fagus Sylvatica and Norway Spruce. Gray dots in the background: precipitation and temperature from studies. With increasing diameter, decay rate declines, temperature has a clear positive effect on decay rate and annual rainfall shows an optimum in between 1100 and 1300 mm.

phases (model 3) or directly by applying a sigmoidal curve (model 4). The sigmoidal model, with one parameter less, explains more variance than model 3, so from a statistical point of view, this would be the most promising type of model. Fig. 2 shows these four models with the underlying data. It can be recognized that the bending of the curve with a fixed initial value (100%) is much more distinct than with a flexible time lag (model 2), which reveals a significant over-estimation of the decay rate after 30–40 years. Some observations in Fig. 2 show that no decay at all takes place in the “slow” phase (according to nomenclature by Yatskov et al., 2003). This, in turn, suggests the use of a sigmoidal function. Such a model, inevitably leads to non-constant decay rate in time. As visible in Fig. 2 there are distinct differences in the decay between coniferous and hardwood species. 4.2. Mixed nonlinear decay model As already stated in Eq. (5), it was tried to explain the constant decay rate by other variables. These variables where not found by applying a stepwise regression, rather they were implemented in the model by hand. The reason for this is the plausibility of the resulting model. It should be parsimonious and only contain variables, which have a physiological explanation. With the help of a residual analysis (plotting residuals against new variables) one can find variables that promise to improve the model. After adding each new variable a likelihood ratio tests (based on a 2 -distribution) was performed, to check whether the new model generates a significant improvement. The best model, found that way, has the following equation: (ˇ0 +bi +ˇ1 I CH +ˇ2 d+ˇ3 I MD +ˇ4 tj +ˇ5 py +ˇ6 p2 ) y ·(t +lag)] ij

Rij = 100%−100% · e[−e

+εij (6)

With an increasing variance, as decomposition proceeds: ˆ

var(εij ) =  2 Rˆ ij2ı And two random variables: bi ∼ N(0, S2 )

εij ∼ N(0,  2 )

With the variables: • ICH : indicator-variable that takes the value 1 if the tree species is coniferous, otherwise hardwood it takes the value 0 • d: mean diameter (cm) • IMD : indicator-variable that takes the value 1 if the study has observed loss of mass, otherwise 0 (loss of density) • tj : mean temperature in July (◦ C) • py : sum of precipitation per year (mm) • lag: starting time of decomposition process in years, before t = 0 • S: tree species (random effect) • ı: factor describing the increasing variance This model explains 79.7% of the variance and is superior to an equivalent multiple regressions model, with same variables. Table 3 contains all parameters associated with the model. The mixed nonlinear model additionally includes the random parameter bi added to the intercept ˇ0 in the exponent, to account for differences between tree species. The dataset contains 48 tree species with respectively 1–31 observations. The estimated random effects are given in Table 4. The variance between the tree species (S2 ) and the variance increasing factor (ı), differ significantly from zero. Therefore, tree species should be taken into account when applying the model, although there is a relatively high residual variance ( 2 = 43.21).

910

J. Zell et al. / Ecological Modelling 220 (2009) 904–912

Table 5 Comparison between some predicted decay rates and observed decay rates in literature. Density (D) or mass loss (M)

Diameter (mean) (cm)

Temperature July (◦ C)

Ann. precipitation (mm)

Included in data?

D

10 20 47.5 12.5 30 50 10 37. 5 15 42.5

10

650

y

Tarasov and Birdsey (2001)

Picea abies Picea abies Picea abies Picea abies Picea abies

D

11.2 11.6 11.9 11.7 12

15.8

760

y

Naesset (1999)

0.027

Picea abies

D

11.3

y

Kahl (2003)

0.041–0.079 0.015–0.038 0.027–0.052

Pinus contorta Picea glauca Abies lasiocarpa

M

20

n

Herrmann and Prescott (2008)

Predicted value

k, literature

Tree species

0.043 0.038 0.029 0.018 0.015 0.012 0.033 0.025 0.037 0.028

0.024–0.058 0.010–0.027 0.017–0.019 0.059 0.005–0.039 0.002–0.024 0.088 0.039 0.071 0.044

Pinus sylvestris Pinus sylvestris Pinus sylvestris Picea abies Picea abies Picea abies Betula pendula Betula pendula Populus tremula Populus tremula

0.025 0.025 0.025 0.025 0.025

0.030 0.035 0.035 0.036 0.029

0.023 0.042 0.062 0.038

15.1

14

1944

638

Study

Ranges of k in Herrmann and Prescott (2008) and Tarasov and Birdsey (2001) due to different phases of decomposition.

The model gives directly conclusions on the expected decomposition rate constant: (ˇ0 +bi +ˇ1 I CH +ˇ2 d+ˇ3 I MD +ˇ4 tj +ˇ5 py +ˇ6 p2 y)

k(S, CH, MD, tj , py ) = e

Some authors suggest calculating the time when 95% of the original mass would have been decomposed (often referred to as turnover time). It can be calculated as t(S, CH, MD, tj , py )0.95 =

3 (ˇ0 +bi +ˇ1 I CH +ˇ2 d+ˇ3 I MD +ˇ4 tj +ˇ5 py +ˇ6 p2 y)

e

Within this model, conifers decompose slower by a factor of exp(ˇ1 ) = 0.627. Or a 10 cm thicker log slows decomposition by exp(ˇ2 ·10) = 0.898. As expected, studies based on mass loss (vs. density loss) show a clearly higher decomposition rate (exp(ˇ3 ) = 1.87). Increasing temperature has an accelerating effect on decomposition, often referred to value called Q10 , which is just the increasing factor of elevating temperature by 10 ◦ C. Q10 is then estimated to be 1.60. The maximum of decomposition concerning annual rainfall is found by setting the first partial derivate of the above equation to zero. It can found by a value of (−ˇ5 /2ˇ6 ) = 1226. No significant differences where found between studies relying on population monitoring compared to studies relying on chronosequences. Nor mass density, a variable associated on the species-level, was found to be significant (for both multiple regression model, as well as mixed model). As the main interest in this model is the prediction of decay rates, in Fig. 3 there are examples to predict mass loss decay rates for Norway Spruce and Fagus Sylvatica, two European tree species, where only few observations exists. Diameter of the piece of deadwood and climatic variables, expressed in mean temperature in July and annual rainfall, are main drivers influencing the decay rate, for a given tree species. The model contains a quadratic relationship between mean annual rainfall and the decay rate, which implies an optimum decay rate according to annual rainfall. Table 5 contains a few examples of comparisons between predicted decomposition rate constants by this model and some studies, in which these rates where estimated. As Tarasov and Birdsey (2001), as well as Herrmann and Prescott (2008) divide phases with different decomposition rates, ranges of k are presented. For Picea abies the estimated values fits well to the values given in the studies, with only one exception. The declining trend

of decomposition rates with increasing diameter can also be found in the measurements presented from Tarasov and Birdsey (2001). As the study from Herrmann and Prescott (2008) is not included in this database, it is remarkably, how well the estimations fit the observed ranges of decomposition rates. Only Picea glauca is purely estimated. The deviation could be due to a low number of observations (n = 3) for this tree species. 5. Discussion In the present model, a constant decay rate over time was assumed. A clearly more flexible description of the overall decay process can be derived assuming a time lag. This means, that in average the process of decay has already started 1.3 years before starting time given in the studies. Maybe the estimated starting time of the decomposition process is erroneous, with a tendency to underestimate the time span. This is reasonable, as dating logs in advanced stages of decomposition is difficult (Laiho and Prescott, 2004). With the new parameter the skewness in the residuals disappears, they loose their tendency over the predicted values. Cleary the given data would recommend to use of phases for the decay processes (Yatskov, 2003) or to model a non-constant decay process (Yin, 1999). As the specification of the biological process of coarse woody debris is not in the aim of this study, rather finding driving variables that influence decomposition rate constant on a broad basis, this assumption can still be maintained. Yin (1999) formulated a decay process that is not constant over time (with a comparable dataset). This approach was not applied in the present paper, for several reasons. Firstly, the point in time of the input has to be known, if a non-constant decay process is used in further applications. This is seldom the case for large-scale estimates of coarse woody debris. Further the assumption of a constant decay rate will promise a simpler application due to its possibility to calculate steady states and its easier handling in modeling of carbon stocks in larger areas, where a point in time of input can only be roughly estimated. The difficult modeling process was the other reason not to use a constant decay rate. With the present dataset it was not possible to estimate a triple exponential model without difficulties of convergence. Multiple regression analysis leads to a good explanation of the constant decay rate over time by other co-variables. Only variables were selected, where a relation to the decay process was

J. Zell et al. / Ecological Modelling 220 (2009) 904–912

physiologically justifiable. One group of these variables comprises climatic variables, such as mean temperature in July and sum of precipitation per year (Mackensen et al., 2003). The other group encompasses tree species-related variables such as conifers or deciduous tree species and the mean diameter of the pieces of deadwood (Harmon et al., 1986) that showed a significant effect on the decay rate. As these variables are easily available, the model can be further used on a broad basis for predictions in ecological models and is advantageous compared to statistical models, which are not applicable in the region of interest (e.g. Yin, 1999), or models which rely on derived decomposition rate constants of a given region or study (Kram et al., 1999; Tinker and Knight, 2001; Verburg and Johnson, 2001; Scheller and Mladenoff, 2004; Zhu et al., 2003; Liski et al., 2005). The mixed model leads to a better representation of the mean decay of the total population, as the size of the groups of the different tree species are weighted and systematic deviations according to 43 tree species can be taken into account. The model is superior to the multiple regression model in all statistical parameters (LRT-Test, AIC, standard error of model). As the decomposition rate in this model is a quadratic function of mean annual rainfall, it also implies an optimum of decay for this variable. As can be seen in Fig. 3, these curves show an optimum within the range of 1100–1300 mm rainfall (exactly 1226 mm), which is in good accordance with the findings from Mackensen et al. (2003). They describe an optimal range between 1200 and 1300 mm. This is in line with the findings from Progar et al. (2000), which detected a declining decomposition rate with higher moisture content. On the other hand, slower decomposition can be suggested by dryness. Although the relation between log diameter and decomposition rates is discussed controversial in the literature, here a clearly negative effect on the decomposition rate is estimated. The physiological reason for this is the ratio between surface:volume, as larger ratios ease the access for decomposers and diminishing barriers for exchange processes (gas and water). Theoretically the Q10 for fungal respiration has been found to be a factor of 2–3 (Deverall, 1965). However, here we estimated this value to be only 1.6, but this value is derived from studies in situ, not in the laboratory. Furthermore the difference between mass and density losses is huge and highly significant (on a 95% confidence interval). It can be explained as losses due to fragmentation are not observed in density based studies. This is also supported from Harmon et al. (2000) and was empirically shown in Yin (1999). Including other random effects in the model (we investigated the site of the studies) did not lead to an improvement. Without the random effect on the tree species-level, the study site lead to a significant improvement compared to fixed effect model, however in this model all other climatic variables related to the study site had to be removed from the model. It was considered to be more advantageous to leave the effects that could be quantified as physiologically well-explained parameters in the model, instead of considering them as random deviations related to the study site.

6. Conclusions Overall, the development of general decay model for deadwood with a constant decay rate is advantageous against using the outcomes of case-studies, if one wants to calculate e.g. carbon flows in this pool over time. The complex process of decomposition will behave differently in ecosystems (9 studies from Europe and 30 studies from North-America), however we can suppose, that this process is determined by its physical environmental conditions and by tree-specific attributes. Using the mixed model, we are able to

911

quantify the effect on the tree species-level, for which not explicit parameters have to be estimated. With this rather small amount of parameters in a comprehensive decay model, the few observations of remaining deadwood of different studies are used effectively. Clearly, when we have a look at the residuals of our mixed model their could be other driving factors, further explaining the deviations, which are not only visible by grouping them to tree species, but also by grouping them to the site of the studies. Within the given variables and observations no more conclusions can be drawn from a statistical point of view, but the way how this analysis was set up, has the potential of further applications: to bring the complex process of decomposition to a unified mathematical expression and being able to use the information of all observations in an effective way. Acknowledgments The authors are grateful to Thomas Maschler, who helped enormously in the literature review. The study was only possible with funding from the German Federal Ministry of Food, Agriculture and Consumer Protection. References Abbott, D.T., Crossley Jr., D.A., 1982. Woddy litter decomposing following clearcutting. Ecology 63, 35–42. Agren, G.I., Bosatta, E., 1998. Theoretical Ecosystem Ecology. Understanding Element Cycles. Cambridge University Press, Cambridge, p. 234. Albrecht, L., 1991. The importance of dead wood in the forest. Forstwissenschaftl. Centralbl. 110 (2), 106–113. Barber, B.L., van Lear, D.H., 1984. Weight loss and nutrient dynamics in decomposing woody lobolly pine logging slash. Soil Sci. Soc. Am. J. 48, 906–910. Berg, B., McClaugherty, C., 2003. Plant Litter. Springer, Berlin, Heidelberg, New York, Hong Kong, London, Milan, Paris, Tokyo, p. 286. BMELV, 2007. Bundeswaldinventur - alle Ergebnisse und Berichte. http://www.bundeswaldinventur.de. Boddy, L., Swift, M.J., 1984. Wood decomposition in an abandoned beech and oak coppiced woodland in SE England. Holarctic Ecol. 7, 229–238. Böttcher, H., Freibauer, A., Obersteiner, M., Schulze, E.-D., 2008. Uncertainty analysis of climate change mitigation options in the forestry sector using a generic carbon budget model. Ecol. Model. 213, 45–62. De Vries, B.W., Kuyper, T.W., 1988. Effect of vegetation type on decomposition rates of wood in Drenthe, The Netherlands. Acta Bot. Neerl. 37 (2), 307–312. Deverall, B.J., 1965. The physical environment for fungal growth. 1. Temperature. In: Ainsworth, G.C., Sussman, A.S. (Eds.), The fungi. I. The fungal cell. Academic Press, New York, pp. 543–550. Edmonds, R.L., 1987. Decomposition rates and nutrient dynamics in small-diameter woody litter in four forest ecosystems in Washington, U.S.A. Can. J. For. Res. 17, 499–509. Edmonds, R.L., Vogt, D.J., Sandberg, D.H., Driver, C.H., 1986. Decomposition of Douglas-fir and red alder wood in clear-cuttings. Can. J. For. Res. 16, 822–831. Erickson, H.E., Edmonds, R.L., Peterson, C.E., 1985. Decomposition of logging residues in Douglas-fir, western hemlock, Pacific silver fir, and ponderosa pine ecosystems. Can. J. For. Res. 15, 914–921. Fahey, T.J., 1983. Nutrient dynamics of aboveground detritus in Lodgepole Pine (Pinus contorta ssp. latifolia) Ecosystems, Southeastern Wyoming. Ecol. Monogr. 53, 51–72. Fahey, T.J., Hughes, J.W., Pu, M., Arthur, M.A., 1988. Root decomposition and nutrient flux following whole-tree harvest of northern hardwood forest. Forest Sci. 34, 744–768. Fahey, T.J., Stevens, P.A., Hornung, M., Rowland, P., 1991. Decomposition and nutrient release from logging residue following conventional harvest of Sitka Spruce in North Wales. Forestry 64 (3), 289–301. Fogel, R., Cromack, K., 1976. Effect of habitat and substrate quality on Douglas fir litter decomposition in western Oregon. Can. J. Bot. 55, 1632–1640. Gosz, J.R., 1980. Biomass distribution and production budget for a nonaggrading forest ecosystem. Ecology 61 (3), 507–514. Gosz, J.R., Likens, G.E., Bormann, F.H., 1973. Nutrient release from decomposing leaf and branch litter in the Hubbard Brook Forest, New Hamshire. Ecol. Monogr. 43 (2), 173–191. Harmon, M.E., Franklin, J.F., Swanson, F.J., Sillins, P., Gregory, S.V., Lattin, J.D., Anderson, N.H., 1986. Ecology of coarse woody debris in temperate ecosystems. Adv. Ecol. Res. 15, 133–302. Harmon, M.E., Krankina, O.N., Sexton, J., 2000. Decomposition vectors: a new approach to estimating woody detritus decomposition dynamics. Can. J. For. Res. 30, 76–84. Herrmann, S., Prescott, C.E., 2008. Mass loss nutrient dynamics of coarse woody debris in three Rocky Mountain coniferous forests: 21 year results. Can. J. For. Res. 38, 125–132.

912

J. Zell et al. / Ecological Modelling 220 (2009) 904–912

IPCC, 2003. Good Practice Guidance for Land Use, Land-Use Change and Forestry, p. 307. Kahl, T., 2003. Abbauraten von Fichtentotholz (Picea abies (L.) Karst.) - Bohrwiderstandsmessungen als neuer Ansatz zur Bestimmung des Totholzabbaus, einer wichtigen Größe im Kohlenstoffhaushalt mitteleuropäischer Wälder. FriedrichSchiller-Universität Jena, Jena, p. 103. Kaipainen, T., Liski, J., Pussinen, A., Karjalainen, T., 2004. Managing carbon sinks by changing rotation length in European forests. Environ. Sci. Pollut. 7, 205–219. Kollmann, F., 1982. Technologie des Holzes und der Holzwerkstoffe. Bd. 1. SpringerVerlag, Berlin, Heidelberg, New York, p. 1050. Kram, P., Santore, R.C., Driscoll, C.T., Aber, J.D., Hruska, J., 1999. Application of the forest–soil–water model (PnET-BGC/CHESS) to the Lysina catchment, Czech Republic. Ecol. Model. 120, 9–30. Kruys, N., Jonsson, B.G., Stahl, G., 2002. A stage-based matrix model for decay-class dynamics of woody debris. Ecol. Appl. 12, 773–781. Laiho, R., Prescott, C.E., 1999. The contribution of coarse woody debris to carbon, nitrogen, and phosphorus cycles in three Rocky Mountain coniferous forests. Can. J. For. Res. 19, 1592–1603. Laiho, R., Prescott, C.E., 2004. Decay and nutrient dynamics of coarse woody debris in northern coniferous forests: a synthesis. Can. J. For. Res. 34, 763–777. Liski, J., Palosuo, T., Peltoniemi, M., Sievänen, R., 2005. Carbon and decomposition model Yasso for forest soils. Ecol. Model. 189 (1–2), 168–182. Mackensen, J., Bauhus, J., Webber, E., 2003. Decomposition rates of coarse woody debris—a review with particular emphasis on Australian tree species. Aust. J. Bot. 51 (1), 27–37. Mattson, K.G., Swank, W.T., Waide, J.B., 1987. Decomposition of woody debris in a regenerating, clear-cut forest in the Southern Appalachians. Can. J. For. Res. 17, 712–721. Miller, W.E., 1983. Decomposition rates of aspen bole and branch litter. Forest Sci. 29 (2), 351–356. Moore, T.R., 1984. Litter Decomposition in a Subarctic Spruce-Lichen Woodland, Eastern Canada. Ecology 65 (1), 299–308. Müller-Using, B., 2005. Totholzdynamik eines Buchenbestandes im Solling. Berichte des Forschungszentrums Waldökosysteme, Reihe A (193). Naesset, E., 1999. Decomposition rate constants of Picea abies logs in southeastern Norway. Can. J. For. Res. 29 (3), 372–381. Olson, J.S., 1963. Energy storage and the balance of producers and decomposers in ecological systems. Ecology 44 (2), 322–331. Palosuo, T., Liski, J., Trofymow, J.A., Titus, B.D., 2005. Litter decomposition affected by climate and litter quality—testing the Yasso model with litterbag data from the Canadian intersite decomposition experiment. Ecol. Model. 189 (1–2), 183–198.

Pinheiro, J.C., Bates, D.M., 2000. Mixed-Effects Models in S and S-PLUS. Springer, New York, Berlin, Heidelberg, p. 528. Progar, R.A., Schowalter, D., Freitag, C.M., Morell, J.J., 2000. Respiration from coarse woody debris as affected by moisture and saprotroph functional diversity in Western Oregon. Oecologia 124, 426–431. Ratkowsky, D.A., 1990. Handbook of Nonlinear Regression Models. Marcel Dekker Inc., New York und Basel, p. 239. Rice, M.D., Lockaby, B.G., Stanturf, J.A., Keeland, B.D., 1997. Woody debris decomposition in the Atchafalaya River Basin of Louisiana following Hurricane Disturbance. Soil Sci. Soc. Am. J. 61, 1264–1274. Scheller, R.M., Mladenoff, D.J., 2004. A forest growth and biomass module for a landscape simulation model, LANDIS: design, validation, and application. Ecol. Model. 180, 211–229. Shorohova, E.V., Shorohov, A.A., 2001. Coarse woody debris dynamics and stores in a boreal virgin spruce forest. Ecol. Bull. 49, 129–135. Swift, M.J., Healey, I.N., Hibberd, J.K., Sykes, J.M., Bampoe, V., Nesbitt, M.E., 1976. The decomposition of branch wood in the canopy and floor of a mixed deciduous woodland. Oecologia 26, 139–149. Tarasov, M.E., Birdsey, R.A., 2001. Decay rate and potential storage of coarse woody debris in the Leningrad Region. Ecol. Bull. 49, 137–147. Taylor, B.R., Prescott, C.E., Parsons, W.J.F., Parkinson, D., 1991. Substrate control of litter decomposition in four Rocky Mountain coniferous forests. Can. J. Bot. 69, 2242–2250. Thürig, E., Palosuo, T., Bucher, J.B., Kaufmann, E., 2005. The impact of windthrow on carbon sequestration in Switzerland: a model-based assessment. Forest Ecol. Manage. 210, 337–350. Tinker, D.B., Knight, D.H., 2001. Temporal and spatial dynamics of coarse woody debris in harvested and unharvested lodgepole pine forests. Ecol. Model. 141, 125–149. Verburg, P.S.J., Johnson, D.W., 2001. A spreadsheet-based biogeochemical model to simulate nutrient cycling processes in forest ecosystems. Ecol. Model. 141, 185–200. Yatskov, M., Harmon, M.E., Krankina, O.E., 2003. A chronosequence of wood decomposition in the boreal forests of Russia. Can. J. For. Res., 1211–1226. Yin, X., 1999. The decay of forest woody debris. numerical modeling and implications based on some 300 data cases from North America. Oecologia 121, 81– 98. Zhu, Z., Arp, P.A., Meng, F., Bourque, C.P.-A., Foster, N.W., 2003. A forest nutrient cycling and biomass model (ForNBM) based on year-round, monthly weather conditions, part I: assumption, structure and processing. Ecol. Model. 169, 347–360.