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 Laikipia County, Kenya, 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 Laikipia Conservancies Association (LCA) to support investment decision in a landscape restoration project that aims to sequester carbon in rangelnads.

Materials and methods

Study background

The LCA project 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 Laikipia. The goal of LCA 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 and security.

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 in the Borana example.
#>                      variable distribution upper median lower
#> 1                     n_years        const  40.0     NA  40.0
#> 2                  general_CV    tnorm_0_1   0.4     NA   0.2
#> 3        mean_annual_rainfall      posnorm 750.0     NA 400.0
#> 4    grazing_intensity_borana   tnorm_0_1    0.8     NA   0.6
#> 5       fire_frequency_borana        const   0.0     NA   0.0
#> 6    percentage_lignin_borana      posnorm  20.0     NA   5.0
#> 7 percentage_cellulose_borana      posnorm  40.0     NA  25.0
#> 8      sand_percentage_borana      posnorm  35.0     NA  10.0
#> 9         baseline_soc_borana      posnorm  28.9     NA  13.2

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){
  
#Borana Conservancy
maximum_aboveground_production<-vv((0.84*mean_annual_rainfall-27.5)*
                                   (1.33-0.0075*sand_percentage_borana),
                                   general_CV,n=n_years)
proportion_leaf_biomass<-vv(0.6+0.24*grazing_intensity_borana,
                            general_CV,n=n_years)
leaf_area_index<-vv(proportion_leaf_biomass/0.6-0.015*
                 exp(4.6*grazing_intensity_borana),general_CV,n=n_years)
grazer_modified_production<-leaf_area_index*
                            maximum_aboveground_production
estimate_aboveground_productivity<-maximum_aboveground_production*
                                   proportion_leaf_biomass/0.6-0.015*
                                   exp(4.6*grazing_intensity_borana)
belowground_production<-vv(917.4-0.763*mean_annual_rainfall,
                           general_CV,n=n_years)
lignocellulosic_content_borana<-vv((percentage_lignin_borana/
                                percentage_cellulose_borana),general_CV,n=n_years)
plant_derived_carbon<-0.45*(lignocellulosic_content_borana*
                      grazer_modified_production*
                      (1-grazing_intensity_borana)*(1-fire_frequency_borana)+
                      (lignocellulosic_content_borana+0.05)*
                      belowground_production)
dung_derived_carbon<-lignocellulosic_content_borana*0.45*
                     grazing_intensity_borana*grazer_modified_production
wetdays<-vv((0.0004*mean_annual_rainfall-0.025)*240,
            general_CV,n=n_years)

#Interpreting SOC at equilibrium to be equal to baseline SOC
maximum_microb_resp_rate<-wetdays*(0.7+0.3*sand_percentage_borana/100)*
                          (0.00044*baseline_soc_borana-0.579)
change_soc_borana<-(plant_derived_carbon+dung_derived_carbon-
                    maximum_microb_resp_rate)/100

#Applying a sigmoid function to adjust for soc sequestration potential overtime   
change_soc_borana_adjusted<-gompertz_yield(max_harvest=change_soc_borana,
                            time_to_first_yield_estimate= 1,
                            time_to_second_yield_estimate= 5,
                            first_yield_estimate_percent=0.20*100,
                            second_yield_estimate_percent=0.50*100,
                            n_years,
                            var_CV = 0,
                            no_yield_before_first_estimate = TRUE)

return(list(change_soc_borana_adjusted=change_soc_borana_adjusted))
}

Simulation of carbon stock changes

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) over a 40 year period, is as shown below

Monte Carlo outputs
variable X0. lower_bound X10. X25. X50. X75. X90. upper_bound X100. mean
y.change_soc_borana_adjusted1 0.10 0.17 0.22 0.30 0.39 0.49 0.56 0.66 0.77 0.39
y.change_soc_borana_adjusted2 0.13 0.23 0.29 0.41 0.53 0.66 0.76 0.91 1.04 0.53
y.change_soc_borana_adjusted3 0.17 0.29 0.37 0.52 0.68 0.85 0.97 1.15 1.34 0.68
y.change_soc_borana_adjusted4 0.21 0.36 0.46 0.64 0.84 1.04 1.20 1.42 1.63 0.83
y.change_soc_borana_adjusted5 0.24 0.42 0.54 0.76 0.98 1.22 1.40 1.66 1.94 0.98
y.change_soc_borana_adjusted6 0.28 0.48 0.61 0.86 1.12 1.39 1.60 1.89 2.20 1.12
y.change_soc_borana_adjusted7 0.31 0.53 0.68 0.96 1.25 1.55 1.78 2.12 2.46 1.24
y.change_soc_borana_adjusted8 0.34 0.58 0.74 1.04 1.35 1.68 1.94 2.31 2.67 1.35
y.change_soc_borana_adjusted9 0.36 0.63 0.80 1.12 1.46 1.80 2.09 2.49 2.87 1.45
y.change_soc_borana_adjusted10 0.38 0.67 0.85 1.18 1.54 1.92 2.20 2.61 3.05 1.54
y.change_soc_borana_adjusted11 0.40 0.69 0.88 1.24 1.61 2.00 2.31 2.74 3.17 1.61
y.change_soc_borana_adjusted12 0.42 0.72 0.91 1.29 1.68 2.07 2.40 2.87 3.29 1.67
y.change_soc_borana_adjusted13 0.43 0.74 0.94 1.32 1.73 2.13 2.49 2.92 3.40 1.72
y.change_soc_borana_adjusted14 0.44 0.77 0.97 1.36 1.77 2.20 2.53 3.00 3.48 1.77
y.change_soc_borana_adjusted15 0.45 0.78 0.99 1.38 1.81 2.23 2.57 3.05 3.56 1.80
y.change_soc_borana_adjusted16 0.46 0.79 1.00 1.41 1.84 2.27 2.63 3.08 3.60 1.83
y.change_soc_borana_adjusted17 0.46 0.80 1.02 1.42 1.87 2.31 2.64 3.13 3.63 1.85
y.change_soc_borana_adjusted18 0.47 0.81 1.03 1.44 1.89 2.33 2.68 3.20 3.67 1.87
y.change_soc_borana_adjusted19 0.47 0.81 1.03 1.45 1.89 2.34 2.71 3.20 3.73 1.89
y.change_soc_borana_adjusted20 0.47 0.82 1.04 1.47 1.91 2.36 2.73 3.23 3.74 1.90
y.change_soc_borana_adjusted21 0.48 0.82 1.04 1.47 1.93 2.37 2.73 3.25 3.75 1.91
y.change_soc_borana_adjusted22 0.48 0.83 1.06 1.48 1.92 2.39 2.75 3.25 3.80 1.92
y.change_soc_borana_adjusted23 0.48 0.84 1.06 1.48 1.92 2.39 2.76 3.27 3.81 1.93
y.change_soc_borana_adjusted24 0.48 0.84 1.06 1.48 1.94 2.40 2.77 3.31 3.83 1.93
y.change_soc_borana_adjusted25 0.48 0.84 1.06 1.49 1.94 2.42 2.77 3.30 3.83 1.94
y.change_soc_borana_adjusted26 0.49 0.84 1.07 1.50 1.94 2.41 2.78 3.32 3.86 1.94
y.change_soc_borana_adjusted27 0.49 0.84 1.07 1.50 1.94 2.41 2.78 3.32 3.84 1.95
y.change_soc_borana_adjusted28 0.49 0.84 1.07 1.50 1.96 2.43 2.79 3.33 3.88 1.95
y.change_soc_borana_adjusted29 0.49 0.84 1.07 1.50 1.95 2.42 2.79 3.33 3.81 1.95
y.change_soc_borana_adjusted30 0.49 0.84 1.07 1.50 1.96 2.43 2.79 3.31 3.86 1.95
y.change_soc_borana_adjusted31 0.49 0.85 1.07 1.51 1.96 2.43 2.79 3.34 3.82 1.95
y.change_soc_borana_adjusted32 0.49 0.84 1.07 1.50 1.96 2.43 2.80 3.31 3.83 1.95
y.change_soc_borana_adjusted33 0.49 0.84 1.07 1.50 1.96 2.43 2.79 3.33 3.84 1.95
y.change_soc_borana_adjusted34 0.49 0.85 1.07 1.50 1.96 2.43 2.81 3.31 3.87 1.96
y.change_soc_borana_adjusted35 0.49 0.85 1.07 1.51 1.96 2.43 2.81 3.31 3.86 1.96
y.change_soc_borana_adjusted36 0.49 0.84 1.08 1.50 1.96 2.43 2.80 3.32 3.84 1.96
y.change_soc_borana_adjusted37 0.49 0.84 1.07 1.50 1.96 2.43 2.80 3.33 3.84 1.96
y.change_soc_borana_adjusted38 0.49 0.85 1.07 1.51 1.96 2.43 2.81 3.31 3.88 1.96
y.change_soc_borana_adjusted39 0.49 0.84 1.08 1.51 1.95 2.43 2.79 3.35 3.91 1.96
y.change_soc_borana_adjusted40 0.49 0.84 1.07 1.51 1.97 2.42 2.80 3.34 3.83 1.96

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
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.
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

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/