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

Peatlands have often been neglected in Earth system models (ESMs). Where they are included, they are usually represented via a separate, prescribed grid cell fraction that is given the physical characteristics of a peat (highly organic) soil. However, in reality soils vary on a spectrum between purely mineral soil (no organic material) and purely organic soil, typically with an organic layer of variable thickness overlying mineral soil below. They are also dynamic, with organic layer thickness and its properties changing over time. Neither the spectrum of soil types nor their dynamic nature can be captured by current ESMs.

Here we present a new version of an ESM land surface scheme (Joint UK Land Environment Simulator, JULES) where soil organic matter accumulation – and thus peatland formation, degradation and stability – is integrated in the vertically resolved soil carbon scheme. We also introduce the capacity to track soil carbon age as a function of depth in JULES and compare this to measured peat age–depth profiles. The new scheme is tested and evaluated at northern and temperate sites.

This scheme simulates dynamic feedbacks between the soil organic material and its thermal and hydraulic characteristics. We show that draining the peatlands can lead to significant carbon loss, soil compaction and changes in peat properties. However, negative feedbacks can lead to the potential for peatlands to rewet themselves following drainage. These ecohydrological feedbacks can also lead to peatlands maintaining themselves in climates where peat formation would not otherwise initiate in the model, i.e. displaying some degree of resilience.

The new model produces similar results to the original model for mineral soils and realistic profiles of soil organic carbon for peatlands. We evaluate the model against typical peat profiles based on 216 northern and temperate sites from a global dataset of peat cores. The root-mean-squared error (RMSE) in the soil carbon profile is reduced by 35 %–80 % in the best-performing JULES-Peat simulations compared with the standard JULES configuration. The RMSE in these JULES-Peat simulations is 7.7–16.7 kg C m⁻³ depending on climate zone, which is considerably smaller than the soil carbon itself (around 30–60 kg C m⁻³). The RMSE at mineral soil sites is also reduced in JULES-Peat compared with the original JULES configuration (reduced by ∼ 30 %–50 %). Thus, JULES-Peat can be used as a complete scheme that simulates both organic and mineral soils. It does not require any additional input data and introduces minimal additional variables to the model. This provides a new approach for improving the simulation of organic and peatland soils and associated carbon-cycle feedbacks in ESMs.

1 Introduction

Peatlands are extremely carbon-dense ecosystems, occupying only around 3 % of the land surface but storing up to 30 % of the vast soil carbon stock (Frolking et al., 2011). High-latitude peatlands alone store more than 400 Gt C (Hugelius et al., 2020), and tropical peatland carbon is thought to be more than 100 Gt C (Dargie et al., 2017). This carbon stock has accumulated over millennia – approximately 10 000 years since the Last Glacial Maximum – but can be released very quickly if the peatland becomes dry or otherwise loses its function (Maljanen et al., 2010; Tiemeyer et al., 2016). This has been taking place across the world's peatlands over the last ∼170 years due to land use conversion for agriculture, leading to additional greenhouse gas emissions (Leifeld and Menichetti, 2018). Climate change may also lead to drying or shifts in vegetation that drive carbon loss in currently functional peatlands (Swindles et al., 2019; Dieleman et al., 2015). In addition, peat fires are increasing in severity under climate change (e.g. Scholten et al., 2021). Thus, this carbon stock is both large and vulnerable.

It is therefore vital that we include peatlands in Earth system models (ESMs) that are used to make projections of future climate change, including feedbacks within the global carbon cycle (Loisel et al., 2021). However, none of the models in the recent 6th Coupled Model Intercomparison Project included a representation of peatlands (Arora et al., 2020).

Peatlands can display both vulnerability and resilience via a suite of autogenic feedbacks (Waddington et al., 2015), with self-restoring properties that allow them to persist in conditions where they would not form today but with the potential for rapid carbon losses if they are pushed beyond their resilience threshold. In particular, the physical characteristics of peat can change over time – often in response to changes in the water table or permafrost thaw (see Frolking et al., 2011) – and this in turn influences the hydrological dynamics. Up to a certain point, peat that is more decomposed holds water better. Thus, if a peatland water table drops and peat starts to decompose, the peat that is more decomposed leads to increased water-holding capacity and can bring the water table back up again, leading to resilience. On the other hand, if the peat drainage or decomposition is more severe, it can cross a threshold where it loses the ability to maintain its water table, leading to rapid carbon loss and further degradation of the soil structure. This threshold is shown for example in Wang et al. (2021, Fig. 7), where the soil characteristics change dramatically above a threshold bulk density of 200 kg m⁻³ (0.2 g cm⁻³).

The global land surface schemes that do simulate peatland carbon stocks (Qiu et al., 2018; Bechtold et al., 2019; Müller and Joos, 2020) do not simulate the interplay of processes that leads to the self-sustaining and threshold-type behaviours. Thus, the vulnerability of the carbon stocks in such models cannot be properly simulated. In particular, while modellers have prescribed thermal and hydraulic properties for organic soils (Beringer et al., 2001; Lawrence and Slater, 2008; Chadburn et al., 2015b; Guimberteau et al., 2018), they do not let these parameters vary dynamically as the carbon in the soil changes – for example, an organic layer might decompose substantially during the course of a simulation, and therefore its thermal and hydraulic properties should also change. Occasionally models have simulated such a coupling with a limited set of parameters (Koven et al., 2009), but none have produced a fully coupled version.

Since these dynamics are driven by changes in the vertical structure of the soil organic matter, it is important to resolve the vertical profile of soil carbon (as opposed to a scheme where the soil carbon is treated as a single “box”, e.g. Joint UK Land Environment Simulator (JULES)-CN in Wiltshire et al., 2021). Previous studies have shown that the standard vertically resolved soil carbon scheme in land surface models fails to recreate soil carbon profiles at sites with peat or a thick organic layer (Chadburn et al., 2017). Essentially, the models are not able to accumulate peat on top of the soil column since the soil layers are not allowed to grow or shrink, and thus carbon is continually added to the top soil layer, which contains an unrealistically high carbon content, and the high carbon concentration does not extend far enough into deeper soil layers.

Specialised peat models such as DigiBog and the Holocene Peat Model (HPM) (Baird et al., 2012; Frolking et al., 2010) vertically resolve peatland structure by tracking the carbon that is added each year by treating it as a separate layer added on top of the soil column. This results in a very large number of layers that would be computationally unwieldy for global modelling. It is also only applicable to peatlands and does not provide the functionality to model the continuous transitions between mineral and organic soils (both in time and space).

In this paper we present a new scheme that resolves these issues, allowing vertical accumulation of peat and dynamic coupling between thermal and hydraulic soil properties. This scheme is implemented and demonstrated in the JULES land surface model for northern and temperate sites. However, the new methods and relationships we use in this model can be used to improve other land surface schemes.

2 Model description

2.1 Overview of standard JULES

JULES is the land surface model used in the UK Earth System Model (UKESM) (Sellar et al., 2019). It is a community model that represents the surface energy balance, heat and water fluxes, snowpack dynamics, vegetation dynamics, soil biogeochemistry, and carbon and nitrogen fluxes (Best et al., 2011; Burke et al., 2017a; Clark et al., 2011; Harper et al., 2016; Wiltshire et al., 2021). As well as being used in UKESM, JULES takes part in multimodel analyses such as the Inter‐Sectoral Model Intercomparison Project (Rosenzweig et al., 2017) and the Global Carbon Project (Friedlingstein et al., 2019; Saunois et al., 2020) and has been used to make global projections, for example, of future hydrology, permafrost thaw, and carbon and methane emissions and their climate feedbacks (Burke et al., 2017b; Chadburn et al., 2015b; Comyn-Platt et al., 2018; Gedney et al., 2019)

JULES includes a vertically resolved soil carbon scheme (Burke et al., 2017a), although this has not yet been used in the Earth system model configuration. The scheme is based on Roth-C (Jenkinson, 1990; Jenkinson and Coleman, 1999), with the carbon pools of the Roth-C model simulated separately for each soil layer. Some vertical processes have been added, such as a diffusive mixing, which represents bioturbation and/or cryoturbation (see Burke et al., 2017a, for details). This soil carbon scheme has more recently been coupled to a vertically resolved nitrogen model described in Wiltshire et al. (2021). In this paper we build on this vertically resolved soil carbon–nitrogen scheme in JULES. Note that all the simulations in this paper use the same branch of JULES. We generically refer to any configuration with the new peat functionality enabled as “JULES-Peat”, which is further sub-divided into different simulations.

2.2 Modification to decomposition functions

As part of the development of the new, peat-enabled version of JULES we improved the response of soil carbon decomposition both to soil moisture and nitrogen availability. These changes were made based on well-known principles. Firstly, microbial activity drops to zero in completely dry conditions (Yan et al., 2018). Secondly, respiration in anaerobic conditions is known to be no higher than 20 % of the maximum rate in aerobic conditions (Schuur et al., 2015). Finally, when microbes lack nitrogen, they tend to decompose plant litter faster in order to “mine” for nitrogen (Craine et al., 2007); this is in contrast to the original scheme introduced by Wiltshire et al. (2021) in which the decomposition of litter is inhibited when nitrogen is in short supply.

The decomposition of soil carbon in JULES is calculated as follows: for each soil carbon pool (C_p kg m⁻², where p denotes the pool number), the turnover rate when the nitrogen in the system is not limiting (R_p,pot) is given by

$$\begin{array}{}\text{(1)}& R_{p,\mathrm{pot}}=k_p{C}_p{F}_T\left(T_{\mathrm{soil}}\right)F_{\mathit{\theta }}\left(\mathit{\theta }\right)F_v\left(v\right),\end{array}$$

where k_p is a fixed constant for each pool in s⁻¹ (Clark et al., 2011). The functions of temperature (F_T(T_soil)) and moisture (F_θ(θ)) depend on the temperature (T_soil, K) and moisture content (θ, fraction of saturation) of the soil. The function F_v(v) depends on the vegetation cover fraction (v) (Clark et al., 2011). When the vertically resolved soil carbon scheme is used, there is an additional multiplier, exp (−zτ_resp), where z is depth in the soil and τ_resp represents an additional decay of carbon decomposition rate with depth.

F_θ is a function of the soil moisture. The standard version of JULES uses the following function, which is also shown in green in Fig. 1:

$$\begin{array}{}\text{(2)}& F_{\mathit{\theta }}\left(\mathit{\theta }\right)=\left\{\begin{array}{ll}\mathrm{1}-\mathrm{0.8}(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{o}})& \text{for }\mathit{\theta }>{\mathit{\theta }}_{\mathrm{o}}\\ \mathrm{0.2}+\mathrm{0.8}\left(\frac{\mathit{\theta }-{\mathit{\theta }}_{min}}{{\mathit{\theta }}_{\mathrm{o}}-{\mathit{\theta }}_{min}}\right)& \text{for }{\mathit{\theta }}_{min}<\mathit{\theta }\le {\mathit{\theta }}_{\mathrm{o}}\\ \mathrm{0.2}& \text{for }\mathit{\theta }\le {\mathit{\theta }}_{min}\end{array}\right.,\end{array}$$

where ${\mathit{\theta }}_{\mathrm{o}}=\mathrm{0.5}(\mathrm{1}+{\mathit{\theta }}_{\mathrm{w}})$, θ_w is the wilting point water content as a fraction of saturation and ${\mathit{\theta }}_{min}=\mathrm{1.7}{\mathit{\theta }}_{\mathrm{w}}$. F_θ takes a value between ∼ 0.6–0.85 in saturated conditions; i.e. the decomposition rate in saturated conditions is between 60 %–85 % of its maximum rate. However, in reality, aerobic respiration essentially stops in saturated conditions, and anaerobic respiration takes place instead, with a rate less than 20 % of the maximum aerobic respiration rate (Schuur et al., 2015). The fact that decomposition is suppressed under saturated conditions is key to the formation of peat. Therefore, we modified the decomposition function so that it takes a value of 0.2 when the soil moisture is saturated. We also changed the behaviour of this function under dry conditions, since there are a number of studies available that indicate the shape of this function (Moyano et al., 2012, 2013; Yan et al., 2018), which should increase in a close-to-linear manner from zero decomposition rate at zero soil moisture content.

Figure 1 Original (green) and updated (blue) function of moisture used to determine soil respiration (F_θ; Eq. 3).

In addition, for undecomposed organic soils specifically, critical and wilting point soil moisture can be very small due to the large pore spaces and thus low capillary suction (the critical point can be as low as 10 % saturation). The formulation of θ_o, on the other hand, limits the optimum soil moisture content for respiration to a minimum of 50 % saturation, which can be up 5 × higher than the critical point. The critical point is defined by a capillary suction of 3.36 m. Moyano et al. (2013) show in their Fig. 3b that the respiration response to soil moisture reaches a maximum at around this value. They show a moisture response curve for a high carbon soil in their Fig. 3a. The curve reaches a maximum at around 30 %–40 % saturation and stays at a high value until ∼ 75 % saturation, in contrast to the original formulation in JULES, which reaches its maximum only at the point θ_o. Therefore, for soil layers in which the critical soil moisture is lower than θ_o, we set the soil respiration to reach its maximum at the critical soil moisture content and remain at its maximum value until the original maximum θ_o, resembling the “high C content” curve in Moyano et al. (2013, Fig. 3a). We therefore define a “lower” θ_o, ${\mathit{\theta }}_{\mathrm{o},\mathrm{l}}=\min ({\mathit{\theta }}_{\mathrm{crit}},{\mathit{\theta }}_{\mathrm{o}})$. The old and new functions are shown in Fig. 1.

The new function used in JULES-Peat is therefore

$$\begin{array}{}\text{(3)}& F_{\mathit{\theta }}\left(\mathit{\theta }\right)=\left\{\begin{array}{ll}\mathrm{1}-\mathrm{0.8}\left(\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{o}}}{\mathrm{1}-{\mathit{\theta }}_{\mathrm{o}}}\right)& \text{for }\mathit{\theta }>{\mathit{\theta }}_{\mathrm{o}}\\ \mathrm{1}& \text{for }{\mathit{\theta }}_{\mathrm{o},\mathrm{l}}<\mathit{\theta }\le {\mathit{\theta }}_{\mathrm{o}}\\ \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{o},\mathrm{l}}}& \text{for }\mathrm{0}\le \mathit{\theta }\le {\mathit{\theta }}_{\mathrm{o},\mathrm{l}}\end{array}\right..\end{array}$$

In the standard version of JULES, in situations where nitrogen is limiting, the decomposition of the litter carbon pools (decomposable plant material, DPM; and resistant plant material, RPM) is reduced. This is because the more decomposed pools have a higher nitrogen content – or lower C : N ratio – and therefore to decompose the litter carbon into the BIO (biomass) and HUM (hummus) pools requires a source of nitrogen – and nitrogen is thus “immobilised”. Thus, the decomposition terms R_pot for DPM and RPM pools are multiplied by a factor F_N, which is given by

$$\begin{array}{}\text{(4)}& F_N={\displaystyle \frac{\left(M_{\mathrm{BIO}}+M_{\mathrm{HUM}}-I_{\mathrm{BIO}}-I_{\mathrm{HUM}}\right)\mathrm{\Delta }t+N_{\mathrm{in}}}{\left(D_{\mathrm{DPM}}+D_{\mathrm{RPM}}\right)\mathrm{\Delta }t}},\end{array}$$

where N_in is the total soil inorganic N pool in kg [N] m⁻², M_p and I_p are mineralisation and immobilisation of nitrogen, respectively, from pool p in kg [N] m⁻² s⁻¹, and Δt is the time step. D_DPM and D_RPM are the net demand associated with decomposition of each of the litter pools:

$$\begin{array}{}\text{(5)}& D_p=I_{p,\mathrm{pot}}-M_{p,\mathrm{pot}},\end{array}$$

where p is either RPM or DPM; see Wiltshire et al. (2021) for details.

However, in reality the microbes would continue decomposing the litter pools in order to access the nitrogen for their own survival. They would not be able to transform all of the decomposed carbon into biomass due to lack of nitrogen, but the carbon would decompose and would simply be released to the atmosphere as carbon dioxide, i.e. their carbon use efficiency reduces (Manzoni et al., 2012). Therefore, instead of modifying the litter decomposition rate with the factor F_N, we modify the fraction of decomposed carbon that is released to the atmosphere vs. stored in the soil. This means that the limitation term has to take a slightly different form. The new function is as follows:

$$\begin{array}{}\text{(6)}& F_N={\displaystyle \frac{\begin{array}{c}\left(M_{\mathrm{BIO}}+M_{\mathrm{HUM}}-I_{\mathrm{BIO}}-I_{\mathrm{HUM}}+M_{\mathrm{DPM}}\right.\\ \left.+M_{\mathrm{RPM}}\right)\mathrm{\Delta }t+N_{\mathrm{in}}\end{array}}{\left(I_{\mathrm{DPM}}+I_{\mathrm{RPM}}\right)\mathrm{\Delta }t}},\end{array}$$

The fraction of decomposed carbon that stays in the soil (rather than being released to the atmosphere), in other words the carbon use efficiency, is then multiplied by F_N for the DPM and RPM pools. While nitrogen was not a focus of this study, the need for this modification became apparent once the soil column was allowed to expand with addition of plant litter. This led to an unrealistic positive feedback in which litter carbon was not decomposed due to lack of nitrogen availability, meaning that as litter was added to the layer it took up an ever larger volume, eventually pushing the more nitrogen-rich pools out of the layer completely (further down the column), resulting in zero nitrogen availability and forming unrealistically thick litter layers with no turnover.

2.3 Change of soil column height

Chadburn et al. (2017) showed that the typical soil profile simulated by ESM land surface schemes with vertically resolved soil carbon (JULES and ORCHIDEE) displays a smooth decline with depth that resembles a mineral soil profile. However, in highly organic soils the soil carbon concentration typically increases with depth to a certain point before beginning to decline (Harden et al., 2012). This is because the density of the organic material in the surface is usually lower than in the deeper soil, so there is simply less material altogether in the surface layers, and therefore less carbon. The organic material in deeper layers has a higher density: in part because it becomes compressed by soil/water above it, and in a large part because it is generally more decomposed.

The crucial missing factor in global models (e.g. JULES, ORCHIDEE, CLM) is that the models do not account for the volume that is added to the soil when organic material is added via plant litter, or (conversely) the reduction in volume when organic material decomposes. This means that as well as being unable to simulate the typical profile of a peatland (soil carbon increasing with depth near the surface), unrealistically high carbon contents in surface layers are often simulated; see the original JULES version, shown as red lines in Fig. 5.

Figure 2 Functions used in JULES-Peat. Vertical dashed lines show the range of data that were used to fit the functions. These correspond to the minimum and maximum bulk densities for organic material that we derived for use in JULES (ρ_dpmrpm and ρ_biohum; Sect. 2.3). The additional literature data for saturated hydraulic suction and the Clapp–Hornberger exponent shown in purple were derived from the following papers: Londra (2010), Rezanezhad et al. (2012), Da Silva et al. (1993), Weiss et al. (1998), Päivänen (1973), Boelter (1964), Rydén et al. (1980) and Schwärzel et al. (2006).

In JULES-Peat, the profile of litter inputs into each soil layer and decomposition of soil carbon in the layer is calculated as in Burke et al. (2017a), except that the modified decomposition function is used (Sect. 2.2). When these increments come to be applied to the soil carbon profile, however, the thickness of the soil layers is now recalculated based on the volume of organic matter added or removed. We calculate the change in layer thickness by prescribing a density to each carbon pool, using a higher density for the decomposed carbon pools than the litter carbon pools. After addition or removal of carbon in a given time step, the new effective thickness dz_eff,n of soil layer n relative to the initial layer thickness dz_n is given by

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathrm{d}{z}_{\mathrm{eff},n}& =\mathrm{d}{z}_n+\left(\mathrm{d}{C}_{\mathrm{dpm},n}+\mathrm{d}{C}_{\mathrm{rpm},n}\right)/\left(f_{\mathrm{c}}{\mathit{\rho }}_{\mathrm{dpmrpm}}\right)\\ & +\left(\mathrm{d}{C}_{\mathrm{bio},n}+\mathrm{d}{C}_{\mathrm{hum},n}\right)/\left(f_{\mathrm{c}}{\mathit{\rho }}_{\mathrm{biohum}}\right),\end{array}\end{array}$$

where ρ_dpmrpm and ρ_biohum are the bulk densities associated with the carbon pools (in kg m⁻³), f_c is the fraction of organic material that is carbon and dC_i,n is the increment in carbon pool i in soil layer n. We picked the density of the litter pools, ρ_dpmrpm, to be the lowest density that is typically measured for peat (where DPM and RPM are the two litter carbon pools in JULES), and the density of the more decomposed carbon pools (ρ_biohum) to be the highest density that is typically observed for peat (where BIO and HUM are the more decomposed carbon pools in JULES). Thus, the bulk density of organic material in any given soil layer will fall somewhere between these two extreme values given that each layer typically contains all four carbon pools (albeit in different ratios). The values we chose were 35 kg m⁻³ as the minimum, ρ_dpmrpm, and 210 kg m⁻³ as the maximum, ρ_biohum. These match well with commonly quoted literature values (e.g. Chambers et al., 2010) and were derived from the 5th and 95th percentiles of the bulk densities in the global peat core dataset described in Sect. 3. These limits are shown by vertical dashed lines in Fig. 2. The maximum bulk density of 210 kg m⁻³ corresponds well to the threshold bulk density for peat functioning in, e.g. Wang et al. (2021). We relate the bulk density of the organic material to the carbon content by assuming that f_c=0.56 or, in other words, that 56 % of the organic matter is carbon, which was also based on the 95th percentile of the percentage of carbon in the peat core dataset from Gallego-Sala et al. (2018) (see Sect. 3) and is consistent with the range of observations, e.g. in Chambers et al. (2010).

The new layer thicknesses are labelled as “effective” layer thicknesses (dz_eff). In order to avoid technical difficulties and potential numerical problems with variable soil layer thickness (e.g. if surface layers become very thick), the soil carbon profile is then interpolated back onto the original soil layers.

In order to interpolate the carbon profile, the carbon quantities (C+dC; kg m⁻²) are first transformed to carbon densities (Cden, in kg m⁻³) by dividing them by the layer thicknesses, dz.

Following this, the interpolation of the effective carbon density on the effective layers back into the original layers depends on whether the centre of the original layer is above or below the centre of the effective layer and is calculated as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\mathrm{Cden}}_{n,i}& ={\mathrm{Cden}}_{\mathrm{eff},n,i}+{\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}\left({\mathrm{Cden}}_{\mathrm{eff},n+\mathrm{1},i}-{\mathrm{Cden}}_{\mathrm{eff},n,i}\right)\\ & \times {\displaystyle \frac{z_{\mathrm{eff},n}-z_n}{z_{\mathrm{eff},n}-z_{\mathrm{eff},n+\mathrm{1}}}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left({\mathrm{Cden}}_{\mathrm{eff},n-\mathrm{1},i}-{\mathrm{Cden}}_{\mathrm{eff},n,i}\right)\\ & \times {\displaystyle \frac{z_{\mathrm{eff},n}-z_n}{z_{\mathrm{eff},n}-z_{\mathrm{eff},n-\mathrm{1}}}},\end{array}\end{array}$$

when the original soil layer depth z_n is deeper than the effective soil layer depth z_eff,n, and

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mathrm{Cden}}_{n,i}& ={\mathrm{Cden}}_{\mathrm{eff},n,i}+{\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}\left({\mathrm{Cden}}_{\mathrm{eff},n-\mathrm{1},i}-{\mathrm{Cden}}_{\mathrm{eff},n,i}\right)\\ & \times {\displaystyle \frac{z_{\mathrm{eff},n}-z_n}{z_{\mathrm{eff},n}-z_{\mathrm{eff},n-\mathrm{1}}}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left({\mathrm{Cden}}_{\mathrm{eff},n+\mathrm{1},i}-{\mathrm{Cden}}_{\mathrm{eff},n,i}\right)\\ & \times {\displaystyle \frac{z_{\mathrm{eff},n}-z_n}{z_{\mathrm{eff},n}-z_{\mathrm{eff},n+\mathrm{1}}}},\end{array}\end{array}$$

when the original soil layer depth z_n is shallower than the effective soil layer depth z_eff,n, where z is the centre of each soil layer and i indicates the carbon pool (DPM, RPM, BIO or HUM).

Mathematically, this represents an approximated second-order Taylor expansion of the function Cdens(z) around the point z_eff but with a particular choice regarding the second-order derivative. In order to preserve the vertical structure of the soil, the second-order derivative is assumed to be around z+δz and z−δz, and thus if the gradient of the function changes sharply it will not be smoothed out. This means that a peat layer will not end up being numerically smeared into the rest of the profile. This is explained in detail with equations in the Supplement. Briefly, we used a simple test model where soil carbon inputs and outputs are given prescribed input and turnover rates, we account for the expansion and contraction of the soil column when carbon is added or removed and tested the method of interpolating back onto the original soil layers (i.e. as used in JULES). We ran this simple (and thus much quicker to run) script with very high resolution soil layers to see what the “true” solution for the soil carbon profile would be. We therefore confirmed that our choice of second-order derivative gave the best approximation of the true solution when a lower-resolution soil was used. For details, see the Supplement. Figure S2 in the Supplement shows that in the chosen scheme there does still appear to be some “smearing” in the deeper layers, which are thicker, but using a smaller interval δz^′ leads to numerical instability in the thinner surface layers. Thus, for future development a scheme where δz^′ depends on the soil layer thickness could be considered.

While the carbon profile is interpolated onto the original soil layers in order to keep the layer thicknesses constant, the thickness of the deepest soil layer is updated in order to track the overall change of soil column height. The minimum thickness for that layer is taken from the soil layer thicknesses specified at runtime, and this corresponds to the layer thickness when there is no carbon in the soil. This base layer is extended based on the total volume of the carbon pools in the soil column, and this extension is considered to be the surface elevation. However, the extra thickness of the bottom layer is neglected when calculating fluxes of heat and water and only applies when calculating carbon and nitrogen stocks and fluxes (which are conserved during layer adjustments). We decided that the complexity of modifying the heat and water calculations to account for the variable-thickness base layer was not worth the added complexity, given that the fluxes at the base of the soil column (8 m depth in this study) are generally very small.

2.4 Simulating the age profile

Peat age was simulated following a similar method to Burke et al. (2017a), where the fraction of old carbon is traced throughout the simulation. During each update of the soil carbon pools, the age of each carbon pool is tracked and the weighted average of the soil age is taken for each carbon pool in each soil layer.

Each soil carbon pool in each soil layer is assigned an age, A, at the start of the simulation, which currently is either zero on initialising the spin-up or can be initialised from an existing simulation. Each time the carbon pools are updated, the age of each soil carbon pool in each layer is increased by the time step length. These values are then modified as the soil carbon pools are updated, either due to input of fresh carbon from litter (which has an age of zero and therefore reduces the age of the soil pool), due to mixing of carbon between two layers in which the ages are different, or due to input of carbon into BIO and HUM pools from other pools via decomposition. The general formula to update the age (A) for carbon pool C_i (kg m⁻²), with an increment of carbon $C_i\to C_i+\mathrm{\Delta }{C}_i$, is as follows:

$$\begin{array}{}\text{(10)}& A_{C_i}\to {\displaystyle \frac{A_{C_i}{C}_i+A_{\mathrm{\Delta }{C}_i}\mathrm{\Delta }{C}_i}{C_i+\mathrm{\Delta }{C}_i}}.\end{array}$$

ΔC_i includes both incoming and outgoing fluxes from the pool. For the outgoing fluxes in ΔC_i, we assume that A is the same as for the C_i pool. For an incoming litter flux we assume that A is zero, and incoming fluxes from other pools naturally take the age value from the corresponding pool.

If the soil height accumulation is switched on, the age then must be interpolated back onto the original soil layers as described in Sect. 2.3. We use the same interpolation method as for the soil carbon, which is given in Eqs. (8) and (9).

2.5 Coupling between properties and C concentration

In order to dynamically update the soil physical characteristics, we assume that the physical properties of the organic material in the soil are a function of its bulk density. The bulk density that we simulate in JULES-Peat depends on how decomposed the soil carbon is. More highly decomposed organic matter has a higher bulk density, and its properties change as it decomposes. Notably, the hydraulic conductivity becomes much lower as bulk density (or decomposition) increases, which is included in other peat models (Frolking et al., 2010; Young et al., 2017), but we also fitted relationships between the bulk density and other key physical characteristics, namely porosity, saturated hydraulic suction, the Clapp–Hornberger exponent and thermal conductivity. For the heat capacity we assumed that the heat capacity of the organic material does not significantly change with decomposition status, and therefore we used (1 − porosity, i.e. what fraction of the organic material is solids) multiplied by the heat capacity value of solids of 2.5×10⁶ J K⁻¹ m⁻³, which we took from Beringer et al. (2001).

Since different carbon pools have different bulk densities (see Sect. 2.3), we first calculate the bulk density of the combined organic material in each soil layer, i.e.

$$\begin{array}{}\text{(11)}& {\mathit{\rho }}_{\mathrm{org},n}={\displaystyle \frac{\mathrm{1}}{f_{\mathrm{c}}}}\sum _i{C}_{n,i}/\left(f_{\mathrm{org},n}\mathrm{d}{z}_n\right),\end{array}$$

where f_org,n is the volumetric fraction of organic matter in the soil layer given by

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{f}_{\mathrm{org},n}& ={\displaystyle \frac{\mathrm{1}}{f_{\mathrm{c}}}}\left((C_{n,\mathrm{1}}+C_{n,\mathrm{2}})/{\mathit{\rho }}_{\mathrm{dpmrpm}}\right.\\ & \left.+(C_{n,\mathrm{3}}+C_{n,\mathrm{4}})/{\mathit{\rho }}_{\mathrm{biohum}}\right)/\mathrm{d}{z}_n.\end{array}\end{array}$$

Recent studies have shown that bulk density of peat shows strong relationships with its thermal and hydraulic properties (Liu and Lennartz, 2019; O'Connor et al., 2020). We combined data from these recent syntheses with additional values from the literature in order to get the best estimate of the relationships, which we show in Fig. 2. We fitted the relationships between bulk density and the other physical characteristics of peat using this combined dataset. Fitting was done using orthogonal least squares after normalising the data so that both variables being fitted had the same range of values. For the saturated hydraulic conductivity, the two available datasets showed markedly different relationships (see Fig. 2b), and thus we did not combine these but instead used only the data from Liu and Lennartz (2019) since this was a global synthesis as opposed to the O'Connor et al. (2020) data which were from a single region. The data in Liu and Lennartz (2019) also agreed better with other data, such as Wang et al. (2021), and the original values used for organic soils in JULES, originally given in Dankers et al. (2011) (also shown in Fig. 2). In addition, the fit for porosity was forced to pass through 1 at a bulk density of zero as a physical constraint (this was achieved by modifying the normalisation factor until the intercept of the fit was exactly 1). The additional literature data for saturated hydraulic suction and the Clapp–Hornberger exponent shown in purple in Fig. 2 were derived from the following papers: Londra (2010), Rezanezhad et al. (2012), Da Silva et al. (1993), Weiss et al. (1998), Päivänen (1973), Boelter (1964), Rydén et al. (1980) and Schwärzel et al. (2006).

Specifically, we relate the following soil properties to bulk density:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{sat},\mathrm{org}}=\exp \left(\mathrm{0.023}{\mathit{\rho }}_{\mathrm{org}}-\mathrm{5.08}\right),\text{(14)}& {\displaystyle }{b}_{\mathrm{org}}=\mathrm{0.0304}{\mathit{\rho }}_{\mathrm{org}}+\mathrm{1.53},\text{(15)}& {\displaystyle }{K}_{\mathrm{sat},\mathrm{org}}=\exp \left(-\mathrm{0.0532}{\mathit{\rho }}_{\mathrm{org}}-\mathrm{6.63}\right),\text{(16)}& {\displaystyle }{\mathit{\theta }}_{\mathrm{sat},\mathrm{org}}=\mathrm{1}-{\mathit{\rho }}_{\mathrm{org}}/\mathrm{1260},\text{(17)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{org}}=\mathrm{0.06},\text{(18)}& {\displaystyle }{\mathrm{hcap}}_{\mathrm{org}}=\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{6}}{\mathit{\rho }}_{\mathrm{org}}/\mathrm{1260},\end{array}$$

where Ψ_sat,org is soil matric suction at saturation (m), b_org is the Clapp–Hornberger exponent (unitless), K_sat,org is hydraulic conductivity (in units of m s⁻¹; note that JULES uses kg m⁻² s⁻¹ so this is multiplied by 1000 for use in JULES), θ_sat,org is the volumetric soil moisture at saturation (m³ m⁻³), λ_org is the dry thermal conductivity (W m⁻¹ K⁻¹) and hcap_org is the dry heat capacity (J m⁻³ K⁻¹). If a soil layer is not 100 % organic then we combine these calculated organic parameters with the properties of the underlying mineral soil, following Chadburn et al. (2015a). The remaining hydraulic parameters θ_wilt and θ_crit, which are the volumetric soil moisture at the wilting point and critical point, respectively (defined in terms of hydraulic suction), are functions of the other parameters and are recalculated when the other parameters are updated (see, e.g. Chadburn et al., 2015a).

3 Simulations and evaluation data

We used a large suite of simulations at 24 sites that have been use for JULES development and evaluation in Chadburn et al. (2017, 2020), Nakhavali et al. (2018), Smith et al. (2021) and Gao et al. (2022), along with Scotty Creek (Helbig et al., 2016, 2017a, b), Pleistocene Park (Euskirchen et al., 2017b), Imnavait (Euskirchen et al., 2017a), and Eight Mile Lake (Celis et al., 2019). We were able to include the four new sites that have not been used in previous studies with JULES due to additional data becoming available. The sites are fairly evenly distributed between tundra, boreal and temperate climate zones; see Table 1. Some of the sites, namely Abisko, Seida and Imnavait, provided data from different landscape types, resulting in 29 simulations in total. The climate forcing data were prepared as described in Chadburn et al. (2017) using WFD and WFDEI (Weedon et al., 2011, 2014) corrected with local climate data from the sites and covers the period 1901–2018 inclusive. The simulations were spun up for 10 000 years using repeated climate forcing data from 1901–1910.

Jammet et al. (2017)Chadburn et al. (2017)Jammet et al. (2017)Chadburn et al. (2020)Drewer et al. (2010)Gao et al. (2022)Gielen et al. (2010, 2011)Nakhavali et al. (2018)Long et al. (2010)Flanagan and Syed (2011)Gao et al. (2022)Walmsley et al. (2011)Nakhavali et al. (2018)Kittler et al. (2017)Göckede et al. (2019)Chadburn et al. (2020)Nilsson et al. (2008)Sagerfors et al. (2008)Gao et al. (2022)Celis et al. (2019)Kutsch et al. (2010)Schrumpf et al. (2011)Nakhavali et al. (2018)Euskirchen et al. (2017a)Kjellman et al. (2018)Smith et al. (2021)Fortuniak et al. (2021)Gao et al. (2022)Van der Molen et al. (2007)Parmentier et al. (2011)Chadburn et al. (2017)Aurela et al. (2009)Lohila et al. (2010)Chadburn et al. (2020)Moore et al. (2011)Brown et al. (2014)Gao et al. (2022)Euskirchen et al. (2017b)Boike et al. (2019)Chadburn et al. (2015a)Helbig et al. (2016, 2017a, b)Marushchak et al. (2013)Biasi et al. (2014)Chadburn et al. (2020)Zhang et al. (2020)Gao et al. (2022)Boike et al. (2018)Chadburn et al. (2017)Peichl and Arain (2006)Peichl et al. (2010)Nakhavali et al. (2018)Valach et al. (2021)Miller et al. (2008)Miller and Fujii (2010)Chadburn et al. (2020)Elberling et al. (2008)Chadburn et al. (2017)

Table 1 Sites used in the suite of JULES simulations. References are both for site data and for simulations of these sites with JULES.

Table 2 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 spin-up and on during the main run. The “Decomp. function” refers to changing from the original to the new decomposition function shown in Fig. 1, and we also changed a switch that was using the total soil moisture instead of the unfrozen soil moisture, and thus the new decomposition function is a function of the unfrozen soil moisture (which is more realistic since frozen water is not available for microbes to use). Where the Initial C is given as “Peat”, we initialise all spin-ups with the spun up profile for Auchencorth from JULES-Peat-W, with ∼ 1.5 m of peat, and otherwise initialise the model with zero soil carbon.

The initial JULES simulation (JULES, Table 2) is based on the configuration in Chadburn et al. (2020) but now has 20 soil layers extending to around 7.9 m, with thicknesses given in Table S3. This was originally derived from the JULES-ES configuration (see https://jules.jchmr.org/content/core-configurations, last access: 16 July 2021), with extra processes added to enhance the simulation, particularly for high latitudes. For example, an extra plant functional type (PFT) is included to represent Arctic grass, based on C₃ grass with temperature optimum adjusted to grow in colder climates, layered soil carbon and nitrogen are simulated, and a bedrock column is included below the soil to simulate heat conduction. Starting from this baseline simulation, we then activated the new processes in JULES-Peat described in Sect. 2. See the Code and Data Availability section for the full configuration.

For most of the simulations the standard TOPMODEL-based large-scale hydrology scheme was used, which calculates the lateral flow of water from each grid cell based on the topographic index information of the grid cell (Gedney and Cox, 2003). In order to simulate a wetter site, for example a topographically controlled peatland, which would essentially be a wetter fraction of the grid cell than the grid cell average, we simply set the lateral flow to zero (JULES-Peat-W and JULES-Peat-W10, Table 2). Neither of these hydrological scenarios are necessarily expected to be realistic for the sites. The aim was to test the response of the model to wetter vs. drier conditions. How to simulate peatland hydrology realistically is a challenge and will be addressed in future work (see also Bechtold et al., 2019). In JULES-Peat-W-drain and JULES-Peat-W10-drain (Table 2), the lateral flow is set to zero during spin-up but switched back on during the main run to approximate drainage.

The list of simulations is shown in Table 2. JULES-Peat-B and JULES-Peat-B10 are baseline JULES-Peat simulations. The impact of switching off the lateral flow (JULES-Peat-W and JULES-Peat-W10) and initialising the sites with peat (JULES-Peat-i and JULES-Peat-i10) – instead of letting the carbon build up from zero – was then tested. Initialising with peat soil tests whether the model retains a peat layer at sites where peat was not able to form from scratch, which would indicate self-regulating functions of peat in JULES-Peat. To initialise with peat we used the spun-up profile from Auchencorth JULES-Peat-W simulation, since this simulation had formed a thick, 1.7 m layer of peat (see Fig. S12). The only site that formed thicker peat in JULES-Peat was CA_WP1, which formed an extremely thick (5–6 m) peat layer (Table S4, Fig. S12), and thus we chose Auchencorth as a site with a thick but not extreme peat profile. In addition to these simulations, new processes were firstly switched on one by one in factorial until the full “JULES-Peat” was run – these simulations are shown in the Supplement. During this process a few different parameter combinations were tested to make sure the soil carbon profile and the age–depth profile looked realistic. In particular we altered the rate of decay of soil respiration with depth (efolding), τ_resp, and the efolding depth of litter inputs to the soil, τ_lit (called ξ_lit in Wiltshire et al., 2021). A higher value of τ_lit means more of the litter is added to the surface compared to deeper in the soil. We ran JULES-Peat with two different values of τ_lit, 5 and 10, where 5 is the original value in JULES and any simulation with a “10” on the end (Table 2) has τ_lit=10. Note that we did not develop the vegetation module further here, and this should be addressed in future work (Sect. 5.3).

For evaluation we used a globally distributed dataset of peat profiles (Gallego-Sala et al., 2018). We divided these data into major climate zones and selected only those zones that are covered by the JULES simulations: temperate, boreal and tundra. This left 216 sites: 12 tundra, 127 boreal and 77 temperate. We compared these peat profiles against the site simulations divided into the same climate zones. We also used this dataset to derive the values for maximum and minimum bulk density (ρ_dpmrpm and ρ_biohum; see Sect. 2.3) and the fraction of carbon in organic matter (f_c). For this calculation we selected only the data points that were organic rich with minimal mineral content, for which the percentage of carbon (by mass) was higher than 30 %. This left over 24 500 data points where the vast majority of the soil by volume is organic material.

Before comparing the peat core dataset against the JULES simulations, we only remove values where the percent carbon is less than 15 % as we take this as a common definition of organic soil (United States Department of Agriculture, 1975), which can include some mineral material. We use the same definition to assess where JULES-Peat simulates a peat or a mineral soil, noting that to estimate the percent carbon by mass in JULES, we assume the mineral fraction of the soil has a bulk density of 1500 kg m⁻³ (Hossain et al., 2015). We isolate the peat layers from the JULES simulations in order to evaluate comparable soil layers against the observed peat core dataset and select only sites where JULES simulates peat in at least the top four soil layers (i.e. to 39 cm depth, since this is the closest layer to the 40 cm depth specified for defining organic soils; United States Department of Agriculture, 1975).

For further evaluation data we used individual soil carbon profiles from other data sources, which are available for some of the sites in Table 1 (references given in the table).

4 Results

4.1 Representation of mineral soils

Since JULES is a global model, it is important that adding the functionality to represent peat does not degrade model performance for mineral soils. Therefore, we first evaluate the model at mineral soil sites.

Figure 3 Evaluation at mineral sites. Observations are shown with black crosses or lines. The black line in (h) is the median of all sites, with the grey area indicating the interquartile range. Note that the axis ranges are different in (h). The simulations JULES-Peat-B and JULES-Peat-B10 have different values of τ_lit. For information about the sites refer to Table 1, and for details of the JULES simulations refer to Table 2.

In soils where the organic material is a relatively small fraction of the total soil, the original soil carbon scheme is able to perform well, since the expansion of the soil column due to the addition of organic material will be relatively small. Figure 3 shows the model performance at mineral soil sites where measured soil carbon profiles were available. While the individual sites are not well simulated, the general form of the profiles – resembling an exponential decline with depth – are recreated reasonably well in the standard JULES configuration.

Table 3 RMSEs of JULES simulations against various observations. “Temp.” is short for temperate. The final two columns refer to the median of the site-specific observations shown in Figs. 3 and 4 respectively. The best-performing simulation for each column is highlighted in bold.​​​​​​​

Figure 3 also shows two JULES-Peat simulations in light green and dark green (JULES-Peat-B and JULES-Peat-B10, respectively, Table 2). As discussed in Sect. 3, these have different values of τ_lit and are both plausible. In general, the new model version is also able to simulate a profile that resembles a mineral soil despite forming peaty profiles at a few of the sites, especially Hainich and Carlow (where the carbon is overestimated by all configurations of JULES, Fig. 3e–f). Aggregated across all sites, the updated model versions produce a profile with somewhat lower carbon density at the surface compared to standard JULES and less of a decline in carbon with depth (Fig. 3h). The lower carbon density at the surface matches better with observations than the original JULES simulation, but the carbon at depth tends to be overestimated. In terms of root-mean-squared error (RMSE), the aggregated profile is improved in JULES-Peat-B and JULES-Peat-B10 (RMSE 6.2 and 8.1 kg m⁻³) compared to JULES (RMSE 12.1 kg m⁻³); Table 3. Overall we conclude that despite poor model performance at individual sites, the aggregated soil carbon profiles in both JULES and JULES-Peat adequately resemble observed mineral soil profiles (Fig. 3h).

4.2 Model evaluation at peatland sites

We assess the performance of the full JULES-Peat model configuration at the selection of simulated sites for which observed soil carbon profiles are available and organic soils are present (Fig. 4). A few of the individual sites are well simulated, and almost all sites are simulated significantly better in all of the JULES-Peat configurations than they are in JULES (RMSE of median profile 12.9–13.5 kg m⁻³ for JULES-Peat configurations shown and 25.5 kg m⁻³ for JULES; Table 3). Note that two additional JULES-Peat simulations are shown in Fig. 4, JULES-Peat-W and JULES-Peat-W10 (dark blue and light blue lines, respectively), where the lateral water flow is set to zero since we expect that this would lead to a wetter soil and a more realistic simulation of a topographically controlled peatland.

Figure 4 Evaluation at organic sites. Observations are shown with black crosses or lines. The black line in (i) is the median of all sites, with the grey area indicating the interquartile range. Note that the axis ranges are different in (i). In JULES-Peat-W and JULES-PeatW10 (dark blue and light blue lines, respectively) the lateral water flow is set to zero. The simulations with “10” on the end have τ_lit=10. For information about the sites refer to Table 1, and for details of the JULES simulations refer to Table 2.

Figure 5 Evaluation of JULES-Peat carbon profile against the median profile from peat cores. The grey-shaded area shows the interquartile range. Only sites where the model simulates a sufficiently thick organic layer to classify it as peat (> 15 % carbon by mass for ≥ 39 cm; see Sect. 3) are shown. The sites that are shown for each simulation are highlighted in bold in Table S4.

Figure 6 Evaluation of JULES-Peat age–depth profile against the median profile from peat cores. The grey-shaded area shows the interquartile range. Only sites where the model simulates a sufficiently thick organic layer to classify it as peat (> 15 % carbon by mass for ≥ 39 cm, see Sect. 3) are shown. The sites that are shown for each simulation are highlighted in bold in Table S4.

We then evaluate the soil carbon and age–depth profiles in JULES and JULES-Peat against the global dataset of peat cores (described in Sect. 3) (Gallego-Sala et al., 2018) in Figs. 5 and 6. These figures show simulated soil carbon profiles at sites where the model simulates a carbon percentage by mass of more than 15 % for at least 39 cm (see Sect. 3) and compares these against the median of the equivalent data (percent carbon > 15 %) from the global dataset of peat cores (Gallego-Sala et al., 2018). It is clear that at these peaty sites, the original JULES configuration simulates a carbon density that is too high in the surface layers and too low in the deeper soil (red lines in Fig. 5).

In order to quantify the total improvement in the various JULES-Peat simulations compared with the original JULES configuration we take the RMSE between the median soil carbon profiles for each climate zone, shown in Table 3. The best-performing version (shown in bold) has a RMSE that is reduced by 35 % for temperate peatland sites (from 23.0 to 14.9 kg m⁻³, Table 3) and by almost 80 % for boreal peatland sites (RMSE reducing from 37.3 to 7.7 kg m⁻³, Table 3). We see that the age at the soil surface was typically too high in the original configuration of JULES (Fig. 6, red lines). In JULES-Peat, the age at the soil surface is better simulated, age throughout the profile is generally realistic (solid blue and green lines in Fig. 6, mostly falling within the interquartile range of the observations), and RMSE in age is reduced by 32 % (temperate sites) and 56 % (boreal sites) for the configurations that perform best in terms of carbon profile (JULES-Peat-W and JULES-Peat-B10 respectively, Table 3).

Figure 7 Evaluation of JULES-Peat bulk density profile against median profile from peat cores. The grey-shaded area shows the interquartile range. Only sites where the model simulates a sufficiently thick organic layer to classify it as peat (> 15 % carbon by mass for ≥ 39 cm, see Sect. 3) are shown. The sites that are shown for each simulation are highlighted in bold in Table S4. Model profiles are cut off because only the organic layers are plotted.

We also evaluated the bulk density profiles against the same peat core dataset, and the results are shown in Fig. 7. The bulk density in JULES-Peat tends to start at realistic values at the surface but increases too quickly with depth. It is remarkable how little the observed bulk density at boreal and temperate peatland sites varies with depth compared with the tundra sites (Fig. 7), although this may be related to the larger sample size of boreal and temperate sites (127 and 77 vs. 12 for tundra) leading to a more “smoothed” profile.

4.3 Impact of environmental and initial conditions on peat profiles

The best simulations in JULES-Peat are generally those where drainage is impeded to make the soil wetter (JULES-Peat-W; see the bold numbers in Table 2). In particular, JULES-Peat-W simulates a more realistic carbon density profile than the other model configurations for the temperate peatland sites (see the dark blue lines in Fig. 5). For the boreal sites, JULES-Peat-W10 accurately captures the gradient of the soil carbon profile in the top 50 cm of soil (light blue lines in Fig. 5). The age–depth profiles in JULES-Peat-W and JULES-Peat-W10 also correspond most closely to the median measured age–depth profiles. Since peat generally forms in wetter places, the fact that the simulations without lateral water flow out of the soil (JULES-Peat-W and JULES-Peat-W10) compare best against observations from peat soils is expected if the model realistically forms more peat in wetter soils.

In simulations whose name ends with “10” (e.g. JULES-Peat-W10, shown in light blue) the distribution of litter inputs into the soil is more weighted towards the surface (τ_lit=10 as opposed to 5, Table 2). We generally find that using the lower value of τ_lit matches better with the data for the temperate peatlands and that the higher value is better for boreal peatlands (see the RMSE values for carbon profile in Table 3). This suggests that τ_lit should depend on plant functional type, which it would in reality (shallower-rooting plants would deposit more of their litter nearer the surface) and implies that smaller or shallower-rooting PFTs should grow in colder regions, as is indeed the case (Jackson et al., 1996).

JULES-Peat was sometimes able to accumulate peat starting with no carbon in the soil at time = 0 at one of the tundra sites (Iskoras) in selected configurations (JULES-Peat-B and JULES-Peat-W, see Fig. S12). These simulations form a relatively thin organic layer at this site. However, when the simulations were initiated with an existing peat profile instead of zero soil carbon in JULES-Peat-i (see Table 2), the peat profile was maintained at Iskoras and continued to accumulate (Fig. S12 and Table S4). This then forms a realistic carbon density profile and a reasonably realistic age–depth profile, shown by dashed cyan lines in Figs. 5a and 6a. There are multiple reasons why the model may not accumulate much peat at tundra sites, including a lack of representation of more favourable paleoclimate conditions during spin-up and the simulation of soils that are too dry (Smith et al., 2021). Nonetheless, the Iskoras simulation by JULES-Peat-i indicates that when the model does simulate peat at a tundra site it can form a realistic profile (Fig. 5a, dashed cyan line). Similarly, the bulk density profile simulated for the tundra sites is realistic for the JULES-Peat-i simulation (dashed cyan line in Fig. 7).

4.4 New processes in JULES-Peat

We initially introduced the additional processes and parameter changes that were incorporated in JULES-Peat into simulations one by one to test each process before running the full model. These simulations are shown in the Supplement (Sect. S2). The key development that allows the shape of the carbon profiles to be more realistic is accounting for the volume of organic material added and removed from the soil column (Eq. 7) as described in Sect. 2.3. In particular, this process enables more carbon to reach the deeper soil and makes the carbon density in the surface lower since it expands when plant litter (which has low density; ρ_dpmrpm=35 kg m⁻³) is added. These differences are clear in all of the JULES-acc and JULES-Peat simulations in comparison to the original JULES (Figs. S3 and 5). The majority of the reduction in RMSE is already achieved by adding in this process alone (reduced from 23.0 to 10.4 kg m⁻³ at temperate sites and 37.3 to 17.3 kg m⁻³ for boreal sites; Table S2).

While it does not significantly reduce the RMSE by itself, modifying the moisture function to suppress decomposition when saturated (Sect. 2.2, Eq. 3) allows more peat to form in wetter areas, which is a crucial factor in simulating realistic peatland distribution and future dynamics. Reducing drainage makes the simulation worse for mineral soil sites, which is exactly what would be expected (wetter soil ⇒ more peat forms), and it almost universally improves carbon profiles for organic sites (compare JULES-Peat-B and JULES-Peat-B10 against JULES-Peat-W and JULES-Peat-W10 in Table 3).

The other major new process introduced is that JULES-Peat simulates its own soil characteristics (Eqs. 13–18, Sect. 2.5). While we do not have measured profiles of soil characteristics to compare against, we can compare the soil thermal and hydraulic parameters simulated by JULES-Peat against those prescribed in JULES. Comparisons are shown in the Supplement for key parameters: K_sat, Ψ_sat, θ_sat and b for the simulations JULES-Peat-B, JULES-Peat-B10 and JULES-Peat-i10 (Figs. S7–S10). For some sites there is a good correspondence between the simulated and prescribed parameters, while for others there are significant differences, but all simulated profiles behave sensibly.

Simulating these soil properties dynamically leads, in many cases, to a thicker organic layer (compare JULES-accC and JULES-accC10 with JULES-Peat-B and JULES-Peat-B10 in Table S4) and more soil organic carbon (Figs. S5 and S6). This increase in carbon results in profiles that are significantly more realistic for organic sites (Fig. S6; RMSE reduced from ∼ 18–20 to 13–14 kg m⁻³) and marginally less realistic for mineral soil sites (Fig. S5; RMSE increased from ∼ 6–7 to 6–8 kg m⁻³). This suggests a self-reinforcing feedback of peat accumulation with soil characteristics, e.g. if peat accumulation has started, it is more likely to continue, and can lead to various important dynamics, which are discussed in the following sections.

5 Discussion

5.1 Drainage, subsidence and feedbacks between hydrology and soil carbon

For the simulations where the lateral flow of water was switched off during spin-up in order to simulate a wetter and more “peaty” soil (JULES-Peat-W and JULES-Peat-W10), we ran an alternative realisation of the 20th century where the lateral flow was switched back on (JULES-Peat-W-drain and JULES-Peat-W10-drain). At many sites the lateral flow is negligible in any case and this does not significantly affect the results. However, at some sites the “wet” (-W) simulation maintains a water table near the surface, whereas the baseline simulation (JULES-Peat-B/JULES-Peat-B10) does not, and the drained simulations therefore experience a major drop in water table and a subsequent degradation of the peat and a drop in the surface elevation (subsidence). Figure 8 shows the four sites for which the change in water table is most pronounced in JULES-Peat-W-drain (water table change given in Table 4).

Figure 8 Carbon and surface height dynamics following drainage: JULES-Peat-W and JULES-Peat-W-drain are both spun up with a wet soil column, but JULES-Peat-W-drain has lateral flows switched back on at the start of the historical simulation, which is shown here. JULES-Peat-B is also shown for comparison, showing that carbon does not accumulate at these sites to the extent that it does when the sites are wet (compare JULES-Peat-B to JULES-Peat-W). Subsidence rates (in cm yr⁻¹) over the first 40 years are indicated in the panel labels. Water table depth in each simulation is given in Table 4.

Table 4 Median water tables corresponding to the simulations shown in Fig. 8.

Liu et al. (2020) tracked the surface subsidence rate over time following drainage of peatlands for two different land use practices – forestry and agriculture. They typically see a very high subsidence rate of around 3–10 cm yr⁻¹ in the first few years after drainage. Following this, a more steady subsidence rate of 0.5–2 cm yr⁻¹ can be observed for the next few decades. In the sites in Fig. 8, we see a surface subsidence rate more in line with the longer-term subsidence rate of 0.5–2 cm yr⁻¹ (e.g. Auchencorth loses > 40 cm in the first 40 years and 60 cm in around 80 years; see Fig. 8). The lack of the initial very rapid subsidence suggests that there may be some processes missing in JULES, for example the mechanical raising and lowering of the peat surface by as much as tens of centimetres as the water table fluctuates, known as bog breathing (Howie and Hebda, 2018). However, the long-term subsidence rate is at least of the right order of magnitude. After an initial period of subsidence lasting around 50 years, the drop in surface height stabilises or slows. The carbon loss behaves in a similar way, although the slowing of carbon loss is not as pronounced. In these test simulations (noting that the method of “draining” the sites is a proof of concept and is not based on reality), up to 17 kg C m⁻² is lost from these sites, which would represent a highly significant addition of carbon to the atmosphere (on the order of tens of gigatons of carbon globally) if it took place across a significant fraction of the world's peatlands.

Both positive and negative feedbacks exist within peatland ecosystems (Waddington et al., 2015). JULES-Peat is able to capture some of the key feedbacks by simulating dynamic soil properties. Firstly, a negative (damping) feedback takes place following drainage in which the peat compacts and becomes more resistant to water flow, thus re-wetting the soil to some extent; see Fig. 9. There is a strong correlation between decreased hydraulic conductivity and increased soil moisture: Pearson correlation between soil moisture and log (K_sat) of −0.94 for Auchencorth and −0.87 for Carlow (in the top 3 m) using monthly data points for individual layers for the whole simulation. In these particular simulations, this effect was not strong enough to bring the water table back to the surface by the end of the simulation, but in some test simulations this effect was observed. In reality, it is rare for drained peatlands to self-restore, but it does occasionally happen (Alice Milner, personal communication, 2021​​​​​​​; Morag Angus, personal communication, 2021). On the other hand, a similar mechanism can lead to a positive (amplifying) feedback during spin-up, where the accumulation of peat leads to a lower drainage rate and thus further accumulation of peat. This is seen most strongly at Auchencorth and Carlow, which are the only sites from the UK and Ireland that were simulated (Fig. 10). After sufficient peat formation (in particular after being initialised with peat in JULES-Peat-i and JULES-Peat-i10; see Table 2), these sites are able to gradually raise their water tables over the course of the spin-up, while accumulating more and more carbon; see Fig. 10. It is significant that the UK and Ireland sites are the only ones where this mechanism takes place, since this is where the majority of the world's blanket bogs are found – peatlands that are able to maintain themselves autogenically without topographic controls. This indicates that JULES-Peat should be capable of simulating blanket bogs, which are currently missing in global peatland models due to their reliance on TOPMODEL-based wetland fraction to determine peatland area (Müller and Joos, 2020).

Figure 9 Response of the soil profile following drainage in JULES-Peat-W-drain.

Figure 10 Hydrological feedback between peat accumulation and water table level. Note that the spin-up time is noted in thousands of years so the total spin-up period is 10 000 years. Vertical profiles of K_sat and carbon are shown at the end of the spin-up.

5.2 Multiple steady states

Since there are feedbacks in JULES-Peat between soil physics and soil carbon that can be self-reinforcing, the “end state” of the model spin-up now depends on the initial conditions. In mathematical terms there can be more than one equilibrium state. This also implies the existence of tipping points where the system can “tip” from one state to another under sufficient forcing (Ritchie et al., 2021). Essentially, there may be some sites at which initialising the model with peat allows it to further accumulate peat, but initialising the model with mineral soil maintains a stable mineral soil profile. Practically, this means that peat can exist outside of climatic conditions where it would form from scratch, and thus this is very important when we consider the disturbance of existing peatlands for which the original state could be impossible to recover under current climates.

Figure 11 Demonstration of how peat accumulation rate can depend on initial conditions. Profiles of carbon and age are shown at the end of spin-up. Note that the ages below a certain depth at Twitchell are not relevant due to the carbon pools being almost zero. Note also that the spin-up time is given in thousands of years (so the total spin-up period is 10 000 years).

In order to test this we compare the simulations where peat is initialised vs. not initialised in Fig. 11. We show two sites where peat only accumulates when it is initialised (Abiskomire_noSnowCor and Iskoras), a site where peat always accumulates (CA_WP1), and a site where peat never accumulates and a mineral soil always forms (Twitchell). This behaviour is apparent in the age–depth profiles (bottom row of Fig. 11). The age is initialised with the existing age–depth profile of the peat for the runs that are initialised with peat (JULES-Peat-i and JULES-Peat-i10), hence the ages overall can be higher, since the standard runs start from age zero with no carbon. However, when the model converges to a single steady state, the age profile also starts to converge, at least at the surface (Twitchell and CA_WP1 in Fig. 11). In contrast, at sites where peat only accumulates when the model is initialised with peat, the age at the surface actually becomes lower in the simulations where it was initialised with an existing age–depth profile than in the simulations where it was initialised at zero. This indicates that carbon is accumulating more quickly when the model is initialised with peat (-i and -i10) (Abiskomire_noSnowCor and Iskoras in Fig. 11).

It is interesting to note that the sites where highly distinct steady states are simulated are palsa mires in the sporadic permafrost zone. These sites are on the cusp of the permafrost–non-permafrost transition. Thawing of permafrost peatlands has been shown to increase the carbon accumulation rate in some cases (Turetsky et al., 2007), and, interestingly, the simulations with high peat accumulation rate at Abiskomire_noSnowCor and Iskoras (JULES-Peat-i and JULES-Peat-i10) simulate a thawed soil, whereas the simulations with little peat accumulation (JULES-Peat-B and JULES-Peat-B10) simulate permafrost (Fig. S11).

Table 5 Dominant vegetation type and fraction at the end of the spin-up for the sites and runs shown in Fig. 11. EG stands for evergreen.

The different steady states also appear to be associated with the presence of different vegetation types; see Table 5. The sites that develop very different carbon profiles when initialised with peat (Abiskomire_noSnowCor and Iskoras, JULES-Peat-i and JULES-Peat-i10) also develop a different dominant vegetation type, for example needle-leaf evergreen trees instead of Arctic grass in Abiskomire_noSnowCor (Table 5). This interaction between vegetation and soil carbon highlights the importance of further developing the vegetation model to represent peatland vegetation (see Sect. 5.3 for further discussion). It is worth noting that for a site that simulates peat accumulation, the “steady-state” condition can be a constant rate of carbon accumulation rather than a constant quantity of carbon, i.e. peatlands are never in equilibrium, which differs from the standard definition of steady state that is currently used when spinning up land surface models.

5.3 Next steps for modelling global peatlands: hydrology and vegetation

Simulating landscape-level peatland hydrology is a major challenge. This work has taken a step forward by enabling peat soils to react appropriately to long-term changes in hydrology and to include some of the key feedbacks that the soils then have on the water flows. We have also recently developed an improved methane emissions scheme (Chadburn et al., 2020). However, to model peatlands within the landscape globally, including any methane emissions, the distribution of water around the landscape must be taken into account.

For instance, the majority of peatlands globally are topographically controlled. This means they are found in flat, lowland areas (Sheng et al., 2004; Martini, 2006). The standard way of modelling groundwater in ESMs does not explicitly model these areas but simulates a “grid cell average” soil moisture. This means that at typical resolutions, the saturated lowland areas where peat forms would be less than the size of a grid cell, and thus saturated conditions would never be explicitly simulated, and secondly, even with a high enough resolution to resolve peatlands, there is no mechanism for lowland grid cell soils to receive water from the surrounding uplands. Therefore, a key step is to explicitly model different hydrological regimes and features within the landscape. The simplest way to do this is via a tiling approach. Bechtold et al. (2019) found that they were able to recreate the hydrological dynamics in the majority of the world's peatlands by using a tile that was entirely hydrologically disconnected from the rest of the grid cell. A further step would be to simulate the hydrological connection between the tiles in the grid cell, as this is necessary to simulate not only some existing peatlands but also peatland initiation (Väliranta et al., 2017).

Furthermore, the within-soil-column hydrology is not well modelled for organic soils in JULES. While peatlands in JULES-Peat are able to maintain a water table through the hydraulic characteristics of the peat itself, the water table generally does not sit as close to the surface as the observed water table in intact peatlands, which is around 10 cm (Evans et al., 2021) (see, e.g. Carlow in Fig. 10), although it can occasionally reach 10 cm when lateral flow is set to zero (Table 4). In cold regions, a representation of ponding can be necessary to prevent too much snowmelt from running off and leaving the soil too dry (Smith et al., 2021). On top of this, the hydraulic behaviour of mosses, which form a primary component of high-latitude and temperate peatlands, is very different from that of vascular plants. Mosses do not extract water from the soil and it essentially only evaporates from the surface, which could very well lead to a raised water table (Van Breemen, 1995). Thus the inclusion of a moss PFT and its unique functions in land surface models like JULES should be a priority, and several models have indeed done this (e.g. Porada et al., 2016; Chaudhary et al., 2017; Shi et al., 2021).

It is clearly important to adequately represent the features of peatland vegetation. As well as the hydrological behaviour of mosses, it will be crucial to include an appropriate distribution of plant litter inputs to the soil (see the difference between the simulations with different values of τ_lit, e.g. in Fig. 5), an appropriate recalcitrant litter fraction (for example, mosses are more recalcitrant than grass and therefore more likely to lead to peat accumulation), and suppression of the growth of non-wetland vegetation such as trees under saturated conditions (this is not included in JULES and is necessary to simulate the mossy peat that is found in northern latitudes, since larger vegetation would otherwise outcompete the mosses). In addition, the input of carbon to the peat is determined by the net primary productivity of the ecosystem, and thus this is a key quantity to evaluate when developing peatland-appropriate PFTs.

Finally, the JULES-Peat model configuration has not yet been tested in tropical peatlands, which differ from northern peatlands in terms of hydrology and vegetation and have only recently gained attention in the modelling community (Kurnianto et al., 2015; Apers et al., 2021). There is a clear need for more focused study of tropical peatlands, given their large spatial extent and carbon stock and the potential impacts of their ongoing drainage (Dargie et al., 2017; Mishra et al., 2021). Some of the key principles behind peat dynamics are universal (for example, suppression of decomposition in wet soils, dynamic growth of the soil surface), but model parameters such as those in the relationships used to determine soil characteristics may need to be updated for tropical peat (e.g. Eqs. 13–18).

6 Conclusions and outlook

We have demonstrated a new scheme integrated in an ESM land surface model that can simulate both peat and mineral soils depending on site conditions and that can simulate dynamic transitions from peat to mineral soil or vice versa. The new model configuration, which we call JULES-Peat, includes some key ecohydrological feedbacks that take place in peat soils. At some sites, whether or not peat accumulates depends on the initial conditions.

The model performs well by all metrics that we compared it against, and it can now simulate a soil profile that resembles peat for the first time in JULES. As well as simulating mechanisms that determine the (in)stability and resilience of peatlands for the first time, this model has the potential to simulate blanket bogs, which current global peatland models are unable to do (Müller and Joos, 2020). We noted when designing the interpolation scheme (Sect. S1) that the interpolation can lead to some “smearing” of the carbon profile in the deeper soil, and indeed the model does not simulate sharp transitions between peat layers and underlying mineral soil that can often be seen in reality. Thus some improvement to the interpolation scheme may still be possible, which could also lead to improved physical soil characteristics. It should also be noted that the JULES soil layers need to be set at a sufficiently high resolution to be capable of resolving such a transition.

As outlined above in Sect. 5.3, major challenges remain around appropriately modelling peatland vegetation and large-scale hydrology, as well as a need to test the model for tropical peatlands. It may also be necessary to model microtopography and/or ponding in order to simulate soil hydrology correctly (Smith et al., 2021). Since models individually tackle different parts of this problem, the next steps will inevitably involve combining existing schemes, or at least concepts, for simulating vegetation, large-scale hydrology and microtopography with the soil dynamics simulated here in JULES-Peat (e.g. Bechtold et al., 2019; Porada et al., 2016; Shi et al., 2021), along with the latest methane emissions schemes (e.g. Chadburn et al., 2020).

Peatlands are of utmost importance in terms of mitigating climate change, both as carbon sinks and as potentially very large carbon sources that may exacerbate climate change (Leifeld and Menichetti, 2018). Modelling global peatlands and their dynamics should therefore be a priority for land surface and Earth system modelling.