Dr Alasdair Sykes 2021-04-20
This repository contains inputs, code and outputs for soil carbon
modelling work performed to accompany the Legumes Translated project.
The model input data compares a set of control rotations to a series of
alternative legume-based rotations; model outputs detail expected soil
carbon stock change as a result of these different management
strategies. All code in the root directory is reproducible with the
contents of this repository; external data read-in is handled in the
raw-data
directory.
Soil carbon modelling is performed using an R package-based implementation of the IPCC (2019) Tier 2 steady-state soil carbon model for national level greenhouse gas reporting. Full details of the model methodology can be found here; documentation of the corresponding R package implementation can be found here.
Spatial data used in the modelling process is (1) climate data from the CRU project, available open-source here here, and (2) soil data from the SoilGrids project, available open-source here. Spatial definition of the study sites is based on NUTS 2017 Level 2 data, with corresponding shapefiles available here.
-
01-spatial-data-processing.R
Processing of the study spatial data, resulting in study site-specific outputs of climate and soil data. -
02-model-data-assembly.R
Reading of the raw supplied data files and construction of the raw model data. -
03-model-input-processing.R
Conversion of the assembled raw data input model inputs and creation of the soil carbon model runs using thesoilc.ipcc
package. -
04-model-output-processing.R
Summarisation of the model outputs into summary datasets and plots. -
05-model-output-analysis.R
Exploratory analysis primarily aimed at drawing out the driving factors behind C stock differences between treatment and control rotations.
-
model-data Data directory containing model input and output datasets.
-
model-output-summaries Data directory containing summarised model outputs.
-
raw-data Data directory containing raw spatial data, raw study data as supplied, and data translation scripts.
-
spatial-data Data directory containing processed spatial data.
-
intermediate-data Output data directory containing human-readable combined raw and model input data files in .csv format. Not used in analysis.
For each study site, the spatial data was cropped and masked against the relevant NUTS level 2 shapefile, and arithmetic mean was calculated for mean monthly temperature, precipitation, potential evapotranspiration, and soil sand fraction. The spatial climate data represented 20-year monthly mean conditions. These conditions were repeated for the duration of the model runs.
Data referencing crop management and manure application practices was
taken as supplied and matched up to analogous crop and manure categories
in the relevant IPCC (2019) methodology. The soilc.ipcc
package was
used to build model inputs accordingly, and to compute the model runs.
For each study site and control-treatment pair (control = no legumes in the rotation, treatment = legumes in the rotation) the model was run three times; once for equilibrium control conditions, once for equilibrium treatment conditions, and once for a transition from control to treatment conditions. For each model run, the model was run in for a period of 50 years with averaged conditions (climate, soil and organic inputs) followed by a period of 500 years in the relevant rotation. For the transition from control to treatment, the model was run in in mean control conditions, calculated for year 1–50 in the control rotation, and transitioned to treatment conditions for years 51–500.
Model results were summarised to calculate a mean soil carbon stock for the final 50 years of equilibrium conditions in control and treatment rotations, and to find a 20- and 50-year mean carbon stock change post transition in the transitional run.
Results from the model transitional runs are presented in Fig 1.
The summary statistics for the model runs are as follows:
Control | Treatment | Control stocks (tonnes C ha-1) | Treatment stocks (tonnes C ha-1) | Annual stock change, 20-year (tonnes C ha-1) | Annual stock change, 50-year (tonnes C ha-1) |
---|---|---|---|---|---|
01_-leg | 01_+leg | 55.62 | 45.28 | -0.08 | -0.05 |
02_-leg | 02_+leg | 54.58 | 49.38 | -0.05 | -0.02 |
03_-leg | 03_+leg1 | 52.43 | 48.13 | -0.01 | -0.04 |
03_-leg | 03_+leg2 | 52.43 | 45.59 | -0.07 | -0.05 |
04_-leg | 04_+leg | 64.65 | 52.71 | -0.10 | -0.06 |
05_-leg | 05_+leg | 65.14 | 73.72 | 0.09 | 0.05 |
06_-leg | 06_+leg | 48.14 | 59.21 | 0.08 | 0.06 |
07_-leg | 07_+leg | 72.42 | 64.37 | -0.04 | -0.02 |
08_-leg | 08_+leg | 79.87 | 63.15 | -0.16 | -0.10 |
09_-leg | 09_+leg | 53.26 | 53.56 | 0.00 | 0.00 |
10_-leg | 10_+leg | 70.21 | 61.81 | -0.10 | -0.06 |
11_-leg | 11_+leg | 84.19 | 73.98 | -0.12 | -0.07 |
12_-leg | 12_+leg | 93.17 | 74.65 | -0.10 | -0.10 |
13_-leg | 13_+leg | 92.42 | 71.11 | -0.13 | -0.14 |
14_-leg | 14_+leg | 80.54 | 70.32 | -0.06 | -0.05 |
15a_-leg | 15a_+leg1 | 45.20 | 42.79 | -0.02 | -0.01 |
15a_-leg | 15a_+leg2 | 45.20 | 45.70 | -0.01 | 0.00 |
15b_-leg | 15b_+leg1 | 38.46 | 32.81 | -0.05 | -0.03 |
15b_-leg | 15b_+leg2 | 38.46 | 32.45 | -0.05 | -0.03 |
16_-leg | 16_+leg | 77.83 | 41.29 | -0.21 | -0.16 |
17_-leg | 17_+leg1 | 38.13 | 37.19 | -0.01 | -0.01 |
17_-leg | 17_+leg2 | 38.13 | 39.64 | 0.01 | 0.01 |
18_-leg | 18_+leg1 | 100.25 | 113.68 | 0.10 | 0.04 |
18_-leg | 18_+leg2 | 100.25 | 93.96 | -0.15 | -0.07 |
18_-leg | 18_+leg3 | 100.25 | 93.96 | -0.15 | -0.07 |
18_-leg | 18_+leg4 | 100.25 | 96.09 | -0.14 | -0.06 |
18_-leg | 18_+leg5 | 100.25 | 45.49 | -0.47 | -0.29 |
18_-leg | 18_+leg6 | 100.25 | 104.99 | -0.09 | -0.02 |
19a_-leg | 19a_+leg | 55.94 | 60.34 | 0.07 | 0.04 |
19b_-leg | 19b_+leg | 34.84 | 36.91 | 0.01 | 0.01 |
20_-leg | 20_+leg | 77.66 | 56.30 | -0.21 | -0.14 |
Some analysis has been performed in order to determine the driving factors behind C stock response to the change in rotations. Initially, correlation between the key model input and output variables for treatment (legume) rotations was computed (Fig. 2).
Fig. 2. Correlation plot for treatment rotation variables. Variabile
total_y
indicates the size of the total soil C pool in tonnes per
hectare.
While this plot shows the driving factors behind equilbrium soil C for
the treatment variables, these treatments are also following some
potentially quite different control rotations. To further interrogate
the data, a similar correlation-based approach was taken, with the focus
on the net difference (treatment - control
) between the model input
and output variables. Note that some environmental variables are
identical for control and treatment, and that therefore differential
correlation was not possible to compute.
Interpretation of Fig. 3 suggests that differences in equilibrium carbon
stocks are largely driven by differences in total carbon inputs
c_input
, differences in nitrogen and lignin fractions n_frac
and
lig_frac
, and to a lesser extent differences in tillage practices
till_fac
. Resulting calculated model variables alpha
, k_a
and
k_s
carry this through to the final total_y
(tonnes C
ha-1) estimate.
Fig. 4 shows the major differences in model variables where results following transition to treatment rotations were positive vs. negative for resulting soil carbon stock changes.
Fig. 4. Boxplot showing differences between model variables for positive and negative soil carbon responses.
Interpretation of Fig. 4 confirms the indication of the correlation plots, but suggests the loss of soil C, where present, is primarily due to differences in organic matter inputs to the soil from the different rotations.
To aid further interpretation/analysis of the model results, an
intermediate data directory has been added (intermediate-data/
) which
contains human readable combined raw and model input data files.
All queries relating to this work should be addressed to the author at [email protected].