|
Thomas E. Hilinski
Department of Soil and Crop Sciences
Colorado State University
January 2001
An important component in calculating budgets of terrestrial soil carbon is an estimate of total soil organic carbon. The CENTURY model simulates organic C in the top 20 cm of soil only, so a method for projecting the depth distribution of organic C is necessary for estimates of the profile total. This document describes the implementation of the exponential depth distribution of organic C in CENTURY soil organic matter model. This model replaces the original CENTURY version 5 algorithm for distributing C below the simulation layer. Implementation of this algorithm also required changes to the algorithms for the layered soil, erosion, and deposition.
Replacement of the original C distribution algorithm was necessary in order to overcome problems whose source was in the assumptions underlying the design of the algorithm. These problems affected simulation results in the following ways:
Advantages of the new algorithm include its ability to approximate observed depth distributions of soil C, a direct link between the C distribution and the size of the soil C pools used in CENTURY, and easy calculation of the storage of C over any depth interval. Disadvantages are increased computation time, and error sensitivity to C density profiles that deviate substantially from an exponentially decreasing form.
Description of the vertical distribution of soil organic C as a function which decays exponentially with depth has been described by several authors desiring to model their profile data. For example, Bernoux et al. (1998) compared power and exponential decay functions to Amazon forest soils using regression, and found the exponential model fit best, with R2 values greater than 0.82. Their model is based upon the depth distribution equation of Arrouays and Pelissier (1994) for temperate forest soils. Elzein and Balesdent (1995) used a numerical simulation of organic matter decay incorporating an exponential depth distribution model to successfully simulate organic matter decay and transport in a temperate forest soil. Rosenbloom (1997) and Rosenbloom et al. (2001) used a similar form to extrapolate C at depth based upon the potential maximum C at the surface, in a model of soil erosion. Harden et al. (2001) used the same form as Rosenbloom (1997) to distribute input of plant and decomposition C to the soil over the profile depth.
CENTURY dynamically simulaties soil organic C (SOC) in the top 20 cm of the soil. SOC is divided into several pools depending upon turnover rate of SOC and the type of material the C is in, i.e., plant lignin, microbial material, or soil particulate C. Below the simulation layer is a set of identical pools whose size is dependent upon (1) their initial values specified in the site parameters for the simulation, and (2) in the course of a simulation, the depth distribution of C, determined by the algorithms specified here.
Soil organic C is distributed in the physical soil according to an exponential density distribution, expressed in a equation modified from Bernoux et al. (1998) and Rosenbloom (1997):
C(z) = Cb + (C0 - Cb) * exp(-K * z)
|
Eq. 1 |
where
C(z)
|
C density (g cm-3) at depth z (cm, z ³ 0), the maximum C density. | ||
Cb
|
C density at bottom of profile (g cm-3), the minimum C density. Profile depth = zmax. | ||
C0
|
C density at the surface (g cm-3) | ||
K
|
scale constant (cm-1) |
The behaviour of the scaling constant K has the primary effect on the distribution of C and the total profile C. The scaling constant K is the primary control of the total profile C. K decreases with increasing profile depth and decreasing C density at the profile bottom (figure 1). The total profile C content, determined by the area under the curve, is inversely proportional to K.
|
|
|
|
| Figure 1. K vs. profile depth. | Figure 2. C density depth distribution for different values of K. |
In order to implement this equation, the algorithm needs to find the values of K, C0, and Cb repeatedly, using the information provided by CENTURY: total organic C from the simulation layer lower layer pools. The simulation layer is at 0-20 cm and the lower layer is below this depth. With deposition the simulation layer can extend to 30 cm depth. Erosion can subsequently reduce this to no less than 20 cm. The profile depth must be as least 20 cm, and may be as deep as is needed for estimation of profile C storage.
The density distribution equation has three unknowns in only one equation, so a few assumptions are needed to provide further information and generate two additional equations. First, assume the density at ½ of the profile depth is 10 percent of the surface density. Second, assume the density at the bottom is 10 percent of the average lower layer density. Two additional relationships are needed. Layer C density multiplied by layer thickness gives the total C for the layer. The total C for the profile is provided by CENTURY as the sum of simulation C and lower layer C.
Now we have enough information to build three equations for three unknowns:
C(zmax) - (0.1 * CENTURY lower layer C density) = 0
|
Eq. 2a | |
C(z = 5 cm) - 1.5 * CENTURY simulation layer density = 0
|
Eq. 2b | |
å(C(z) * z) - total profile C = 0
|
Eq. 2c |
Equation 2a is the difference in estimated and known values of Cb, and is solved for Cb. Equation 2b is the difference in estimated and known values of C0, and is solved between for C0. Equation 2c is the difference in estimated and known values of the total profile C (g cm-2), and is solved for parameter K. Equation 2b uses a depth of 5 cm in order to avoid abnormally high estimates of the surface C density.
Integrating the density equation (equation 1) over two depths (z1 < z2) gives functions for the amount of C in a layer, and for the total profile C. These functions are shown in equations 3 and 4.
Layer C = [ exp(-K*z1) - exp(-K*z2) ] * (C0 - Cb) / K + (z2 - z1) * Cb
|
Eq. 3 | |
Profile C = [ 1 - exp(-K * zmax) ] * (C0 - Cb) / K + zmax * Cb
|
Eq. 4 |
Iterative numerical methods require initial estimates of the variables for which a solution is being sought. Here, C0 is set to 50 percent higher than the simulation layer C density. Cb is set to 10 percent of the lower layer density. To obtain an initial estimate of scale constant K, we assume the density at half the profile depth is 10% of surface density. Plugging this assumption into equation 1 and solving for the K provides the following equation:
K = -2.0 / zmax * ln (Ratio)
|
Eq. 5 |
where
Ratio = (0.1 * C0 - Cb) / (C0 - Cb)
|
Multidimensional root finding normally requires a numerical technique such as the Newton-Raphson method, however these methods work best when the initial guesses are close to the root, otherwise the methods can diverge. I tested the Newton-Raphson method with these equations and found it to be unstable. Instead, I developed an algorithm which assumes the estimates of C0 and Cb provide reasonable values (equations 2a and 2b) which remain fixed, and searches iteratively for a value of K, which used to calculate the total profile C (equation 4), matches the known profile total (the sum of CENTURY's simulation layer and lower layer C).
Potential sources of error using this algorithm include the following:
A test of this algorithm is available using field data from a recent study of erosion on agricultural lands. Harden et al. (1999) examined eroded and non-eroded sites in Mississippi, within the Mississippi River basin, with the intent to quantify the net loss or storage of organic C loss from erosional-depositional environments. Using data from their non-eroded Baptist Cemetery soil profile, the solution of these equations (calculated independently of CENTURY) provides the parameter values: K = 0.0625 cm-1, C0 = 0.0587 g cm-3, Cb = 2.20e-4 g cm-3. Graphs of the data and the modified Rosenbloom equation using these parameters is shown in figures 3 and 4. The data and the model compare favorably. Total profile C is identical for data and model, as it should be since the algorithm fits for total profile C.
|
|
| Figure 3. C density by soil layer for the Baptist Cemetery profile, comparing data with model. | Figure 4. The depth distribution of C for the Baptist Cemetery profile, comparing data with model. |
The object used by CENTURY to manage its soil, TCenturySoil, also manages the depth distribution of organic C. An initial distribution is calculated prior to the start of the simulation, based upon the sum of the initial soil organic C pools, the starting simulation layer thickness, and the sizes of the lower layer C pools, initialized from site parameters (see the function TCentury::PreliminaryCalcs). At the end of each simulation year, the distribution is recalculated, based upon the mean monthly sum of the soil C pools, and the current values for the simulation layer thickness and the lower layer C pools (see the function TCentury::EndOfYearTasks).
There are new output variables relating to the soil C distribution, shown in the following table:
Variable Name Description Units SOILDEPTH Depth of soil profile. cm LHSOMTC Lower horizon total soil organic C. g m-2 SOMTPC Total soil organic C in the soil profile = SOMSC + LHSOMTC. g m-2
Arrouays, D., Pelissier, P., 1994. Modeling carbon storage profiles in temperate forest humic loamy soils of France. Soil Sci. 157, 185-192.
Bernoux, M., Arrouays, D., Cern, C.C., Bourennane, H., 1998. Modeling vertical distribution of carbon in oxisols of the western Brazilian Amazon (Rondonia). Soil Science 163, 941-951.
Elzein, A., Balesdent, J., 1995. Mechanistic simulation of vertical distribution of carbon concentrations and residence times in soils. Soil Sci. Soc. Am. J. 59, 1328-1335.
Harden, J.W., Sharpe, J.M., Parton, W.J., Ojima, D.S., Fries, T.L., Huntington, T.G., Dabney, S.M., 1999. Dynamic replacement and loss of soil carbon on eroding cropland. Gbl. Biogeochem. Cycles, 13, 885-901.
Harden, J.W., Fries, T.L., Pavich, M., 2001. Evolution of hillslopes in Iowa and implications for carbon cycling. Gbl. Biogeochem. Cycles, (in press).
Rosenbloom, N., 1997. A hillslope evolution model for the coupled prediction of soil texture and ecosystem dynamics. Ph. D. Dissertation, Univ. of Colorado, Boulder. 177 pp.
Rosenbloom, N., Doney S.C., Schimel D.S, 2001. Geomorphic evolution of soil texture and organic matter in eroding landscapes. Biogeochem. (in press).
Copyright 2001 Thomas E. Hilinski. All rights reserved.