Introduction

Soil organic carbon (SOC) in grasslands and savannas represents one of the largest reservoirs of carbon on earth (Lal et al., 2007). However, grassland management practices, such as use of fire, fertilization, re-vegetation and restoration, etc., that can lead to net sequestration of carbon are still poorly understood (Ritchie, 2014). In Maasai Mara, herbivores can have dramatically different effects on soil organic carbon (SOC), both positive and negative, depending on soil type, precipitation, plant species composition and grazing intensity (McSherry & Ritchie, 2013). To model soil carbon dynamics in grasslands we use the SNAP model (Ritchie, 2014) developed within a decisionSupport package (Luedeling & Gohring, 2017; Luedeling & Whitney, 2018). The model can make predictions of soil carbon stocks based on mean rainfall, plant lignin and cellulose, soil texture, fire history, and grazing intensity conditions over an extended period of time. We use Maasai Mara ecosystem as a case study to support decision making for One Mara Carbon Project (OMCP).

Materials and methods

Study background

The One Mara Carbon Project (OMCP) is an ecosystem-level conservation initiative that aims to protect the ecological resources and provide economic and social co-benefits for local communities in the Maasai Mara. The goal of OMCP is to reduce emissions, restore the degraded landscape, inform grassland management options, provide alternative livelihoods for communities, build resilient conservancies, increase biodiversity conservation efforts, and enhance tourism.

Decision support package

The decisionSupport () function in the R package of the same name requires two inputs:

  1. An input_table (in .csv format) specifying the names and probability distributions for all variables used in the analysis.To generate the input_table, we first describe the input variables for the SNAP model. The SNAP model (initially developed using data from Serengeti National Park) requires five input variables: mean annual rainfall (mm/yr); grazing intensity (1-(grazed biomass/ungrazed biomass)) (McNaughton, 1985); fire frequency (number of fires recorded/number of years over which fires were monitored); lignin and cellulose content (proportion) in aboveground herbaceous vegetation, and sand percentage. These variables constitute the input_table as shown below.
Monte Carlo inputs
Description lower upper distribution variable
Number of years to run the simulation (n) 10.00 10.00 const n_years
Coefficients for introducing annual variation 0.20 0.40 tnorm_0_1 general_CV
Carbon markets discount rate (%) 8.00 13.00 posnorm carbon_discount_rate
Price per ton of carbon ($) 8.00 20.00 posnorm price_ton_carbon
Mean annual rainfall (mm/yr) 450.00 900.00 posnorm mean_annual_rainfall
Grazing intensity (1-(grazed biomass/ungrazed biomass)) 0.25 0.90 tnorm_0_1 grazing_intensity
Fire frequency (n/yrs) 0.25 0.75 tnorm_0_1 fire_frequency
Lignin and cellulose content (proportion) 0.15 0.35 tnorm_0_1 lignin_and_cellulose_content
Sand percentage (%) 25.00 70.00 posnorm sand_percentage
Sampling depth (cm) 30.00 30.00 const sampling_depth
Soil organic carbon concentration (g C kg-1) 0.45 0.54 posnorm soil_organic_cab_conc
Bulk density (kg m-3) 1.48 1.54 posnorm bulk_density
Conservancy interventions cost 1000.00 2500.00 posnorm interventions_cost

To generate the input data for the variables (including those with missing information), the variable distributions are first described by a 90% confidence interval, which are specified by lower (5% quantile) and upper (95% quantile) bounds. To achieve this, project proponents undergo a calibration training (Hubbard, 2014) on how to generate variable estimates, for which they are 90% confident that the actual values lie within the provided ranges. This helps to include variables with missing information in the analysis. The process is usually backstopped by subject matter experts, literature review and actual field measurements. For the shapes of variable distribution, const describes variables that are constant throughout, norm describes variables with normal distribution, tnorm_0_1 describes variables with a truncated normal distribution that can only have values between 0 and 1 (useful for probabilities) and posnorm describes variables with normal distribution truncated at 0 (only positive values allowed).

  1. An R_function () that predicts decision outcomes based on the variables named in input_table. This function is customized by the user to address a particular decision problem. For this analysis, the function grazing_soil_carbon_dynamics () describing the causal relationships between the SNAP model variables (Ritchie, 2014) was used.

Loading data in R

We could simply start translating SNAP model causal relationships into a mathematical decision model now. However, the model function will be designed to make use of variables provided to it externally (random numbers drawn according to the information in the input_table). To do this effectively, we need to load the data in R and define sample values for all variables and test pieces of the function code during the model development process. This is accomplished with the following helper function make_variables ().

make_variables<-function(est,n=1)
{ x<-random(rho=est, n=n)
    for(i in colnames(x)) assign(i,
     as.numeric(x[1,i]),envir=.GlobalEnv)
make_variables(estimate_read_csv("Soil_carbon_dynamics.csv"))
}

This function isn’t included in the decisionSupport package, because it places the desired variables in the global environment. This isn’t allowed for functions included in packages on R’s download servers. Applying make_variables to the input_table (with default setting n=1) generates one random number for each variable, which then allows us to easily test the code we’re developing.

Mathematical model development

The following is our draft version of the SNAP model, with explanations of all causal relationships (we will format this better once we agree on the final model)

grazing_soil_carbon_dynamics <- function(x, varnames){
  
#Eq.(1)
maximum_aboveground_production<-vv((0.84*mean_annual_rainfall-27.5)*
                                   (1.33-0.0075*sand_percentage),
                                   general_CV,n=n_years)
#Eq. (4) 
proportion_leaf_biomass<-vv(0.6+0.24*grazing_intensity,
                            general_CV,n=n_years)
#Eq.(3)
leaf_area_index<-vv(proportion_leaf_biomass/0.6-0.015*
                    exp(4.6*grazing_intensity),general_CV,n=n_years)
#Eq.(2)
grazer_modified_production<-leaf_area_index*
                            maximum_aboveground_production
#Eq.(5)
estimate_aboveground_productivity<-maximum_aboveground_production*
                                   proportion_leaf_biomass/0.6-0.015*
                                   exp(4.6*grazing_intensity)
#Eq. (6)
belowground_production<-vv(917.4-0.763*mean_annual_rainfall,
                           general_CV,n=n_years)
#Eq.(7)
plant_derived_carbon<-0.45*(lignin_and_cellulose_content*
                            grazer_modified_production*
                           (1-grazing_intensity)*(1-fire_frequency)+
                           (lignin_and_cellulose_content+0.05)*
                            belowground_production)
#Eq.(8)
dung_derived_carbon<-lignin_and_cellulose_content*0.45*
                     grazing_intensity*grazer_modified_production
#Eq.(9)
wetdays<-vv((0.0004*mean_annual_rainfall-0.025)*240,
             general_CV,n=n_years)

# Baseline SOC calculations
baseline_carbon_stocks<-vv(sampling_depth*soil_organic_cab_conc*
                           bulk_density*(1-sand_percentage/100),
                           general_CV,n=n_years)
#Eq.10
maximum_microb_resp_rate<-vv(wetdays*(0.7+0.3*sand_percentage/100)*
                            (0.00044*baseline_carbon_stocks-0.579),
                             general_CV,n=n_years)
#Eq.11
change_carbon_stocks<-plant_derived_carbon+dung_derived_carbon-
                      maximum_microb_resp_rate

carbon_revenues<-change_carbon_stocks*price_ton_carbon-interventions_cost

carbon_stocks_NPV<-discount(carbon_revenues,carbon_discount_rate,calculate_NPV=TRUE)

return(list(change_carbon_stocks=change_carbon_stocks
            cashflow_carbon_stocks_NPV=carbon_revenues,
            carbon_project_NPV=carbon_stocks_NPV))
}

Simulation of carbon stock changes and revenues

To run the model, the grazing_soil_carbon_dynamics () function, along with the data from the input_table, are fed into the Monte Carlo simulation (MC) function (Luedeling and Whitney, 2018; Luedeling et al. 2021) to conduct the full analysis. See the vignette guides on applying the mcSimulation () function (Fernandez et al. 2021). Below is the code we used to perform the Monte Carlo simulation with 10,000 model runs.

SOC_simulation<- mcSimulation(
                 estimate = estimate_read_csv("Soil_carbon_dynamics.csv"),
                 model_function = grazing_soil_carbon_dynamics,
                 numberOfModelRuns = 1e4, # 10,000 runs
                 functionSyntax = "plainNames"
)

Results

A tabulated summary of the results describing the 90% confidence interval of changes in SOC stocks (Mg C ha-1) and revenues (USD), over a 10 year period, is as shown below

Monte Carlo outputs
variable X0. lower_bound X10. X25. X50. X75. X90. upper_bound X100. mean
y.change_carbon_stocks1 4.6e-01 8.5e-01 9.4e-01 1.1 1.2 1.4 1.5 1.6 2.2 1.2
y.change_carbon_stocks2 4.6e-01 8.5e-01 9.4e-01 1.1 1.2 1.4 1.5 1.6 2.2 1.2
y.change_carbon_stocks3 4.6e-01 8.5e-01 9.4e-01 1.1 1.2 1.4 1.5 1.6 2.2 1.2
y.change_carbon_stocks4 4.6e-01 8.5e-01 9.4e-01 1.1 1.2 1.4 1.5 1.6 2.2 1.2
y.change_carbon_stocks5 4.6e-01 8.5e-01 9.4e-01 1.1 1.2 1.4 1.5 1.6 2.2 1.2
y.change_carbon_stocks6 4.6e-01 8.5e-01 9.4e-01 1.1 1.2 1.4 1.5 1.6 2.2 1.2
y.change_carbon_stocks7 4.6e-01 8.5e-01 9.4e-01 1.1 1.2 1.4 1.5 1.6 2.2 1.2
y.change_carbon_stocks8 4.6e-01 8.5e-01 9.4e-01 1.1 1.2 1.4 1.5 1.6 2.2 1.2
y.change_carbon_stocks9 4.6e-01 8.5e-01 9.4e-01 1.1 1.2 1.4 1.5 1.6 2.2 1.2
y.change_carbon_stocks10 4.6e-01 8.5e-01 9.4e-01 1.1 1.2 1.4 1.5 1.6 2.2 1.2
y.cashflow_carbon_stocks_NPV1 -3.7e+03 -2.5e+03 -2.3e+03 -2039.1 -1737.3 -1424.4 -1138.6 -965.1 -68.0 -1734.6
y.cashflow_carbon_stocks_NPV2 -3.7e+03 -2.5e+03 -2.3e+03 -2039.1 -1737.4 -1424.4 -1138.6 -965.0 -68.0 -1734.6
y.cashflow_carbon_stocks_NPV3 -3.7e+03 -2.5e+03 -2.3e+03 -2039.1 -1737.4 -1424.4 -1138.7 -964.9 -67.9 -1734.6
y.cashflow_carbon_stocks_NPV4 -3.7e+03 -2.5e+03 -2.3e+03 -2039.1 -1737.3 -1424.4 -1138.6 -965.0 -68.0 -1734.6
y.cashflow_carbon_stocks_NPV5 -3.7e+03 -2.5e+03 -2.3e+03 -2039.1 -1737.4 -1424.4 -1138.7 -965.0 -68.0 -1734.6
y.cashflow_carbon_stocks_NPV6 -3.7e+03 -2.5e+03 -2.3e+03 -2039.1 -1737.3 -1424.4 -1138.6 -965.1 -68.0 -1734.6
y.cashflow_carbon_stocks_NPV7 -3.7e+03 -2.5e+03 -2.3e+03 -2039.1 -1737.3 -1424.4 -1138.7 -965.1 -67.9 -1734.6
y.cashflow_carbon_stocks_NPV8 -3.7e+03 -2.5e+03 -2.3e+03 -2039.1 -1737.3 -1424.4 -1138.7 -965.1 -67.9 -1734.6
y.cashflow_carbon_stocks_NPV9 -3.7e+03 -2.5e+03 -2.3e+03 -2039.1 -1737.4 -1424.4 -1138.6 -965.0 -67.9 -1734.6
y.cashflow_carbon_stocks_NPV10 -3.7e+03 -2.5e+03 -2.3e+03 -2039.1 -1737.4 -1424.4 -1138.7 -965.0 -67.9 -1734.6
y.carbon_project_NPV -2.7e+04 -1.7e+04 -1.6e+04 -13615.7 -11498.5 -9470.0 -7525.2 -6408.0 -458.9 -11561.7

Also shown, is a graphical representation of the distributions of carbon stock changes in the first year of implementing the project.

Figure 1: Projected changes in soil organic carbon stocks in the first year of project implementation

Model sensitivity analysis

To analyse model sensitivity we perform a Partial Least Squares Regression (PLSR) analysis of Monte Carlo simulation results. We do this by applying a post-hoc analysis to the mcSimulation() outputs with the plsr.mcSimulation() function in the decisionSupport package. We use the function to determine the Variable Importance in the Projection (VIP) score and coefficients of a Projection to Latent Structures (PLS) regression model. This function uses the outputs of the mcSimulation() selecting all the input variables from the grazing_soil_carbon_dynamics () decision analysis function in the parameter object and then runs a PLS regression with an outcome variable defined in the parameter resultName. This resulted in the graph shown below.

Figure 2: Important variables as determined by VIP analysis of PLS regression models. The green bar indicate a positive correlation of the uncertain variable with the outcome variable. Gray bars indicate variables that did not meet the threshold of the model sensitivity analysis.

Discussion and conclusion

The model results projected here still need refining using data collected in each conservancy. These efforts will help inform project proponents on the most viable intervention options. We also need to include exponential equations where baseline carbon stocks are very low, to avoid getting negative values for microbial respiration in eq. 10. The model doesn’t include risk variables to account for losses due to drought and heat extremes. However, these factors often lead to reduced SOC inputs from vegetation and increased exposure of soils to wind and water erosion (Okin et al., 2018). Increasing frequency of other extremes, such as intense rainfall events, can also lead to surface sealing preventing moisture infiltration, reduced plant growth and ultimately a reduction in soil carbon stocks. These factors need to be sufficiently captured to mimic real world scenarios, but the model has been validated with recent studies in the Northern Rangeland Trust and previous studies in Serengeti National Park (Ritchie, 2014). These studies resulted in the development of the methodology for measuring SOC stock changes in grasslands (VM0032 registered under the VERRA standards).

References

Lal R, Follett F, Stewart BA, Kimble JM. 2007. Soil carbon sequestration to mitigate climate change and advance food security. Soil Science 172:943956 https://doi.org/10.1097/ss.0b013e31815cc498

Ritchie, M. E. (2014). Plant compensation to grazing and soil carbon dynamics in a tropical grassland. PeerJ, 2, e233. https://doi.org/10.7717/peerj.233

McSherry M, Ritchie ME. 2013. Effects of grazing on grassland soil carbon density: a global review. Global Change Biology 19:13471357 https://doi.org/10.1111/gcb.12144.

Luedeling, E., Gohring, L. (2017). R package decisionSupport: Quantitative support of decision making under uncertainty. https://CRAN.R-project.org/package=decisionSupport

Luedeling, E., Whitney, C. W.(2018). R package decisionSupport: Controlled burns in conifer forests. https://cran.r-project.org/web/packages/decisionSupport/vignettes/wildfire_example.html

McNaughton SJ. 1985. Ecology of a grazing ecosystem: The Serengeti. Ecological Monographs 55:259294 https://doi.org/10.2307/1942578.

Luedeling, Eike, Lutz Goehring, Katja Schiffers, Cory Whitney, and Eduardo Fernandez. 2021. decisionSupport: Quantitative Support of Decision Making Under Uncertainty. http://www.worldagroforestry.org/.

Fernandez, Eduardo, Cory Whitney, and Eike Luedeling. 2021. “Applying the mcSimulation Function in decisionSupport.” https://cran.r-project.org/web/packages/decisionSupport/vignettes/example_decision_function.html.

Okin, G. S., Sala, O. E., Vivoni, E. R., Zhang, J., and Bhattachan, A. (2018). The Interactive Role of Wind and Water in Functioning of Drylands: What Does the Future Hold. BioScience 68 (9), 670–677. https://doi.org/10.1093/biosci/biy067

Ritchie, M. E. (2015). VM0032 Methodology for the Adoption of Sustainable Grasslands through Adjustment of Fire and Grazing, v1.0. Verified Carbon Standard. https://verra.org/methodology/vm0032-methodology-for-the-adoption-of-sustainable-grasslands-through-adjustment-of-fire-and-grazing-v1-0/