A new approach to simulate peat accumulation, degradation and stability in a global land surface scheme (JULES vn5.8_accumulate_soil)

stability in a global land surface scheme (JULES vn5.8_accumulate_soil) Sarah E. Chadburn1, Eleanor J. Burke2, Angela V. Gallego-Sala3, Noah D. Smith1, M. Syndonia Bret-Harte4, Dan J. Charman3, Julia Drewer5, Colin W. Edgar4, Eugenie S. Euskirchen4, Krzysztof Fortuniak6, Yao Gao7, Mahdi Nakhavali3, Włodzimierz Pawlak6, Edward A.G. Schuur8, and Sebastian Westermann9 1Department of Mathematics, University of Exeter, Exeter, UK 2Met Office Hadley Centre, Exeter, UK 3Geography Department, University of Exeter, Exeter, UK 4Institute of Arctic Biology, University of Alaska Fairbanks, USA 5UK Centre for Ecology & Hydrology, Bush Estate, Penicuik, Scotland, UK 6Department of Meteorology and Climatology, University of Lodz, Poland 7Finnish Meteorological Institute, Helsinki, Finland 8Center for Ecosystem Science and Society, Northern Arizona University, Flagstaff, USA 9Department of Geosciences, University of Oslo, Oslo, Norway Correspondence: Sarah Chadburn (s.e.chadburn@exeter.ac.uk)


Testing the interpolation method
Here we show an experiment to test the method of interpolating soil carbon back onto the original model layers after the soil column height is changed to account for changes in soil carbon. We use a second order Taylor approximation: where C is the soil carbon density (kgm −3 ) and δz = z −z ef f , where z is the position of the original soil layer that we want to interpolate back on to, and z ef f is the effective layer thickness after incrementing the soil carbon in a given model timestep. This setup is illustrated in Figure S1 (black and red annotations). Figure S1 illustrates an example soil layer, which we call layer d.
The first order derivative (second term on RHS of Equation 1) is approximated as the gradient of the line between either the points d and d + 1, or the points d − 1 and d, depending on whether δz is positive or negative (respectively). i.e.
dC(z ef f,d ) dz or dC(z ef f,d ) dz These are indicated in green on Figure S1 The second order derivative (third term on RHS of Equation 1) at z ef f,d is the gradient of the gradient at z ef f,d , and is approximated as the difference in gradients (grad plus and grad minus ) above and below z ef f,d , divided by an interval δz , i.e.
It is the choice of δz that is interesting. Naively, we would take δz as 0.5(z ef f,d+1 − z ef f,d−1 ) (see delta1 on Figure S1). However, by testing different sizes of δz we found that taking δz = 2|δz| (see delta2 on Figure S1) results in the closest approximation to the true solution for the soil carbon profile. We tested this by running a toy model in R. This was a simplified 'toy' version of the RothC carbon model used in JULES, with just two carbon pools with different turnover times, one of which receives external input of carbon and has a shorter turnover time (analogous to the litter pools) and the other receives carbon from pool 1 and has a longer turnover time (analogous to the biomass and humous pools). We input the equivalent of 200 gCyr −1 , a 'typical' NPP value globally, and turn over pool 1 at a rate of 1 −4 day −1 and pool 2 turns over at a rate 5× smaller, which is equivalent to around 27 years and 137 years as turnover times, which are reasonable values. Litter input is distributed exponentially through the profile with an efolding depth (equivalent to τ lit ) of 5m −1 , as in JULES. Figure S1: Schematic of interpolation method. Figure S2: Testing interpolation method for low resolution soil carbon profile. Note that 'hi res 2nd order' (green) uses delta2 and is therefore the high resolution equivalent of the magenta line. The blue and purple lines are coincident. Figure S2 shows the results. First of all, we tried running this model with a high spatial resolution of 1cm layers. With this high resolution, the first order derivative approximates the true function well, as adding the second order term has little impact on the simulated carbon profile. In the lower resolution simulation (layer thicknesses as in JULES), using only the first order term the sharp peak in the profile is notably 'smeared' out, and using delta1 (see Fig. S1) in the second order term does not improve this, i.e. the blue and purple lines in Figure S2 are on top of each other. Using delta2 (see Fig. S1) on the other hand captures the peak much better (magenta line on Figure S2). There is still some smearing apparent in the deeper layers, but using a smaller value of δz (delta2 ×0.3) begins to lead to numerical instabilities near the surface as shown by the 'zig zag' appearance of the red line towards the top of the profile. Thus we chose delta2 as the pragmatic 'best' solution, although we may consider improving this in future work.

Testing individual processes
As well as the simulations shown in the main manuscript, initially the scheme to account for the space taken up by soil carbon -which allows a peat layer to build up in the profile (Section 2.3 in main manuscript) -was switched on with fixed soil properties (i.e. without the interaction between carbon and soil properties). We called this version 'JULES-acc', with acc standing for accumulatesee Supplementary Table S1. Figure S3 shows these simulations in which the soil column is allowed to grow with the addition of organic material (see Section 2.3 in the main manuscript and supplementary Table S1 for details), along with the original version of JULES. Note that in these simulations, the soil parameters are not dynamic and we use the fixed soil thermal/hydraulic parameters that were originally derived for the sites. Thus we demonstrate the impact of accounting for the volume of the organic material ( Figure  S3). We can see that the soil carbon profiles are much more realistic in JULES-acc than in JULES, and almost all of the JULES-acc simulations shown fall within the interquartile range of the observational dataset.
In order to quantify the improvement we take the root mean squared error (RMSE) between the median soil carbon profiles for each climate zone, shown in Table S2. The version with our chosen decomposition function and parameters, JULES-accC, has a RMSE that is reduced by 40% for temperate peatland sites and by almost 80% for boreal peatland sites (RMSE reducing from 37.3 kgm −3 to 8.2 kgm −3 , Table S2).
We see that the age at the soil surface was typically too high in the original version of JULES ( Figure  S4). In the new version with soil accumulation enabled (JULES-accum), the age at the soil surface is better simulated, and age throughout the profile ranges from a good match with the observations (falling within the interquartile range) to being somewhat too low, particularly for temperate sites in JULES-accB. This motivated the change of parameter τ resp from JULES-accB to JULES-accC which leads to greater carbon ages at depth, and generally the accC(10) simulations fall more within the observed quartiles in terms of age than do JULES-accB. We chose τ resp = 2 since this had been shown in [KRS + 13] to give the best correspondence with soil carbon age-depth profiles in their simulation with the Community Land Model.
The difference between the two JULES-accC simulations is that in JULES-accC10 the distribution of litter inputs into the soil is more weighted towards the surface (τ lit = 10 as opposed to 5, Table  S1). This may make the model more likely to accumulate an organic layer at the surface. While the JULES-accC10 simulations don't compare as well to the large dataset of peat cores ( Figure S3, Table  S2), when compared against individual site observations (Supplementary Figure S6, Organic Sites in Table S2), this simulation is actually much better, so we tested both versions for the JULES-Peat simulations.  Table S1: JULES simulations conducted. Note that T means 'True' (or the process is switched on) and F means 'False' (process switched off). F→T mean that the process was switched off during spinup and on during the main run. The 'Decomp. function' refers to changing from the original to the new decomposition function show in Figure 1F in the main manuscript, and we also changed a switch that was using the total soil moisture instead of the unfrozen soil moisture, so that the new decomposition function is a function of the unfrozen soil moisture (more realistic, since frozen water is not available for microbes to use). Where the Initial C is given as 'Peat' we initialise all spinups with the spun-up profile for Auchencorth from JULES-Peat-W, with ∼1.5m of peat, and otherwise initalise the model with zero soil carbon. Simulation name Accumulate Decomp. fn τ resp τ lit Soil dyn. Lat. flow Init.   Figure S5: Evaluation of simulated soil carbon profile at mineral soil sites.