1. Introduction

Human activities have caused an increased atmospheric concentrations of carbon dioxide (CO₂) and other heat-trapping gases. The average annual increase rose from 1.7 ppm between 1993 and 2003 to 2.1 ppm over the 2003-2012 decade. The concentration of atmospheric C rose to 393.31 ppm by March 2013 [1].

Land use change has had a significant impact on global C stocks with cultivation reported to cause significant depletion of organic matter and releasing carbon dioxide (CO₂) into the atmosphere [2]. Some of the major causes of CO₂ release from the earth to the atmosphere are deforestation and degradation which are driven mostly by agricultural expansion and shifting cultivation [3], production of charcoal and fuel wood [4], legal and illegal timber logging [5] construction and wild fires. Agricultural activities that release C from the soil into the atmosphere include tillage and other forms of soil disturbances facilitating gaseous exchange between the soil and the atmosphere. Soil disturbances also enable the incorporation of plant materials into the soil [6].

Although agricultural activities have been identified as major sources of CO₂ emissions, it is also possible to have agricultural activities that are adapted to reverse these negative effects and promote soil carbon sequestration in addition to other benefits of food security and sustainability of systems [7]. Soil C sequestration can offer a valuable offset for greenhouse gas emissions in agriculture, forestry and other land uses (AFLOU) [8]. Lal [9] stated that developing countries have the greatest potential for soil C sequestration since most of the soils are highly degraded and therefore below C saturation. Agricultural practices that reduce emissions and promote C storage can be designed to achieve sustainable land use but the need for simple, rapid monitoring methods cannot be over emphasised. A number of scientific evidence gaps are linked to the accuracy of C accounting, ascribed to a lack of data and uncertainties related to C storage and C flux models [10].

Soil carbon models have been developed for prediction and provision of information on the rate of soil C sequestration or loss. Several models have been designed and reviewed by Smith et al. [11]. Among these are RothC, DNDC, CENTURY and DAISY all based on conceptual carbon pools with different turnover rates. The RothC model is among the models which have been identified by FAO [12] as a widely applicable easy to use model. The model has been applied in arable and non- arable soils [13-16] in many parts of the world although with limited applicability in Southern Africa.

Changes in the rates of soil C sequestration due to changing environmental or management factors, can take several years to become apparent [17]. Therefore, future impact of agriculture activities and associated land use change can only be predicted by the use of observations in combination with models which can provide a means of evaluating the changing practices in the future. Currently, there is limited information on the potential for C sequestration in sandy and clayey soils of Zimbabwe especially on smallholder farmers’ fields despite the importance of SOC in soil biological activities which in turn affect soil physical and chemical properties.

In this study, we test the applicability of the RothC model for predicting SOC in agricultural and natural forest soils of Zimbabwe. The model was run on seven selected common land management practices encompassing tillage intensities and fertility regimes. The objectives of the study were to: (1) compare estimates of soil C changes from measured and simulated C in tillage and fertility treatments, (2) assess the potential of different land management practices on future soil C storage (3) assess the sensitivity of Roth C model to temperature rise and (4) compare equilibrium levels from RothC with equilibrium levels estimated using the Langmuir equation.

2. Materials and Methods

2.1. Study Sites

Modelling of tillage impacts was done on soils representing clayey and sandy soils from Hereford in Bindura district (17° 42΄ S; 31° 44΄ E), Nyarukunda in Shamva district (17° 00΄S; 31° 43΄E) and Murewa district latitude (17° 39' 13" S and longitude 31° 48' 30" E). A description of the sites is given by Mujuru et al. [18] and Mujuru et al. [19]. Briefly, Hereford soils are red clays varying from silty clay loam to clay, with characteristics corresponding to Rhodic Ferralsols [20, 21] and falling into the category of low activity clays [22]. Sandy soils, are derived from coarse granite covering almost 70% of Zimbabwe [23] and are classified as being in the Kaolinitic order, Fersiallitic group under the Zimbabwean soil classification, which corresponds to Ferric Luvisols [23-25] but using IPCC default classes derived from the harmonised world soils database [22] they can be classified them as Arenosols, (>70% sand and <8% clay) and are broadly referred to as sandy soils. Murewa soils are granitic sands (Haplic Arenosols) [21] which are strongly weathered having low levels of available nutrients and low nutrient reserves. These are interspaced with pockets of dolerite intrusions that give rise to small patches of relatively fertile clays (Chromic Luvisols) [20, 25]. Bindura represented clayey soils whereas Shamva represented sandy soils for tillage assessments whereas Murewa represented both clayey and sandy soils under fertility treatments.

2.2. Modelling of Carbon Sequestration Potential

Data for SOC for whole soil and density fractions under tillage and fertility treatments were obtained from Mujuru et al. [18] and Mujuru et al. [19] respectively (Table 1). The RothC model was applied for the prediction of SOC in each land management practice and evaluated with measured values obtained in 2010. The RothC model used in this study is based on monthly time step calculations that simulate SOC turnover over periods ranging from a few years to a few centuries [26, 27]. The model inputs were (a) climatic data (monthly rainfall (mm), monthly evapotranspiration (mm), average monthly mean air temperature (°C), Climate data were obtained from world climate database collection of meteorological data [28]. The ETO calculator [29] was used to estimate potential evapotranspiration based on temperature and rainfall of each specific location (b) Soil data (clay%, initial soil organic carbon (SOC) stock (Mg C ha^-1), depth of the soil layer (30 cm), inert organic matter (IOM) (Table 1) approximated using equation (1) proposed by Falloon et al, [30] because the radiocarbon content was not known.

(1)

and (c) land use and land management data (soil cover, monthly input of plant residues (Mg ha^-1), monthly input of farmyard manure (FYM) (Mg C ha^-1), residue quality factor (decomposable plant material (DPM)/resistant plant material (RPM) ratio) [31]. Soil cover was based on whether the soil is bare or vegetated in a particular month and is indicated as either covered or fallow [31].

Table 1. Inert organic matter (IOM), clay% and SOC stocks in whole soil and density fractions in agricultural lands and natural forest systems (Mean (SD)) to a depth of 30 cm

The plant materials in RothC are subdivided into DPM and RPM, whereas plant debris in soils is found in fLF and oLF and forms of dissolved organic carbon (DOC). SOC is split into four active pools DPM, RPM, microbial biomass (BIO) and humified organic matter (HUM) which decompose by a first-order process, each with its own characteristic rate, and an amount of inert organic matter (IOM) resistant to decomposition (Figure 1). Plant C inputs are assumed to exclusively enter the DPM and RPM in proportions, which are dependent on the source of the plant materials (Figure 1). The process algorithms in RothC are affected by the three climatic factors with each pool decomposing by first order kinetics with characteristic decomposition rates modified by soil moisture, temperature and plant cover.

Figure 1. Structure of the Rothamsted Carbon Model [31]

The model apportions plant litter input between DPM and RPM depending upon the vegetation type. For most agricultural crops and improved grassland, a DPM/RPM ratio of 1.44 was used, i.e. 59% DPM and 41% RPM whereas DPM/RPM ratio of 0.25 was used for deciduous or tropical woodland i.e. 20% as DPM and 80% RPM. In the model, the proportion that goes to CO₂, BIO and HUM is determined by the amount of clay in the soil [31]. Although RPM estimated using the RothC model correlated well with measured LFC in sandy soils and natural forests, LFC in the cropland clayey soils neither correlated with RPM nor DPM in croplands. Because of this discrepancy the model default decomposition rate constants for active compartments were RPM (0.3), DPM (10.0) BIO (0.66) and HUM (0.02) [31].

The model was run to equilibrium (10 000 years) after estimating the required plant carbon input for each treatment distributing the C inputs into compartments (DPM, RPM, BIO and HUM) with different decomposition rates. The assumption implied in this inverse simulation procedure, was that RothC could simulate the dynamic changes in SOC under the specified conditions. The data of C and radiocarbon age in the different compartments obtained in equilibrium mode were used to run the model in short term mode from 1850 to 2050. The three measured fractions are particulate organic matter found outside of aggregates (fLF), POM found within aggregates (oLF), and a mineral-associated fraction (MaHF). Model estimated pools of RPM + DPM + BIO were compared with LF C (fLF + oLF) whilst HUM +IOM were compared with MaHF for each land management system. The RPM, DPM, BIO, LF represent amounts of C input.

The scenarios (A-G) used for modelling were selected with consideration to the current farming systems and included:

[A] - conventional tillage (CF) with maize legume rotation – consists of an ox drawn plough to a depth of 15-20 cm once before planting, residues were removed and the remaining biomass incorporated into the soil during ploughing in the next season,

[B] - Ripping (RP) - minimum tillage with an ox drawn ripper to a depth of 15-20 cm with maize-legume rotation, crop residues were to be retained in the field after harvesting, ground cover of 2.5-3.0 Mg ha^-1 was required in RP,

[C] - Direct seeding (DS) - no tillage using an ox drawn direct seeder with synchronised seeding and fertiliser application with maize-legume rotation, residues were also retained or supplemented to achieve the 2.5-3.0 Mg ha^-1 ground cover and cropping in each tillage system was maize (Zea maize. L) /cowpeas (Vigna unguiculata L. Walp) rotation at Nyarukunda or soy bean (Glycine max L. Merr) at Hereford, each treatment received annual basal fertiliser of 165 kg ha^-1 compound D (i.e. 11kg ha^-1 N, 10 kg ha^-1 P, 10 kg ha^-1 K), which was followed by 69 kg ha^-1 N applied as ammonium nitrate in splits at 4 and 7 weeks after germination [32, 33],

[D] - Natural forest (NF),

[E] - Conventional tillage with continuous maize cropping and no fertility amendments (control),

[F] - Conventional tillage under continuous maize cropping with annual addition of nitrogen fertiliser (N Fert),

[G] - Conventional tillage with continuous maize cropping with a combination of nitrogen fertiliser and cattle manure (N Fert + manure) where ammonium nitrate supplies 100 kg N ha^-1 and cattle manure applied at 5 Mg ha^-1 supplied an equivalent 10 kg P ha^-1 and 0.9 % N each cropping season.

Each scenario was run for sandy and for clayey soil under current and changing temperature conditions. The SOC stock predicted by RothC model were compared with the SOC data for 2010 [18, 19].

2.3. Comparison of Equilibrium SOC Levels from RothC and Langmuir Equation

SOC equilibrium levels were estimated using the Langmuir equation and compared with equilibrium levels estimated by RothC model. The Langmuir equation can be used to evaluate the adsorption of light fraction C (LF C) onto mineral surfaces, and becoming mineral associated heavy fraction C (MaHF C) and is regarded as sequestered C [34]. We assumed that over time, the LF C decomposes and in part becomes adsorbed onto mineral soil particles as the MaHF C. In addition, soil minerals can randomly adsorb LF C until the MaHF has reached C saturation. Therefore, interaction between LF C and soil minerals follows adsorption and desorption processes that can be described using the Langmuir equation. Equation (2) shows the linearization used to fit the data following Yin and Cai [35] and Bolter and Hornberger [36].

(2)

Where MaHFC_max is the maximum adsorption capacity for organic C (equilibrium value for soil organic C in the MaHF) and k is the equilibrium constant. LFC/MaHFC versus LFC yields a linear relationship with slope 1/(MaHFC_max) and intercept 1 / (k MaHFC). The equilibrium level estimated using Langmuir equation was compared with the equilibrium level obtained from RothC model.

2.4. Data Analysis

Relationships between mean values (± SD) for modelled and measured SOC and equilibrium C levels in each land management system were analysed using linear regression. The goodness of fit was tested using correlation coefficient (r) and root mean square error (RMSE) in addition to standard deviation. RMSE shows the percentage term for the total difference between predicted and observed values. The bias value was calculated as Yi – Xi where Yi = measured SOC or equilibrium level and Xi = modelled SOC or equilibrium level. A t-test was conducted to compare measured with modelled C. Significant differences were tested at p≤ 0.05.

3. Results and Discussion

3.1. Relationship between Measured and Modelled SOC Stocks

The RothC model is designed to simulate soil C turnover in the soil using user estimates of the C inputs making the evaluation of the soil C turnover components easier. Inverse simulation techniques allow the determination of input needed to match the observed SOC. Running the RothC in reverse mode estimated the inputs required to attain the SOC stocks in 2010 for the seven scenarios in clayey and sandy soils. Carbon inputs on clayey soils with continuous conventional tillage, under control, N Fert and N Fert + manure in clayey homefields were higher (2.4, 4.5, 5.0 Mg C ha^-1 respectively) than outfields (2.8, 3.8, 4.6 Mg C ha^-1 respectively). Sandy soils followed similar trends of being higher in homefields (0.9, 1.5, 1.6 Mg C ha^-1 respectively) than outfields (0.4, 1.0, 0.9 Mg C ha^-1 respectively). Other models such as CENTURY do not use user defined C inputs as in the RothC model, because plant growth is simulated. The C inputs in RothC model are therefore not affected by simulation errors which can be carried over in further simulations. In sandy soils, the R² values were lower than clayey soils probably because of the short term measurements and associated small SOC storage and capacity to change.

Simulation modelling using the RothC model, showed good agreement between the simulated and observed SOC in treatments on sandy soils while on clayey soils, control, N fertiliser and a combination of N fertiliser and manure showed good agreement with small deviations (Figure 2). Highest measured SOC was found in natural forests on whilst lowest was in control treatments of both tillage and fertility treatments. Measured SOC on clayey soil was NF > RP > CF > N Fert + manure > DS > N Fert > Control whereas the modelled SOC was: NF > N Fert + manure > N Fert > RP > DS > CF > Control. The higher SOC in NF reflects the changes associated with conversion of forests to croplands. It was also clear that the cropping systems are characterised by lower SOC when compared with natural forests on similar soils [37]. Such C losses support the findings of Scholes and Hall [38] who showed that approximately 50% of SOC is lost in the first 20 years following conversion of tropical woodland, grassland or savannah. Such losses can be evident even after 5 years [9].

Figure 2. The relationship between observed C stocks and modelled C stocks from RothC carbon model. * ------ = Sandy soils, ▲▬ = Clayey soils

Although tillage is considered a major cause of SOC loss in agricultural systems [6], the model did not predict a large difference between tilled and untilled sites. Instead, fertility soil amendments, mainly a combination of N fertiliser and manure, can be more beneficial for SOC sequestration. Similar trends in modelled and measured value were found in sandy soils (Table 2).

Table 2. Soil organic carbon (SOC) measured and predicted using site specific soil input values for the 0-30 cm depth in seven land management practices on clayey and sandy soils

Under the current climate scenario, there is potential for additional carbon storage in all land management practices with greater magnitudes on clayey soils. Compared with measured SOC RothC underestimated all SOC values except CF and DS on sandy soils which were overestimated (Table 2). Clayey soils especially Bindura sites had higher deviations (2−34%) than sandy soils which had 2−13% in agricultural systems and 45% in NF. The t-test however, showed no significant difference between modelled and measured values on sandy soils whilst in clayey soils there was significant difference (p=0.01, R² = 0.84).

In cases where the RothC estimates were higher than measured (CF and DS on sandy soils), there is possibility of misrepresentation of clay content coupled with a general lower trend as a function of soil depth which was higher than the recommended 23 cm. RothC estimates lower than measured values could be a reflection of a potential insensitivity of the model to tillage. Although the model was developed in Rothamsted, UK, based on the prevailing farming systems there, in which soil was ploughed and stubble incorporated soon after crop harvest, the simulated SOC matched the observed SOC, suggesting that large fraction of the organic matter inputs was lost probably through crop uptake before entering soil to be sequestrated as SOC.

The regression model combining both clayey and sandy soils performed better than separated analysis y = 0.7153x + 2.3894; R² = 0.854 whereas individual comparisons of observed and modelled C in clayey and sandy soils had R² = 0.709; p<0.02; SE = 4.71 and R² = 0.593; p = 0.05; SE = 5.28 respectively (Figure 2).

3.2. Relationship between Measured SOC Pools and Conceptual RothC SOC Pools

Density fractionations were compared with model outputs of RPM+DPM and HUM+ IOM. Estimated resistant plant material carbon pool (RPM) in the RothC model and measured light fraction carbon (LFC) content were compared. The conceptual pools estimated using the RothC model for each site are shown in Table 3.

Table 3. Values of DPM, RPM, BIO and HUM compartments obtained from RothC model for the year 2010

Although there was a good relationship between light fraction C and RPM (R² = 0.674, p < 0.001, SD= 0.185, SS = 1.25) (Figure 3a) and between MaHF C and HUM (R² = 0.797, p<0.01, SD=3.74, SS = 824.71), they did not match perfectly for substitution in the model but MaHF C and HUM + IOM had a perfect match (R² = 0.98, p< 0.01) (Figure 3b).

Figure 3. Relationship between (a) measured LF C and model RPM and (b) MaHF C vs modelled HUM + IOM in sandy and clayey soils

Annual rates of carbon sequestration in the range of 0.001 – 0.006 Mg ha^-1 yr-¹ and 0.001 - 0.02 Mg ha^-1 yr^-1 in sandy (Figure 4) and clayey soils (Figure 5) respectively could be obtained over the next 40 years (2010-2050). The highest SOC sequestration rate was on clayey soil under combined N Fert + manure and RP (0.020 Mg ha^-1 yr^-1) whereas on sandy soils DS had highest accumulation rate (0.007 Mg ha^-1 yr^-1). In clayey soils mean rate of C accumulation from 2010 to 2050 was 0.011 Mg ha^-1 yr^-1 whereas on sandy soils mean C accumulation rate was 0.003 Mg ha^-1 yr^-1). The rate of C sequestration in these systems is lower than the general figure of 0.1-0.2 Mg ha^-1 yr^-1 proposed for conservation tillage in drylands by Lal [9].

Figure 4. Predicted SOC in seven land management practices under current climate on sandy soils. CF = conventional farming with maize legume rotation, RP = Minimum tillage with maize legume rotation, DS = no tillage with maize legume rotation, NF = Natural forest, Control = Conventional tillage under continuous maize cropping (no fertility amendments), N Fert = conventional tillage under continuous maize cropping (nitrogen fertiliser 100 kg ha-1) and N Fert + manure = Conventional tillage under continuous maize cropping (nitrogen fertiliser + 5t cattle manure)

Figure 5. Predicted SOC in seven land management practices under current climate on clayey soils. CF = conventional farming with maize legume rotation, RP = Minimum tillage with maize legume rotation, DS = no tillage with maize legume rotation, NF = Natural forest, Control = Conventional tillage under continuous maize cropping (no fertility amendments), N Fert = conventional tillage under continuous maize cropping (nitrogen fertiliser 100 kg ha-1) and N Fert + manure = Conventional tillage under continuous maize cropping (nitrogen fertiliser + 5t cattle manure)

When the model was run under changing climate with temperature rise of 1.5°C rise, clayey soils showed similar trends for NF but greater potential for DS than the two fertility treatments (manure and N fertiliser) (Figure 6 and 7) whereas on sandy soils, fertility treatments had higher accumulation than DS and RP but lower than NF.

Figure 6. Fitted and projected SOC stocks on clayey soil up to year 2050 under temperature rise of 1.5°C. CF = conventional farming with maize legume rotation, RP = Minimum tillage with maize legume rotation, DS = no tillage with maize legume rotation, NF = Natural forest, Control = Conventional tillage under continuous maize cropping (no fertility amendments), N Fert = conventional tillage under continuous maize cropping (nitrogen fertiliser 100 kg ha-1) and N Fert + manure = Conventional tillage under continuous maize cropping (nitrogen fertiliser + 5 Mg ha -1 yr-1 cattle manure)

Figure 7. Fitted and projected SOC stocks on sandy soil up to year 2050 under temperature rise of 1.5°C. CF = conventional farming with maize legume rotation, RP = Minimum tillage with maize legume rotation, DS = no tillage with maize legume rotation, NF = Natural forest, Control = Conventional tillage under continuous maize cropping (no fertility amendments), N Fert = conventional tillage under continuous maize cropping (nitrogen fertiliser 100 kg ha-1) and N Fert + manure = Conventional tillage under continuous maize cropping (nitrogen fertiliser + 5 Mg ha-1 yr-1cattle manure)

Although the N fertiliser sites accumulated more C than no till site (DS) other studies found inorganic fertilisers detrimental to C storage e.g. Farage et al. (2007) apart from the costs associated with its acquisition. Increased temperatures result in increased plant production resulting in higher biomass and C. However, soil respiration also increase with increasing temperature but the balance between these two processes will make the difference. In this assessment, modelling with RothC suggests that increased biomass production may be larger than increased decomposition resulting in increased C stocks.

In all cases, modelling of SOC in clayey and sandy soils of Zimbabwe showed that a new steady state will be reached if the current practices are maintained, and so subsequent declines in SOC become relatively small with time. A 1.5°C rise in temperature shows different responses in clayey and sandy soils with the clayey soils having greater C accumulation whereas in sandy soils, modelled C stocks are below current stocks. Clayey soils benefit from increased temperatures whereas sandy soils tend to decline. Under the changing climate scenario the potential for additional carbon storage is limited in all land management practices on sandy soils. The model has shown that under current climatic conditions all systems except the natural forest on clayey soils have reached steady state whereas a 1.5°C rise in temperature causes some of the systems on clayey soils to sequester more C. The results also show that when holding all the other factors constant, the model is sufficiently sensitive to a rise in global temperatures with sandy soils reaching an equilibrium much earlier than clayey soils.

3.3. Comparison of Equilibrium Level Estimation by Langmuir Equation vs RothC Model

The Langmuir equation could not estimate equilibrium levels for the three fertility systems (control, N fertiliser and N fertiliser + manure) due to poor model fit (Table 4) shown by R² values below 0.50. A comparison of equilibrium levels estimated using the RothC model and Langmuir equation shows positive correlation between the two methods on clay soils (R² = 0,87, P<0.01, SE = 4.86) (Figure 8). In clayey soils, measured SOC in CF, DS, N Fert and N Fert + manure were below equilibrium levels whereas RP, control and natural forest reached equilibrium. In sandy soils all systems were below equilibrium levels estimated by RothC as the Langmuir could not be applied to sandy soils due to poor model fit.

Figure 8. Relationship between equilibrium levels estimated by RothC model and the Langmuir equation

Table 4. Equilibrium C levels estimated by Langmuir equation for land management practices on clayey and sandy soils

It is important to note that uncertainties of up to 20% can be found due to insufficient information about soil and climate (Poussart et al., 2004). Therefore, since the data used in this study was obtained from short term experiments (4-9 years), results must be treated with some caution as they may require further ground validation from long-term experimental trials in addition to quantification of model output uncertainty. The only problem is that it may be difficult to find long term experiments for this purpose in the southern African region.

4. Conclusions

RothC model is one of the most widely used models for the estimation and prediction of SOC stock on agricultural and forest land due to its successful evaluations in the past and the generally good availability of the input data required. The model was successfully used for estimating SOC stock under selected land management practices on clayey and sandy soils of Zimbabwe. On the basis of our results, it can be concluded that RothC 26.3 model pool of HUM + IOM is related to the MaHF measured from density fractionation and the model can be used to estimate SOC stock changes on Zimbabwean agricultural and forest soils. There was a good relationship between equilibrium levels estimated by RothC model and those estimated using the Langmuir equation. The model also showed that under current climatic conditions all systems except the natural forest on clayey soils have approached steady state. Holding other factors constant, a 1.5°C rise in temperature causes some of the systems on clayey soils to sequester more C than the current. The model is therefore sufficiently sensitive to a rise in global temperatures with sandy soils reaching an equilibrium much earlier than clayey soils. The modelling approach represents one of the most promising methods for the estimation of SOC stock changes and allowed us to evaluate the changes in SOC in the past period on the basis of measured data. However, since the data used was for short term experiments (4-9 years), results must be treated with some caution as they may require further ground validation from long-term experimental trials in addition to quantification of model output uncertainty.

ACKNOWLEDGEMENTS

We appreciate assistance from Mr. Admore Mureva.