Evaluation of soil carbon simulation in CMIP6 Earth system models

Abstract

The response of soil carbon represents one of the key uncertainties in future climate change. The ability of Earth system models (ESMs) to simulate present-day soil carbon is therefore vital for reliably estimating global carbon budgets required for Paris Agreement targets. In this study CMIP6 ESMs are evaluated against empirical datasets to assess the ability of each model to simulate soil carbon and related controls: net primary productivity (NPP) and soil carbon turnover time (${\mathit{\tau }}_{\mathrm{s}}$). Comparing CMIP6 with the previous generation of models (CMIP5), a lack of consistency in modelled soil carbon remains, particularly the underestimation of northern high-latitude soil carbon stocks. There is a robust improvement in the simulation of NPP in CMIP6 compared with CMIP5; however, an unrealistically high correlation with soil carbon stocks remains, suggesting the potential for an overestimation of the long-term terrestrial carbon sink. Additionally, the same improvements are not seen in the simulation of ${\mathit{\tau }}_{\mathrm{s}}$. These results suggest that much of the uncertainty associated with modelled soil carbon stocks can be attributed to the simulation of below-ground processes, and greater emphasis is required on improving the representation of below-ground soil processes in future developments of models. These improvements would help to reduce the uncertainty in projected carbon release from global soils under climate change and to increase confidence in the carbon budgets associated with different levels of global warming.

1 Introduction

Soil carbon is the Earth's largest terrestrial carbon store, with a magnitude of at least 3 times the amount of carbon contained within the atmosphere . The response of soil carbon to CO${}_{\mathrm{2}}$-induced global warming has the potential to provide a significant feedback to climate change, but this feedback is currently poorly known . Carbon stored within the atmosphere and global soils is exchanged via carbon fluxes as part of the global carbon cycle . The Earth's terrestrial surface has acted as a carbon sink until now , but there is a possibility of a switch to a source during the 21st century, which would accelerate climate change . Due to the significant quantities of carbon stored in soils globally, understanding and quantifying the potential release of carbon from soils is vital if the existing Paris Agreement targets are to be met .

Earth system models (ESMs) are complex numerical models which simulate both climate and carbon cycle processes and are used to make projections of climate change. The latest generation of the Coupled Model Intercomparison Project (CMIP), CMIP6 , includes an ensemble of ESMs, which are used in the most recent Intergovernmental Panel on Climate Change (IPCC) report (AR6) . The relationships between carbon and environmental drivers used in models help to determine the response of the carbon cycle to climate change . Therefore, representing present-day carbon stores and spatial controls realistically within models is key for estimating carbon emission cuts required for Paris Agreement targets .

Present-day soil carbon can be approximately broken down into above-ground and below-ground controls, which influence the spatial distribution of soil carbon stocks . The above-ground control of soil carbon can be considered as the input flux of carbon into the soil from vegetation. Both the amount of carbon from plant and root litter (known as litterfall) and the fraction of this that is converted to longer-lived soil carbon pools will influence the storage of soil carbon. Net primary productivity (NPP) can be used as a proxy for the litterfall flux, where the fluxes are equal when vegetation is in a steady state. The below-ground control of soil carbon can be quantified simply in terms of the soil carbon turnover time (${\mathit{\tau }}_{\mathrm{s}}$), which is defined as the time carbon resides in the soil . ${\mathit{\tau }}_{\mathrm{s}}$ can be considered as a proxy for below-ground controls on soil carbon storage .

In this study, the representation of late 20th century soil carbon stores and these related controls (NPP and ${\mathit{\tau }}_{\mathrm{s}}$) is evaluated in CMIP6 ESMs. Previously, similar studies have been conducted to evaluate soil carbon in the preceding generations of ESMs, for example, and for CMIP5. There are some existing CMIP6 soil-carbon-related studies; for example, evaluated carbon-concentration and carbon-climate feedbacks in 1 % CO${}_{\mathrm{2}}$ yr${}^{-\mathrm{1}}$ forcing simulations, evaluated the representation of permafrost in models, and investigated future soil carbon stocks under specific land use conditions. This study is the first to specifically focus on global and spatial soil carbon and related controls in CMIP6, with a thorough evaluation against empirical datasets and comparison against the preceding CMIP5 ensemble.

2 Methods

2.1 Earth system models

Soil carbon stores and related controls are examined in 11 CMIP6 ESMs , as listed in Table . Throughout the study, comparisons are made with 10 ESMs from the previous CMIP generation (CMIP5, ), as listed in Table . The ESMs included in this study were chosen due to the availability of the required data in the online repository at the time of analysis (https://esgf-node.llnl.gov/search/cmip6/ (last access: 8 April 2022) https://esgf-node.llnl.gov/search/cmip5/, (last access: 12 April 2022)).

Table 1

The 11 CMIP6 Earth system models included in this study and relevant features of their land carbon cycle components .

Earth system	Modelling	Land surface	Nitrogen	No. of live	No. of dead	Temperature	References
model	centre	model	cycle	carbon pools	carbon pools	and moisture
ACCESS-ESM1.5	CSIRO	CABLE2.4	Yes	3	6	Arrhenius
		$+$ CASA-CNP				and hill
BCC-CSM2-MR	BCC	BCC-AVIM2	No	3	8	Hill
						and hill
CanESM5	CCCma	CLASS-CTEM	No	3	2	$Q_{\mathrm{10}}$
						and hill
CESM2	CESM	CLM5	Yes	22	7	Arrhenius
						and increasing
CNRM-ESM2-1	CNRM	ISBA-CTRIP	No	6	7	$Q_{\mathrm{10}}$
						and increasing
GFDL-ESM4	GFDL	LM4.1	No	6	4	Hill
						and increasing
IPSL-CM6A-LR	IPSL	ORCHIDEE	No	8	3	$Q_{\mathrm{10}}$
		branch 2.0				and increasing
MIROC-ES2L	JAMSTEC	MATSIRO	Yes	3	6	Arrhenius
		VISIT-s				and increasing
MPI-ESM1.2-LR	MPI	JSBACH3.2	Yes	3	18	$Q_{\mathrm{10}}$
						and increasing
NorESM2-LM	NCC	CLM5	Yes	22	7	Arrhenius
						and increasing
UKESM1-0-LL	UK	JULES-ES-1.0	Yes	3	4	$Q_{\mathrm{10}}$
						and hill

Table 2

The 10 CMIP5 Earth system models included in this study and relevant features of their land carbon cycle components , including temperature and moisture functions presented in .

Earth system	Modelling	Land surface	Nitrogen	No. of live and dead	Temperature	References
model	centre	model	cycle	carbon pools	and moisture
BNU-ESM	BNU	CoLM $+$	Yes	–	$Q_{\mathrm{10}}$
		BNU-DGVM			and increasing
CCSM4	CCSM	CLM4	Yes	20	Arrhenius
					and increasing
CanESM2	CCCma	CLASS2.7	No	5	$Q_{\mathrm{10}}$
		$+$ CTEM1			and hill
GFDL-ESM2G	GFDL	LM3	No	10	Hill
					and increasing
GISS-E2-R	NASA-GISS	YIBs,	No	12	Increasing
		version 1.0			and increasing
HadGEM2-ES	MOHC	JULES $+$	No	7	$Q_{\mathrm{10}}$
		TRIFFID			and hill
IPSL-CM5A-LR	IPSL	ORCHIDEE	No	7	$Q_{\mathrm{10}}$
					and increasing
MIROC-ESM	JAMSTEC	MATSIRO	No	6	Arrhenius
		$+$ SEIB-DGVM			and increasing
MPI-ESM-LR	MPI	JSBACH	No	6	$Q_{\mathrm{10}}$
		$+$ BETHY			and increasing
NorESM1-M	NCC	CLM4	Yes	20	Arrhenius
					and increasing

Tables  (CMIP6) and (CMIP5) present information about the included ESMs, specifically more details about the associated land surface model (LSM). It should be noted that there are similarities between some of the LSMs – either advances from earlier models or even the same LSM within different ESMs. For example, CESM2 and NorESM2-LM both use the Community Land Model version 5 (CLM5) . For some modelling centres, both the CMIP5 and CMIP6 versions of the models are included, and in these cases direct comparisons can be made to determine changes from CMIP5 to CMIP6. These generationally related CMIP5 and CMIP6 models are CanESM2 and CanESM5, CCSM4 and CESM2, GFDL-ESM2G and GDFL-ESM4, IPSL-CM5A-LR and IPSL-CM6A-LR, MIROC-ESM and MIROC-ES2L, MPI-ESM-LR and MPI-ESM1.2-LR, NorESM1-M and NorESM2-LM, and HadGEM2-ES and UKESM1-0-LL, respectively. The models where only either the CMIP5 or CMIP6 version from the modelling centre was included are BNU-ESM and GISS-E2-R from CMIP5 and ACCESS-ESM1.5, BCC-CSM2-MR, and CNRM-ESM2-1 from CMIP6. A key general change to note is that CMIP6 has more models that include an interactive nitrogen cycle compared with CMIP5: ACCESS-ESM1.5, CESM2, MIROC-ES2L, MPI-ESM1.2-LR, NorESM2-LM, and UKESM1-0-LL in CMIP6 compared with CCSM4 and NorESM1-M in CMIP5 (the CMIP5 model BNU-ESM includes carbon–nitrogen interactions; however, this process was turned off in CMIP5 simulations; ). Additionally, an increased number of soil carbon pools is seen in some CMIP6 models (e.g. CLM5 has 29 carbon pools compared with 20 in CLM4). include a comprehensive overview of the updates seen in the individual CMIP6 models, which is presented in the “Model descriptions” section of the associated Appendix.

include a summary of the temperature and moisture dependencies of soil respiration/decomposition as assumed in the CMIP5 models (see Table 1 in ). The most common representation of the temperature sensitivity of decomposition is the $Q_{\mathrm{10}}$ equation, which is defined by $f\left(T\right)=Q_{\mathrm{10}}^{(T-T_{\mathrm{0}})/\mathrm{10}}$, where $T$ is temperature and $T_{\mathrm{0}}$ is a reference temperature. With the $Q_{\mathrm{10}}$ equation, decomposition increases exponentially with temperature . The majority of the other models used the Arrhenius equation to represent the temperature sensitivity, where the main difference from the $Q_{\mathrm{10}}$ representation is that decomposition levels off at higher temperature levels . Of the remaining models, the GFDL model simulates an increased decomposition with temperature until some optimal temperature above which it decreases (which , defined as a “hill” function) and the GISS model implements a linear increase in respiration to temperature up to a maximum value . The representation of the decomposition sensitivity to soil moisture was found to be represented in two ways amongst the CMIP5 models, where decomposition was assumed either to increase monotonically with increasing soil moisture or less commonly to increase to some optimum moisture level and then decrease (again described as a “hill” function by ). In this study we note that the representation of temperature and moisture functions remains similar from CMIP5 to CMIP6. The $Q_{\mathrm{10}}$ equation remains the most common representation of soil temperature sensitivity in models, followed by the Arrhenius equation and then “hill” functions. Similarly, the most common representation of the sensitivity of soil to moisture in CMIP6 is a monotonically increasing function, followed by “hill” functions of various sorts.

2.2 Defining soil carbon variables

CMIP defines common output variables , which allows for consistent comparison between the models and for cleaner evaluation of models against observational data. These common output variables also allow for consistent comparison between model generations, in this case between CMIP6 and CMIP5. This study focuses on evaluation of near-present-day soil carbon and related controls. Therefore the results presented in this study use the CMIP standard historical simulation (CMIP scenario historical) for both the CMIP6 and CMIP5 analyses. The historical simulation runs from 1850 to 2015 in CMIP6 and from 1850 to 2005 in CMIP5, where the selected dates for each variable (stated below) were chosen to allow for consistent comparison between CMIP5 and CMIP6 and to best match the modelled data to the empirical data.

To evaluate soil carbon, this study uses “Soil Carbon” (CMIP variable cSoil), which represents the carbon stored in soils, and where applicable “Litter Carbon” (CMIP variable cLitter), which represents carbon stored in the vegetation litter. Total soil carbon ($C_{\mathrm{s}}$) is defined as the sum of these soil carbon and litter carbon variables (cSoil $+$ cLitter), where for models that do not report a separate litter carbon pool, the total soil carbon is taken to be simply the cSoil variable. This allows for a more consistent comparison between the models and between the models and empirical data due to differences in how soil carbon and litter carbon are simulated . Modelled $C_{\mathrm{s}}$ is time-averaged between the years 1950 and 2000 of the historical simulation and is considered spatially (kg m${}^{-\mathrm{2}}$) and as global totals (PgC), where global totals are calculated as an area-weighted sum using the model land surface fraction (CMIP variable sftlf). To calculate northern latitude totals, a sum between the latitudes 60 and 90${}^{\circ }$ N is considered.

The CMIP6 ESMs CESM2 and NorESM2-LM have two different variables to represent soil carbon: (1) CMIP variable cSoil, which represents the full vertical soil profile, and (2) CMIP variable cSoilAbove1m, which represents soil carbon in the top 1 m of soil. This is due to the representation of vertically resolved soil carbon in these models, which means that there are separate carbon pools in the model that represent different soil depths . The CMIP variable cSoilAbove1m is used throughout this study to represent soil carbon for the models CESM2 and NorESM2-LM, unless otherwise stated. The use of this variable is to enable a more consistent comparison with both the other CMIP6 models and the CMIP5 models. Therefore, an assumption of a 1 m depth of soil for modelled soil carbon allows for the fairest evaluation, and evaluation is considered against empirical datasets down to a depth of 1 m (see below). However, comparisons with the cSoil variable for both CESM2 and NorESM2-LM are included in Tables  and of the results.

In order to obtain a clean separation between above-ground and below-ground drivers of soil carbon variations, a quasi-equilibrium approximation is made. We begin with the definition of the effective soil carbon turnover time (${\mathit{\tau }}_{\mathrm{s}}$) , which represents the average time carbon resides in the soil:

1$${\mathit{\tau }}_{\mathrm{s}}={\displaystyle \frac{C_{\mathrm{s}}}{R_{\mathrm{h}}}},$$ where $R_{\mathrm{h}}$ is the output flux of carbon from the soil known as the heterotrophic respiration, which is described as the carbon loss from the decomposition by microbes. This definition of the turnover time implicitly neglects other processes which result in soil carbon release but which are not yet routinely included in ESMs (e.g. peat fires or dissolved organic carbon fluxes).

The definition of the effective turnover time (Eq. ) ensures that the soil carbon at any one time is given by $C_{\mathrm{s}}=R_{\mathrm{h}}{\mathit{\tau }}_{\mathrm{s}}$. In an unperturbed steady state (i.e. neglecting disturbances from land use change, fires, insect outbreaks, etc.), there is no net exchange of carbon between land and atmosphere, and therefore $R_{\mathrm{h}}$ is equal to litterfall, known as fallen organic material from plants. When vegetation and soil carbon are close to a steady state, litterfall and $R_{\mathrm{h}}$ are also approximately equal to NPP, where NPP is defined as the net carbon assimilated by plants via photosynthesis minus loss due to plant respiration. In the contemporary period considered in this study, $R_{\mathrm{h}}$ has been found to be well approximated by NPP . This is because the difference between NPP and $R_{\mathrm{h}}$, which represents the net ecosystem productivity (NEP), is a small fraction of the NPP over the historical period (NPP $\approx$ 60 PgC yr${}^{-\mathrm{1}}$; NEP $\approx$ 3 PgC yr${}^{-\mathrm{1}}$). Therefore the present-day soil carbon can be approximated by 2$$C_{\mathrm{s}}\approx \mathrm{NPP}{\mathit{\tau }}_{\mathrm{s}}$$ to a good accuracy. This allows for a clean separation of soil carbon variation into the above (NPP) and below (${\mathit{\tau }}_{\mathrm{s}}$) ground drivers of soil carbon spatial patterns, following the approach of previous published studies .

To evaluate these soil carbon controls on $C_{\mathrm{s}}$, NPP and ${\mathit{\tau }}_{\mathrm{s}}$ are evaluated separately. This study uses modelled “net primary productivity” (CMIP variable npp), which is defined as the mass flux of carbon out of the atmosphere due to NPP on land. NPP is also considered spatially (kg m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$) and as an area-weighted global total flux (PgC yr${}^{-\mathrm{1}}$). By definition, ${\mathit{\tau }}_{\mathrm{s}}$ is defined by Eq. () and therefore is calculated by soil carbon (as defined above) divided by $R_{\mathrm{h}}$. For $R_{\mathrm{h}}$, the variable “heterotrophic respiration” (CMIP variable rh) is used, which is defined as the mass flux of carbon into the atmosphere due to heterotrophic respiration on land, primarily due to the microbial respiration that occurs in the soil and where the units of $R_{\mathrm{h}}$ are the same as those of NPP. The carbon fluxes (NPP and $R_{\mathrm{h}}$) are time-averaged over the period 1995 to 2005 for consistency between the CMIP generations and the empirical datasets. ${\mathit{\tau }}_{\mathrm{s}}$ can be considered on a spatial level or as an effective global ${\mathit{\tau }}_{\mathrm{s}}$, which is defined as the average ${\mathit{\tau }}_{\mathrm{s}}=$ mean($C_{\mathrm{s}}$)/mean($R_{\mathrm{h}}$) (where the mean represents an area-weighted global average). The advantage of defining an effective global ${\mathit{\tau }}_{\mathrm{s}}$ is that it is not dominated by large spatial outlying values. Using either method, the units for ${\mathit{\tau }}_{\mathrm{s}}$ are in years (yr) by definition.

The relationships of $C_{\mathrm{s}}$, NPP, and ${\mathit{\tau }}_{\mathrm{s}}$ with both temperature and soil moisture are also considered. For temperature, the variable “near-surface air temperature” (CMIP variable tas) representing atmospheric temperature at the surface is considered, where the dates 1995 to 2005 were chosen to be consistent with the carbon fluxes. The variable for atmospheric temperature is considered as opposed to soil temperature as equivalent global observational datasets are required for the analysis. For soil moisture, the variable “moisture in the upper portion of the soil column” (CMIP variable mrsos), which is defined as the mass content of water in the soil layer in the upper portion of the soil (0–10 cm depth), is considered, where the dates 1978 to 2000 were considered to match the empirical data. The standard output mrsos is in units of kg m${}^{-\mathrm{2}}$; however, in this study a volumetric soil moisture, referred to as $\mathit{\theta }$, is used to allow for consistent comparison with the benchmark data. $\mathit{\theta }$ is calculated as mrsos divided by the depth of the soil layer in millimetres, which in this case is $\mathit{\theta }=$ mrsos $/\mathrm{100}$. The variable mrsos for soil moisture was considered opposed to the full soil column moisture (CMIP variable mrso) as this better matched the available empirical dataset for soil moisture. It is noted that this represents surface soil moisture and does not match the depth over which soil carbon is evaluated (0–1 m). This is due to deeper soil moisture products not being as readily available due to limitations of remote sensing methods in penetrating deeper ground. It is expected that the surface soil moisture will be related to deeper soil moisture to some extent but will be influenced by different processes. For example, high surface soil moisture after rainfall events could run off and thus not always reach the deeper soil.

2.3 Empirical datasets

2.3.1 Soil carbon

Observational $C_{\mathrm{s}}$ to a depth of 1 m was obtained by combining the empirical Harmonized World Soils Database (HWSD) and Northern Circumpolar Soil Carbon Database (NCSCD) soil carbon datasets, where NCSCD was used where overlap of the datasets occurs. This is a commonly used method when considering empirical soil carbon and has been previously used in multiple studies, such as , , and . This dataset is referred to here as the “benchmark dataset”.

We use the 95 % confidence intervals given by to derive standard deviations about the global mean soil carbon. To do this, the constructed 95 % confidence intervals were used to calculate upper and lower bounds around the mean value. Then, assuming the data are normally distributed, these derived 95 % confidence intervals were halved to obtain confidence intervals equivalent to a standard deviation error on the mean (1412 $\pm$ 215 PgC). The uncertainty analysis completed in is used for the benchmark soil carbon dataset as no quantitative uncertainty has been previously or since defined for the HWSD and NCSCD datasets .

Additionally, the benchmark dataset was compared with empirical estimates found in the literature to improve the robustness and reliability of the evaluation. find that this derived uncertainty is consistent with other empirical estimates of global soil carbon, for example 1576 PgC in , 1220 PgC in , and 1502 PgC in . This study further compares with empirical estimates of 1395 PgC in and 1515 PgC in . These empirical estimates are within 1 standard deviation of the global mean soil carbon given by the benchmark dataset (Table ).

Table 3

Table of global total and northern latitude total (northern latitudes defined as 60–90${}^{\circ }$ N) soil carbon estimates from multiple empirical datasets, for varying soil depths where applicable.

Empirical	Depth	Global	Northern	Reference
dataset		total	latitude
		$C_{\mathrm{s}}$ (PgC)	total
			$C_{\mathrm{s}}$ (PgC)
HWSD $+$ NCSCD	1 m	1412 $\pm$ 215	401 $\pm$ 61
WISE30sec	1 m	1371 $\pm$ 129	314
	2 m	1952 $\pm$ 198	468
S2017	1 m	1966	515
	2 m	3141	893
GSDE	1 m	1682	526
	2.3 m	2593	849
IGBP DIS	1 m	1567	377

Moreover, additional empirical datasets are considered to improve the reliability of the benchmark dataset (Table ). These additional datasets include the following: (1) the World Inventory of Soil property Estimates (WISE30sec) dataset down to a depth of 2 m , which includes a given standard deviation on the global total soil carbon consistent with our derived benchmark uncertainty, (2) the named “S2017” from the soil carbon estimate (1 and 2 m), which uses a data-driven statistical model and the History Database of the Global Environment (HYDE) land use data, (3) the Global Soil Dataset for use in Earth system models (GSDE), which provides estimates for observational soil carbon down to a depth of up to 2.3 m , and (4) the Global Gridded Surfaces of Selected Soil Characteristics (IGBP-DIS) estimate of soil carbon to a depth of 1 m, derived by the Oak Ridge National Laboratory Distributed Active Archive Centre (ORNL DAAC) . These datasets were combined to obtain a mean estimate for observational soil carbon down to a depth of 1 m, where a global total soil carbon value of 1560 $\pm$ 214 PgC was found. This estimate is consistent with our benchmark dataset estimate and further improves the confidence in our benchmark dataset.

Furthermore, the spatial correlation coefficients between these additional datasets and our benchmark dataset are considered, where the following values correspond to the above datasets: (1) 0.554, (2) 0.625, (3) 0.482, and (4) 0.622. Map plots comparing the empirical soil carbon datasets are shown in Fig. . The estimate for northern latitude total soil carbon has greater uncertainties associated with it, where the standard deviation deduced by combining the empirical datasets is $\pm$ 83 PgC. To account for this increased uncertainty in these regions, the deduced standard deviation of $\pm$ 83 PgC is used on the benchmark soil carbon throughout this study, as opposed to the $\pm$ 61 PgC derived using the uncertainty analysis.

2.3.2 Carbon fluxes

To estimate a benchmark NPP, the commonly used MODIS NPP (2000–2010) dataset is used. The MODIS NPP dataset does not have associated uncertainty estimates, so this study estimates a standard deviation error on the benchmark NPP as derived by . The MODIS NPP dataset is found to be consistent with 251 empirical present-day estimates of NPP found in the literature, which used to estimate a global value of 56.2 $\pm$ 14.3 PgC yr${}^{-\mathrm{1}}$ (compared with a derived MODIS mean value of 56.6 yr${}^{-\mathrm{1}}$).

Moreover, due to the limited choice of observationally derived NPP datasets , models can be further evaluated using a benchmark dataset for $R_{\mathrm{h}}$, where $R_{\mathrm{h}}$ is estimated using the CARDAMOM (2001–2010) heterotrophic respiration dataset . The empirical CARDAMOM $R_{\mathrm{h}}$ has associated estimates of error, which were used to derive a standard deviation uncertainty on the empirical average $R_{\mathrm{h}}$ (51.7 $\pm$ 21.8 PgC yr${}^{-\mathrm{1}}$). This study includes map plots comparing the two empirical datasets, which is shown in Fig. . Global totals for $R_{\mathrm{h}}$ are also considered for comparison against NPP, where the CMIP6 and CMIP5 values are also shown in Appendix Tables  and , respectively.

2.3.3 Soil carbon turnover time

To estimate a benchmark ${\mathit{\tau }}_{\mathrm{s}}$, the estimates of observational $C_{\mathrm{s}}$ are divided by an estimate of $R_{\mathrm{h}}$ (see Eq. ). To estimate an uncertainty on the effective global ${\mathit{\tau }}_{\mathrm{s}}$, this study derived upper (${\mathit{\tau }}_{\mathrm{s}}^{+}$) and lower (${\mathit{\tau }}_{\mathrm{s}}^{-}$) bounds based on the derived $C_{\mathrm{s}}$ and $R_{\mathrm{h}}$ uncertainty estimates. The upper bound was calculated using the following: ${\mathit{\tau }}_{\mathrm{s}}^{+}=$ $C_{\mathrm{s}}^{+}$ / $R_{\mathrm{h}}^{-}$, where $C_{\mathrm{s}}^{+}$ is equal to the mean soil carbon plus 1 standard deviation and $R_{\mathrm{h}}^{-}$ is equal to the mean heterotrophic respiration minus 1 standard deviation. The lower bound was calculated using the following: ${\mathit{\tau }}_{\mathrm{s}}^{-}=C_{\mathrm{s}}^{-}$ / $R_{\mathrm{h}}^{+}$, where similarly $C_{\mathrm{s}}^{-}$ is equal to the mean soil carbon minus 1 standard deviation and $R_{\mathrm{h}}^{+}$ is equal to the mean heterotrophic respiration plus 1 standard deviation.

This method gives a large uncertainty bound around the derived mean estimate (27.0${}_{-\mathrm{11}}^{+\mathrm{27}}$ years), so the benchmark data are further compared with empirical estimates. derived an estimate of mean soil carbon turnover of 32 years, using estimates for mean soil carbon pools and mean soil respiration rates. More recently, derived an estimate for the mean global ecosystem carbon turnover time of 23${}_{-\mathrm{4}}^{+\mathrm{7}}$ years, which is a spatially explicit and observation-based estimate. derived an observational uncertainty range on soil carbon turnover time of 18.5 to 45.8 years, which was derived using similar empirical estimates found in the literature. These estimates give more certainty on the values closer to the derived empirical mean value for ${\mathit{\tau }}_{\mathrm{s}}$.

2.3.4 Soil moisture and air temperature

To estimate soil moisture ($\mathit{\theta }$), the Copernicus Climate Change Service (C3S) “Soil moisture gridded data from 1978 to present” dataset (published 25 October 2018) is used, where the years 1978 to 2000 are considered. This dataset is based on the ESA Climate Change Initiative soil moisture and estimates global surface soil moisture from a large set of satellite sensors . The WFDEI meteorological forcing dataset is used to represent observational air temperatures (1995–2005) , where dates are chosen to allow for consistency between CMIP generations. This study includes no uncertainty analysis of the soil moisture and air temperature empirical datasets as these datasets are only used to evaluate spatial correlations between variables and not to evaluate soil moisture and air temperature in the models.

2.4 Regridding

To allow direct comparisons between the empirical data and model output data, the model data were regridded to match the observational grid. In this case, the observational grid is a 0.5${}^{\circ }$ by 0.5${}^{\circ }$ resolution, 720 longitude and 360 latitude grid. The regridding was done using Iris – the community-driven Python package for analysing and visualising Earth science data . The regridding method assumed conservation of mass and used linear extrapolation, where extrapolation points will be calculated by extending the gradient of the closest two points. Moreover, model land masks are used to calculate the fraction of land in each coastal grid cell (CMIP variable sftlf).

2.5 Statistical analysis

It is difficult to evaluate the spatial distributions of modelled soil carbon and related spatial controls against empirical data with a single metric, so the evaluation for both CMIP6 and CMIP5 involves multiple methods. These include coefficients of variation, spatial standard deviations, spatial Pearson correlation coefficients, and root mean square errors (RMSEs). These methods can be combined to give a more thorough evaluation of spatial soil carbon and associated controls in the CMIP6 models compared with the previous generation of CMIP5 models.

The coefficient of variation is defined as the ratio of the ensemble standard deviation (SD) to the ensemble mean in each grid cell. This is used to show the amount of variability amongst the models in the ensemble scaled to the size of the ensemble mean and so represents the variability spatially in the ensemble and shows how much variation is present across the ensemble in specific regions. It is presented as hatching in a map figure (Fig. ), where shaded “hatched” regions show regions of high variability within the ensemble. These regions show areas where there is disagreement in the ensemble as there is a large spread compared with the mean and was defined as being where SD $/$ mean $>\mathrm{0.75}$. The regions where spatial $C_{\mathrm{s}}<\mathrm{5}$ kg m${}^{-\mathrm{2}}$ were discounted as low values of soil carbon are present in these regions.

The spatial standard deviation is a measure of the spread in the data across the globe compared with the mean value. Pearson correlation coefficients ($r$ values) were used as a spatial measure of the linear correlation between the empirical and modelled data, where a high $r$ value (near 1 or $-\mathrm{1}$) represents a high correlation in the data and a low $r$ value (near 0) represents a negligible correlation. RMSEs were used as an absolute measure of the difference between the modelled and empirical data, where the lower the value, the lower the difference error. The RMSE can be considered as the standard deviation of the difference, and it is a measure to show the deviation of the modelled data in relation to the empirical data. These statistical data, spatial standard deviations, Pearson correlation coefficients, and RMSEs can be presented using a Taylor diagram. A Taylor diagram is a graph used to indicate the performance of a model compared with a benchmark, which in this case are the empirical datasets .

3 Results

3.1 Soil carbon stocks: northern latitude underestimations remain in CMIP6

3.1.1 Global total evaluation

Global total soil carbon (in the top 1 m of soil) is shown to vary amongst the ESMs in CMIP6, with a range of 1294 PgC between the models with the lowest and highest values (Table ). The global total soil carbon for 2 (CanESM5 and MIROC-ES2L) out of the 11 CMIP6 models falls within the benchmark soil carbon uncertainty range, 1197–1627 PgC (mean $\pm$ stand deviation). The models with the largest global total soil carbon are CNRM-ESM2-1 (1810 PgC), BCC-CSM2-MR (1770 PgC), and UKESM1-0-LL (1760 PgC), values greater than the benchmark dataset but not the additional empirical datasets (Table ). The models GFDL-ESM4 (516 PgC) and IPSL-CM6A-LR (639 PgC) have the lowest global total soil carbon values in the ensemble, with global totals significantly lower (approximately 50 % less) than the global totals seen in empirical data. It is noted that, if the full soil carbon profile is considered for CESM2 and NorESM2-LM opposed to a depth of 1 m, the global total soil carbon values are increased to 1870 PgC from 991 PgC in CESM2 and to 2430 PgC from 969 PgC in NorESM2-LM.

Table 4

Table presenting global soil carbon values for the 11 CMIP6 models included in this study and the benchmark datasets, including global total $C_{\mathrm{s}}$ in PgC, northern latitude total (90–60${}^{\circ }$ N) $C_{\mathrm{s}}$ in PgC, and the spatial mean value of $C_{\mathrm{s}}$ with the corresponding standard deviation in kg m${}^{-\mathrm{2}}$.

Earth system model	Global total $C_{\mathrm{s}}$ (PgC)	Northern latitude total $C_{\mathrm{s}}$ (PgC)	Mean $C_{\mathrm{s}}\pm$ SD (kg m${}^{-\mathrm{2}}$)
ACCESS-ESM1.5	900	151	5.86 $\pm$ 5.35
BCC-CSM2-MR	1770	575	11.6 $\pm$ 16.6
CanESM5	1500	218	3.87 $\pm$ 6.52
CESM2 (cSoilAbove1m)	991	294	7.05 $\pm$ 16.6
CESM2 (cSoil)	1870	1036	13.8 $\pm$ 51.7
CNRM-ESM2-1	1810	440	12.2 $\pm$ 9.98
GFDL-ESM4	516	163	1.36 $\pm$ 3.43
IPSL-CM6A-LR	639	66.0	4.80 $\pm$ 3.37
MIROC-ES2L	1460	347	9.31 $\pm$ 10.7
MPI-ESM1.2-LR	970	175	6.68 $\pm$ 5.23
NorESM2-LM (cSoilAbove1m)	969	300	2.61 $\pm$ 6.97
NorESM2-LM (cSoil)	2430	1563	6.60 $\pm$ 41.3
UKESM1-0-LL	1760	194	12.0 $\pm$ 10.9
Ensemble mean	1206 $\pm$ 445	266 $\pm$ 139	2.80 $\pm$ 5.15
Benchmark dataset	1412 $\pm$ 215	401 $\pm$ 83	10.7 $\pm$ 9.28

Both the CMIP5 and CMIP6 ensemble mean global totals fall within the benchmark uncertainty range (Tables  and ). The ensemble mean global total soil carbon is found to have decreased in CMIP6 from CMIP5 (1206 $\pm$ 445 PgC vs. 1480 $\pm$ 810 PgC). However, a significant reduction is seen in the associated standard deviation of the ensemble mean global totals in CMIP6 compared with CMIP5 ($\pm$ 445 PgC in CMIP6 from $\pm$ 810 PgC in CMIP5) and a reduced range of global total values (a range of 1294 PgC is seen in CMIP6 as opposed to 2493 PgC in CMIP5). This suggests that, although a significant range in global soil carbon still exists amongst the CMIP6 ESMs, there is an improved consistency between the models seen in CMIP6 compared with the models in CMIP5, although it is noted that this may be a factor of the selection of models included in each ensemble rather than any change in process representation.

Table 5

Table presenting global soil carbon values for the 10 CMIP5 models included in this study and the benchmark datasets, including global total $C_{\mathrm{s}}$ in PgC, northern latitude total (90–60${}^{\circ }$ N) $C_{\mathrm{s}}$ in PgC, and the spatial mean value of $C_{\mathrm{s}}$ with the corresponding standard deviation in kg m${}^{-\mathrm{2}}$.

Earth system model	Global total $C_{\mathrm{s}}$ (PgC)	Northern latitude total $C_{\mathrm{s}}$ (PgC)	Mean $C_{\mathrm{s}}\pm$ SD (kg m${}^{-\mathrm{2}}$)
BNU-ESM	681	135	5.31 $\pm$ 4.55
CCSM4	507	28.1	4.03 $\pm$ 3.24
CanESM2	1540	300	9.16 $\pm$ 9.11
GFDL-ESM2G	1420	635	9.47 $\pm$ 13.2
GISS-E2-R	2150	609	15.9 $\pm$ 20.8
HadGEM2-ES	1080	148	8.19 $\pm$ 6.24
IPSL-CM5A-LR	1350	346	9.77 $\pm$ 7.64
MIROC-ESM	2550	742	20.5 $\pm$ 15.1
MPI-ESM-LR	3000	204	23.5 $\pm$ 14.8
NorESM1-M	538	31.0	3.61 $\pm$ 3.34
Ensemble mean	1480 $\pm$ 810	318 $\pm$ 246	10.5 $\pm$ 6.02
Benchmark dataset	1412 $\pm$ 215	401 $\pm$ 83	10.7 $\pm$ 9.28

It is found from comparing the previous-generation models in CMIP5 with the updated CMIP6 equivalent that multiple models in CMIP6 have lower quantities of soil carbon than in CMIP5, such as GFDL-ESM4 from GFDL-ESM2G, IPSL-CM6A-LR from IPSL-CM5A-LR, MIROC-ES2L from MIRCO-ESM, and MPI-ESM1.2-LR from MPI-ESM-LR. For example, the CMIP5 model MPI-ESM-LR is reported to have the largest soil carbon magnitude amongst the CMIP5 models, with a global total of 3000 PgC (Table ), whereas the updated CMIP6 model MPI-ESM1.2-LR has a reduced global total soil carbon value of 970 PgC amongst the lowest values reported in CMIP6 and below the observationally derived range (Table ). Conversely, these reductions are negated in the ensemble mean by the remaining models, which have greater quantities of soil carbon in CMIP6 compared with their CMIP5 equivalent, such as CanESM5 from CanESM2, CESM2 from CCSM4, NorESM2-LM from NorESM1-M, and UKESM1-0-LL from HadGEM2-ES. For example, the CMIP5 model NorESM1-M is amongst the lowest soil carbon values presented in this ensemble at 538 PgC (Table ), whereas the updated CMIP6 model NorESM2-LM has an increased global total of 969 PgC (down to 1 m) (Table ).

3.1.2 Northern latitude total evaluation

Northern latitude soil carbon (down to a depth of 1 m and where northern latitudes are defined as 60–90${}^{\circ }$ N) is found to be underestimated in CMIP6, with 8 out of the 11 CMIP6 models having lower northern latitude soil carbon values than the derived observational range (Table ). A total of 2 of the 11 CMIP6 models (CNRM-ESM2-1 and MIROC-ES2L) have northern latitude totals that fall within the uncertainty range derived from the benchmark data, 318–484 PgC (mean $\pm$ stand deviation). The CMIP6 models with the greatest northern latitude total soil carbon are BCC-CSM2-MR (575 PgC), CNRM-ESM2-1 (440 PgC), and MIROC-ES2L (347 PgC). The CMIP6 models with the lowest northern latitude soil carbon are IPSL-CM6A-LR (66 PgC), ACCESS-ESM1.5 (151 PgC), GFDL-ESM4 (163 PgC), MPI-ESM1.2-LR (175 PgC), and UKESM1-0-LL (194 PgC), values significantly lower than the totals seen in empirical data.

The northern latitude soil carbon total was also underestimated in CMIP5, with 6 out of the 10 CMIP5 models estimating northern latitude totals lower than the empirical estimates (Table ). The ensemble mean total northern latitude soil carbon is lower in CMIP6 (266 $\pm$ 139 PgC seen in Table ) than in CMIP5 (318 $\pm$ 246 PgC seen in Table ), which is consistent with the global total results; however, both the CMIP5 and CMIP6 mean values fall below the benchmark range. Similarly, as with global soil carbon, a smaller standard deviation on the mean is found for CMIP6 compared with CMIP5, and there is a reduced range in simulated northern latitude total values amongst the CMIP6 models, where despite a large range seen (66 to 575 PgC), an even greater range is seen in CMIP5 (28.1 to 742 PgC). Moreover, improvements are seen amongst models from CMIP5 to CMIP6. For example, the CMIP5 model NorESM1-M had a northern latitude total soil carbon value of 31.0 PgC, which is significantly lower than what is expected based on the benchmark dataset (Table ). However, the updated CMIP6 version of this model, NorESM2-LM, has a northern latitude total soil carbon value of 300 PgC, which is much more in line with the expected observational values (Table ). An improved representation of northern latitude soil carbon is also seen in CESM2 (compared with CCSM4), which has the same land surface model as NorESM2-LM (CLM5 ).

The CMIP6 models with the lowest global total values for soil carbon do not always correspond to the lowest northern latitude values for soil carbon. For example, UKESM1-0-LL global total soil carbon is amongst the highest global totals seen in CMIP6; however, low quantities of soil carbon are seen in the northern latitudes (approximately 10 % of the global total). Conversely, BCC-CSM2-MR, CESM2, GFDL-ESM4, and NorESM2-LM have approximately 30 % of their global total stocks in the northern latitude region, which is consistent with the ratio seen in the benchmark dataset. This result suggests that representing global total soil carbon stocks consistently with the benchmark soil carbon does not imply consistency in the representation of northern latitude soil carbon stocks, and these should be evaluated separately. However, the large uncertainties associated with the empirical datasets for the northern latitudes are noted (Table ).

3.1.3 Spatial evaluation

A lack of consistency in the simulation of soil carbon was found amongst the CMIP5 models, which can be seen in Fig. a, where differences between the empirical and modelled data are shown. Northern latitude soil carbon was found to be underestimated in CMIP5, where areas of blue can be seen in the northern latitudes of the CMIP5 soil carbon map in Fig. a. This underestimation of CMIP5 northern latitude soil carbon is accompanied by significant overestimations seen in mid-latitude soil carbon. Specifically, large quantities of soil carbon which are inconsistent with our benchmark dataset can be seen in the mid-latitude regions in the following CMIP5 models: CanESM2, GFDL-ESM2G, GISS-E2-R, MIROC-ESM, and MPI-ESM-LR. Less significant overestimations are seen in HadGEM2-ES and IPSL-CM5A-LR (Fig. ). Systematic errors remain in the CMIP6 models; however, there are some improvements seen in the spatial simulation of soil carbon from CMIP5. Soil carbon is still underestimated in the northern latitudes, where the areas of blue still remain in the northern latitudes of the CMIP6 soil carbon map in Fig. a, though regions of overestimations in the northern latitudes are also seen amongst the CMIP6 models in BCC-CSM2-MR, CESM2, CNRM-ESM2-1, and NorESM2-LM (Fig. ), but it is noted that this representation might be more consistent with observations if a dataset including deeper soil carbon stocks was considered. CMIP6 shows improvements in the representation of mid-latitude soil carbon, where less of an overestimation is seen in CMIP6 compared with CMIP5 (Fig. a). This overestimation can still be seen in 4 of the 11 CMIP6 models: ACCESS-ESM1.5, CanESM2, MIROC-ES2L, and UKESM1-0-LL; however, the overestimations in CMIP6 are less inconsistent than when compared with CMIP5, and the number of models showing this limitation in CMIP6 has been reduced (Fig. ).

Figure 1

Maps presenting the difference between the modelled and benchmark data for the CMIP5 and CMIP6 ensembles, for (a) $C_{\mathrm{s}}$ (kg m${}^{-\mathrm{2}}$), (b) NPP (kg m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$), and (c) ${\mathit{\tau }}_{\mathrm{s}}$ (years).

[Figure omitted. See PDF]

Figure 2

Maps of the difference in soil carbon ($C_{\mathrm{s}}$) between the historical simulation of each CMIP6 model and the benchmark data.

[Figure omitted. See PDF]

Despite the differences seen in the spatial representation of soil carbon between the individual models in CMIP6, the ensemble mean has more areas of agreement within the ensemble compared with the ensemble mean in CMIP5. This can be seen in Fig. a, where there is less hatching (where hatched shaded areas represent regions of low agreement amongst the models in the ensemble; see Methods) in the CMIP6 map compared with the CMIP5 map. Specifically, ensemble mean soil carbon in CMIP6 has more areas of agreement in the mid-latitude region compared with the CMIP5 ensemble mean, where significant areas of disagreement are seen. This disagreement is likely due to the overestimation which exists in some of the CMIP5 models (Fig. ). Also, a reduction in the area of disagreement is seen in the northern latitudes in CMIP6 compared with CMIP5; however, this remains the region where the most disagreement exists across the generations. It is noted that this is a measure of agreement within the ensemble and not between the models and empirical data and so is dependent on the choice of ensemble members (see Figs.  and for individual model maps).

Figure 3

Ensemble mean maps for (a) $C_{\mathrm{s}}$ (kg m${}^{-\mathrm{2}}$), (b) NPP (kg m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$), and (c) ${\mathit{\tau }}_{\mathrm{s}}$ (years), presented for the CMIP6 ensemble, the CMIP5 ensemble, and the benchmark datasets. The hatched areas are used to show regions of low agreement within the ensemble (SD / mean $>\mathrm{0.75}$) and where regions of low soil carbon ($<$ 5 kg m${}^{-\mathrm{2}}$) have been excluded. Equivalent maps for the individual CMIP6 and CMIP5 models are shown within the Appendix: see Figs.  and for $C_{\mathrm{s}}$, Figs.  and for NPP, and Figs.  and for ${\mathit{\tau }}_{\mathrm{s}}$, respectively.

[Figure omitted. See PDF]

The inconsistency of the simulation of spatial soil carbon in CMIP6 is further evaluated using the spatial standard deviations, the spatial Pearson correlation coefficients, and RMSEs (see Methods), where the Taylor diagram (Fig. a) presents all three statistical assessments. The spatial standard deviation for soil carbon is shown on the radial axis between the standard $x$ and $y$ axes in Fig. a. The range of spatial standard deviations amongst the CMIP6 models sees a slight reduction from the range amongst the CMIP5 models, though significant differences remain. The CMIP6 models CNRM-ESM2-1, MIROC-ES2L, and UKESM1-0-LL best match the spatial standard deviation derived from the benchmark dataset (Tables  and ). It is found that the spatial representation of modelled soil carbon in CMIP6 is poorly correlated with the empirical soil carbon, where the CMIP6 ensemble spatial correlation coefficient with the empirical data is found to be 0.250. The spatial correlation coefficients between the individual CMIP6 and CMIP5 models with the empirical data can also be seen in Fig. a, where the low spatial correlation coefficients are shown by the curved correlation axis. The lowest spatial correlation coefficients amongst the CMIP6 models were $r$ values of 0.104 in IPSL-CM6A-LR and 0.115 in UKESM1-0-LL. The CMIP6 model that was the most spatially consistent with the empirical data is CNRM-ESM2-1, with an $r$ value of 0.630. The CMIP6 ensemble sees a slight reduction in the RMSE compared with the CMIP5 ensemble, suggesting a slight improvement (Fig. a). Significant improvements in the RMSE are seen in MIROC-ES2L from MIROC-ESM and MPI-ESM1.2-LR from MPI-ESM-LR. These results suggest small improvements in the simulation of soil carbon across this CMIP generation; however, the low spatial correlation coefficients and variable RMSEs seen across the models in CMIP6 suggest that inconsistencies with the benchmark data remain.

Figure 4

Taylor diagrams showing the spatial standard deviation (shown by the radial axis between standard $x$ and $y$ axes), the Pearson correlation coefficients (shown by the curved correlation axis), and the RMSE (shown by the grey contours) for the ESMs in both CMIP5 and CMIP6 compared with the benchmark datasets, for (a) soil carbon ($C_{\mathrm{s}}$), (b) NPP, and (c) soil carbon turnover time (${\mathit{\tau }}_{\mathrm{s}}$).

[Figure omitted. See PDF]

Figure 5

Bar charts comparing the root mean squared errors (RMSEs) in CMIP6 and CMIP5 for (a) soil carbon ($C_{\mathrm{s}}$), (b) NPP, and (c) soil carbon turnover time (${\mathit{\tau }}_{\mathrm{s}}$).

[Figure omitted. See PDF]

3.2 Net primary productivity: improved in CMIP6 relative to CMIP5

3.2.1 Global total evaluation

Global total NPP amongst the CMIP6 models appears to be consistent with the benchmark dataset (Table ), where the CMIP6 ensemble mean for NPP is approximately 95 % of the benchmark mean. The CMIP6 ensemble mean global total NPP (53.0 $\pm$ 9.39 PgC yr${}^{-\mathrm{1}}$) is found to be slightly lower than the derived mean benchmark value; however, it is comfortably within the observational uncertainty range (56.6 $\pm$ 14.3 PgC yr${}^{-\mathrm{1}}$). The equivalent values for the CMIP5 models can be seen in Table , where the CMIP5 ensemble total is also found to be within the observational uncertainty range (56.3 $\pm$ 15.4 PgC yr${}^{-\mathrm{1}}$).

Table 6

Table presenting global carbon fluxes and turnover time values for the 11 CMIP6 models included in this study and the benchmark datasets, including global total NPP (PgC yr${}^{-\mathrm{1}}$) and the effective average soil carbon turnover time (years).

Earth system model	NPP (PgC yr${}^{-\mathrm{1}}$)	${\mathit{\tau }}_{\mathrm{s}}$ (years)
ACCESS-ESM1.5	45.6	19.0
BCC-CSM2-MR	51.2	34.1
CanESM5	75.5	18.1
CESM2 (cSoilAbove1m)	43.9	25.8
CESM2 (cSoil)	–	50.4
CNRM-ESM2-1	45.6	41.5
GFDL-ESM4	52.6	11.2
IPSL-CM6A-LR	46.4	14.6
MIROC-ES2L	59.1	24.5
MPI-ESM1.2-LR	58.9	15.4
NorESM2-LM (cSoilAbove1m)	43.5	24.0
NorESM2-LM (cSoil)	–	60.8
UKESM1-0-LL	60.8	28.1
Ensemble mean	53.0 $\pm$ 9.39	23.3 $\pm$ 8.59
Benchmark datasets	56.6 $\pm$ 14.3	27.0${}_{-\mathrm{11}}^{+\mathrm{27}}$

Table 7

Table presenting global carbon fluxes and turnover time values for the 10 CMIP5 models included in this study and the benchmark datasets, including global total NPP (PgC yr${}^{-\mathrm{1}}$) and the effective average soil carbon turnover time (years).

Earth system model	NPP (PgC yr${}^{-\mathrm{1}}$)	${\mathit{\tau }}_{\mathrm{s}}$ (years)
BNU-ESM	44.3	16.6
CCSM4	42.9	14.3
CanESM2	59.0	72.9
GFDL-ESM2G	74.4	57.3
GISS-E2-R	31.0	47.1
HadGEM2-ES	69.1	16.8
IPSL-CM5A-LR	76.6	19.4
MIROC-ESM	47.1	56.8
MPI-ESM-LR	73.5	42.7
NorESM1-M	45.0	34.5
Ensemble mean	56.3 $\pm$ 15.4	37.8 $\pm$ 19.7
Benchmark datasets	56.6 $\pm$ 14.3	27.0${}_{-\mathrm{11}}^{+\mathrm{27}}$

The standard deviation surrounding the CMIP5 ensemble mean is greater than in CMIP6. This reduced standard deviation in CMIP6 is because several of the models have a simulated global total NPP that more closely matches the benchmark NPP global total value compared with the previous CMIP5 generation: GFDL-ESM4 from GFDL-ESM2G, IPSL-CM6A-LR from IPSL-CM5A-LR, MIROC-ES2L from MIROC-ESM, MPI-ESM1.2-LR from MPI-ESM1-M, and UKESM1-0-LL from HadGEM2-ES. The majority of CMIP6 models see a reduction in NPP from the CMIP5 equivalent model, which in general reduces the overestimation of NPP that was seen in the CMIP5 models (Tables  and ). However, this was not the case for CanESM5 from CanESM2, which sees an increase in the magnitude of NPP from CMIP5 to CMIP6, resulting in a consequent overestimation compared with the benchmark data. A reduced range of modelled global total NPP values is also seen in CMIP6 from CMIP5, where the range is reduced from 48.5 PgC yr${}^{-\mathrm{1}}$ in CMIP5 to 32.7 PgC yr${}^{-\mathrm{1}}$ in CMIP6. These results suggest that, overall, the representation of carbon fluxes in CMIP6 ESMs is more consistent than in CMIP5.

3.2.2 Spatial evaluation

Modelled NPP in CMIP6 appears to be spatially more consistent with the empirical data than in CMIP5. This is suggested by Fig. b, where the difference between the modelled and benchmark NPP is shown for both CMIP5 and CMIP6. It can be seen in the CMIP5 map that NPP is overestimated in the tropical regions, specifically in Africa and South-East Asia, and the equivalent CMIP6 difference map shows a clear reduction in this overestimation. This tropical overestimation of NPP, prominent in CMIP5 (Fig. ), is still seen in the CMIP6 models: CanESM5, MPI-ESM1.2-LR, and UKESM1-0-LL. However, this is not seen in the CMIP6 ensemble mean as it is likely negated by underestimations seen in CESM2, CNRM-ESM2-1, and NorESM2-LM (Fig. ). CMIP6 also sees more consistency with the benchmark dataset in the northern and mid-latitude regions compared with CMIP5, where more white areas are seen in the CMIP6 map in Fig. b. An underestimation of NPP is seen in both CMIP5 and CMIP6 on the western side of South America, though unusually high NPP is seen in this region in the MODIS NPP dataset (Fig. ). Moreover, greater agreement amongst the models within CMIP6 is seen compared with the models in CMIP5. This can be seen in Fig. b, where less hatching representing areas of disagreement within the ensemble is seen in CMIP6 compared with CMIP5. Specifically, CMIP6 sees less hatching in the northern latitudes and in the Middle East and south-eastern Europe as well as regions in South America, southern Africa, and Australia (see Figs.  and  for individual model maps).

Figure 6

Maps of the difference in net primary production (NPP) between the historical simulation of each CMIP6 model and the benchmark dataset.

[Figure omitted. See PDF]

The improved empirical consistency of modelled NPP in CMIP6 is also suggested when further evaluated using the same spatial metrics as with soil carbon. Despite a small range remaining in the spatial standard deviations amongst the CMIP6 models (shown by the radial axis in Fig. b), robust improvements in the spatial correlation coefficients (shown by the curved axis in Fig. b) and RMSEs are seen across the ensemble compared with CMIP5 (Fig. b). Notable improvements in the representation of NPP are seen in GFDL-ESM4 compared with GFDL-ESM2G, IPSL-CM6A-LR compared with IPSL-CM5A-LR, and UKESM1-0-LL compared with HadGEM2-ES, with reduced RMSEs seen in each updated model. A general improvement in the spatial correlation coefficients is seen across all the CMIP6 models, where the circle markers (CMIP6 models) in Fig. b have higher correlation values than the cross markers (CMIP5 models). The general improvement has resulted in the CMIP6 ensemble correlation coefficient (0.836) being greater compared with the equivalent CMIP5 value (0.711). The lowest correlations between modelled and observed NPP amongst the CMIP5 models are GISS-E2-R (0.274) and CanESM2 (0.469). The updated version CanESM5 remains the lowest correlation seen in CMIP6 (0.655); however, an improvement in the correlation is seen. The updated version of the GISS model is not included in the CMIP6 ensemble considered in this study, which could be a reason for the increased ensemble mean correlation. However, this effect does not take away from the improvements seen across the CMIP6 models. HadGEM2-ES (0.764) and MPI-ESM-LR (0.764) were the CMIP5 models with the highest correlation with the benchmark NPP, and the updated CMIP6 equivalents of these models remain the models with the greatest correlations, but again improvements in the correlations are seen (0.816 in UKESM1-0-LL and 0.785 in MPI-ESM1.2-LR).

3.3 Soil carbon turnover time: no major improvements in CMIP6 compared with CMIP5

3.3.1 Global evaluation

There are minor improvements suggested in the simulated effective global ${\mathit{\tau }}_{\mathrm{s}}$ amongst select CMIP6 models (Table ) compared with CMIP5 (Table ). The ensemble mean effective global ${\mathit{\tau }}_{\mathrm{s}}$ was overestimated in CMIP5 (37.8 $\pm$ 19.7 years) when compared with the derived mean ${\mathit{\tau }}_{\mathrm{s}}$ using the benchmark datasets (27.0${}_{-\mathrm{11}}^{+\mathrm{27}}$ years), which is reduced to a less significant underestimation in CMIP6 (23.3 $\pm$ 8.59 years), though both the CMIP5 and CMIP6 estimates fall within the observational uncertainty range. The associated ensemble spread in effective mean ${\mathit{\tau }}_{\mathrm{s}}$ is less in CMIP6 compared with CMIP5, with an ensemble standard deviation of approximately 50 % less. A significant range is seen in the effective global ${\mathit{\tau }}_{\mathrm{s}}$ values amongst the CMIP5 models, with a 5-fold difference between the lowest and highest values (Table ). This range is mostly due to large overestimations seen amongst the CMIP5 models, for example in CanESM2, GFDL-ESM2G, and MIROC-ESM. A reduced range is seen amongst the models in CMIP6; however, a 4-fold range still exists between the lowest and highest values (Table ). This reduced range is partly due to reductions in the effective global ${\mathit{\tau }}_{\mathrm{s}}$ values in CMIP6 models compared with the equivalent model in CMIP5, specifically, CanESM5 from CanESM2, GFDL-ESM4 from GFDL-ESM2G, MIROC-ES2L from MIROC-ESM, and MPI-ESM1.2-LR from MPI-ESM-LR, though overestimations do remain in CMIP6, for example in CNRM-ESM2-1, where the slowest effective turnover time was seen. Moreover, the range is also reduced due to improvements seen in models which underestimated ${\mathit{\tau }}_{\mathrm{s}}$ in CMIP5, such as UKESM1-0-LL from HadGEM2-ES and CESM2 from CCSM4.

3.3.2 Spatial evaluation

The comparison of spatial ${\mathit{\tau }}_{\mathrm{s}}$ in CMIP6 with CMIP5 has more varied results compared with simulated NPP. The CMIP5 ensemble showed an underestimation of ${\mathit{\tau }}_{\mathrm{s}}$ in the northern latitudes, which is replaced with an overestimation of ${\mathit{\tau }}_{\mathrm{s}}$ in CMIP6 when compared with the benchmark data (Fig. c). This northern latitude overestimation in the CMIP6 ensemble is a result of the overestimations of ${\mathit{\tau }}_{\mathrm{s}}$ in CESM2 and NorESM2-LM (Fig. ), which dominate in the CMIP6 ensemble mean. It is noted that this result may differ if deeper soil carbon stocks were considered. The northern latitude underestimation of ${\mathit{\tau }}_{\mathrm{s}}$ is still seen within the CMIP6 models: CanESM5, CNRM-ESM2-1, GFDL-ESM4, IPSL-CM6A-LR, MIROC-ES2L, MPI-ESM1.2-LR, and UKESM1-0-LL (Fig. ). An overestimation of mid-latitude ${\mathit{\tau }}_{\mathrm{s}}$ was seen in the CMIP5 models MIROC-ESM and MPI-ESM-LR (Fig. ), which is no longer seen in the equivalent updated CMIP6 models MIROC-ES2L and MPI-ESM1-2-LR. However, an overestimation of mid-latitude ${\mathit{\tau }}_{\mathrm{s}}$ is seen in CMIP6 models BCC-CSM2-MR, CNRM-ESM2-1, and UKESM1-0-LL (Fig. ). The uncertainty in simulated northern latitude ${\mathit{\tau }}_{\mathrm{s}}$ is also apparent in Fig. c, where the hatching shows the lack of agreement within the CMIP6 ensemble in this region. However, more agreement within the CMIP6 ensemble is seen in the same figure in the mid-latitudes and in the tropical regions compared with CMIP5 (see Figs.  and for individual model maps).

Figure 7

Maps of the difference in soil carbon turnover time (${\mathit{\tau }}_{\mathrm{s}}$) between the historical simulation of each CMIP6 model and the benchmark datasets.

[Figure omitted. See PDF]

The simulation of spatial ${\mathit{\tau }}_{\mathrm{s}}$ in CMIP6 is further evaluated against the empirical data with the additional statistical metrics. Modelled ${\mathit{\tau }}_{\mathrm{s}}$ is found to be poorly spatially correlated with empirical ${\mathit{\tau }}_{\mathrm{s}}$ in both the CMIP5 and CMIP6 models (shown by the curved axis in Fig. c). A slight increase in the ensemble mean spatial correlations is seen from CMIP5 (0.188) to CMIP6 (0.267) due to increases seen amongst individual models between CMIP5 and CMIP6: CESM2 from CCSM4, MPI-ESM1.2-LR from MPI-ESM-LR, and NorESM2-LM from NorESM1-M. However, the consistency of modelled ${\mathit{\tau }}_{\mathrm{s}}$ with the benchmark datasets remains low. A particularly large range is seen in the spatial standard deviations of ${\mathit{\tau }}_{\mathrm{s}}$ amongst the CMIP6 models, which is an increased range from CMIP5 (shown by the radial axis in Fig. c). The CMIP6 models with the most extreme overestimations of the spatial standard deviations compared with the derived benchmark value (NorESM2-LM, CESM2, and ACCESS-ESM1.5) are also found to have large RMSEs (Fig. c). Amongst the remaining CMIP6 models, the RMSEs for modelled ${\mathit{\tau }}_{\mathrm{s}}$ remain relatively consistent between CMIP5 and CMIP6.

3.4 Drivers of soil carbon spatial patterns: soil carbon spuriously highly correlated with NPP in CMIP5 and CMIP6

3.4.1 Global drivers

A negligible correlation ($\approx$ 0) is found between the benchmark estimates of soil carbon and NPP, suggesting that soil carbon is not spatially correlated with NPP in the real world. On the other hand, soil carbon and NPP ($C_{\mathrm{s}}$–NPP) were found to be significantly correlated in the models in both CMIP5 and CMIP6. The $C_{\mathrm{s}}$–NPP spatial correlation was found to be greater than 0.5 for 6 out of the 10 CMIP5 ESMs and 8 out of the 11 models in CMIP6 (Fig. a). However, a low spatial correlation is found in the CMIP6 models CESM2 (0.134), NorESM2-LM (0.261), and BCC-CSM2-MR (0.214), values most consistent with the benchmark datasets. The $C_{\mathrm{s}}$–${\mathit{\tau }}_{\mathrm{s}}$ spatial correlations found in the CMIP6 models tend to underestimate the positive correlation seen in the benchmark datasets (Fig. a). The majority of CMIP6 models see a negligible or slightly negative $C_{\mathrm{s}}$–${\mathit{\tau }}_{\mathrm{s}}$ spatial correlation despite a low positive correlation produced by the benchmark datasets. The models BCC-CSM2-MR, MIROC-ES2L, and NorESM2-LM are most consistent with the benchmark $C_{\mathrm{s}}$–${\mathit{\tau }}_{\mathrm{s}}$ correlation.

Figure 8

Scatter plots investigating the relationships between different Pearson correlation coefficients of climate variables, (a) $C_{\mathrm{s}}$–${\mathit{\tau }}_{\mathrm{s}}$ against $C_{\mathrm{s}}$–NPP, (b) NPP–$T$ against NPP–$\mathit{\theta }$, and (c) ${\mathit{\tau }}_{\mathrm{s}}$–$T$ against ${\mathit{\tau }}_{\mathrm{s}}$–$\mathit{\theta }$.

[Figure omitted. See PDF]

The modelled NPP to temperature (NPP–$T$) spatial correlations in CMIP6 are consistent with the positive relationship seen in the benchmark datasets; however, the magnitude of this positive correlation varies amongst the models (Fig. b). The magnitude of the positive NPP–$T$ correlation is underestimated in CanESM5, GFDL-ESM4, and NorESM2-LM but is otherwise relatively consistent amongst the CMIP6 models. Nonetheless, a much greater range in the modelled NPP–$T$ correlations was seen amongst the CMIP5 models, suggesting an improved representation of this relationship in CMIP6. The variation in modelled NPP–$\mathit{\theta }$ correlations remains in CMIP6, with models disagreeing on the sign and magnitude of the correlation of NPP with soil moisture. The modelled NPP–$\mathit{\theta }$ correlation is the most consistent with the benchmark correlations in GFDL-ESM4, MPI-ESM1.2-LR, and UKESM1-0-LL (Fig. b).

It is generally agreed across the models in CMIP6 and CMIP5 that ${\mathit{\tau }}_{\mathrm{s}}$ and temperature ($T$) are negatively correlated, with the exception of MPI-ESM1.2-LR, where a slight positive correlation is seen (Fig. c). This is consistent with the negative ${\mathit{\tau }}_{\mathrm{s}}$–$T$ correlation derived with the benchmark dataset. There is variation amongst the models in the magnitude of the negative correlation, with a significant overestimation seen in CanESM5. A negative correlation is also seen in the ${\mathit{\tau }}_{\mathrm{s}}$–$\mathit{\theta }$ correlation derived with the benchmark datasets. Inconsistencies with this empirical relationship are seen amongst the models in both CMIP5 and CMIP6, with many negligible and positive correlations deduced (Fig. c). The exception is again MPI-ESM1.2-LR, which in this case is the model most consistent with the benchmark ${\mathit{\tau }}_{\mathrm{s}}$–$\mathit{\theta }$ correlation.

3.4.2 Regional drivers

The spatial correlations of modelled $C_{\mathrm{s}}$–NPP are shown to be overestimated at every latitude in both CMIP6 and CMIP5 compared with the equivalent correlations derived from the empirical datasets. It can be seen that the CMIP6 ensemble mean $C_{\mathrm{s}}$–NPP correlation has an even larger positive bias compared with the benchmark correlation than in CMIP5. The empirical data see a reduced $C_{\mathrm{s}}$–NPP correlation in the northern latitudes, whereas a slight but less significant reduction is seen in the models (Fig. a). The spatial correlation between $C_{\mathrm{s}}$ and ${\mathit{\tau }}_{\mathrm{s}}$ is shown to vary against latitude in the empirical datasets, where a greater correlation is seen in the tropical and northern latitude regions, and a negligible correlation is seen at the mid-latitudes (Fig. b). The CMIP6 models simulate the negligible $C_{\mathrm{s}}$ and ${\mathit{\tau }}_{\mathrm{s}}$ seen at the mid-latitudes relatively consistently with the benchmark data, where an improved consistency is seen from CMIP5. However, the CMIP6 models do not simulate the tropical and northern latitude positive $C_{\mathrm{s}}$–${\mathit{\tau }}_{\mathrm{s}}$ correlations, where a negligible modelled correlation remains in these regions. CMIP5 is more consistent with the benchmark correlations than in CMIP6, where a positive modelled correlation $C_{\mathrm{s}}$–${\mathit{\tau }}_{\mathrm{s}}$ is seen (Fig. b).

Figure 9

The latitudinal profiles of the Pearson correlation coefficients between soil carbon and (a) NPP ($C_{\mathrm{s}}$–NPP), (b) soil carbon turnover time ($C_{\mathrm{s}}$–${\mathit{\tau }}_{\mathrm{s}}$), (c) soil moisture ($C_{\mathrm{s}}$–$\mathit{\theta }$), and (d) temperature ($C_{\mathrm{s}}$–$T$).

[Figure omitted. See PDF]

The spatial correlation between modelled soil carbon and soil moisture ($C_{\mathrm{s}}$–$\mathit{\theta }$) is consistent with the correlations seen in the benchmark datasets at every latitude, with an improvement seen in the tropical correlation patterns in CMIP6 compared with CMIP5 (Fig. c). Both the CMIP5 and CMIP6 ensembles span the benchmark $C_{\mathrm{s}}$–$\mathit{\theta }$ correlation, though large model ranges in the $C_{\mathrm{s}}$–$\mathit{\theta }$ sensitivity are seen across all latitudes. However, there is a reduced ensemble spread in the $C_{\mathrm{s}}$–$\mathit{\theta }$ correlation from CMIP5 to CMIP6 at low- and mid-latitudes. An overestimation of the negative $C_{\mathrm{s}}$–$T$ correlation seen in the benchmark datasets is present in both the CMIP5 and CMIP6 models, except the high latitudes (Fig. d). This modelled $C_{\mathrm{s}}$–$T$ correlation is particularly underestimated at the lower tropical latitudes, where a greater positive correlation is seen here in the benchmark datasets. Figure d suggests a slight improvement in the modelled tropical $C_{\mathrm{s}}$–$T$ correlation in CMIP6 and a worsening of modelled $C_{\mathrm{s}}$–$T$ at the high latitudes than in CMIP5 when compared with the $C_{\mathrm{s}}$–$T$ correlations in the benchmark datasets.

4 Discussion

4.1 Soil carbon stocks

4.1.1 Global total soil carbon

Simulating global soil carbon stocks that are consistent with empirical data is required to produce reliable projections of future soil carbon storage and emission . This study deduces a CMIP6 ensemble mean global total soil carbon of 1206 $\pm$ 445 PgC (Table ), using regridded model resolutions (see Methods). It is noted that states a CMIP6 ensemble of 1553 $\pm$ 672 PgC; however, the full soil carbon profile is considered for CESM2 and NorESM2-LM, as opposed to a depth of 1 m considered in this study. Additionally, this study deduces a comparable CMIP5 ensemble mean global soil carbon value of 1480 $\pm$ 810 PgC (Table ) using equivalent dates in the historical simulation (1950–2000). state an ensemble mean soil carbon value of 1520 $\pm$ 770 PgC in CMIP5; however, the study includes the models BCC-CSM1.1, CESM1-CAM5, and INM-CM4, which are missing from the analysis in this study due to data availability. present a CMIP5 ensemble mean soil carbon value of 1502 $\pm$ 798 PgC, but this calculation includes multiple model versions (for example, LR and MR) from the same modelling centre in their ensemble. A caveat of this evaluation study is the non-independent nature of CMIP ESMs, where for example CESM2 and NorESM2-LM share the same LSM. Additionally, the ensembles included here do not necessarily represent all models that exist within each CMIP generation. However, the evaluation completed here allows for general improvements in the simulation of soil carbon stocks and fluxes between the CMIP5 and CMIP6 generations to be noted and key areas for future model development to be highlighted.

Despite a suggestion of a reduced spread in model estimates of global total soil carbon within CMIP6 relative to CMIP5, discrepancies remain in the consistency of these estimates with the observations between the two CMIP generations. It should also be noted that CMIP6 does not simply contain updated versions of every model in CMIP5: some new models are included and some CMIP5 models not included in CMIP6. These factors together with the uncertainty associated with empirical datasets have resulted in no robust conclusion being drawn on the improvement of soil carbon simulation in CMIP6 compared with CMIP5. Due to the potential significant feedback that exists between soil carbon and global climate, this lack of consistency reduces our confidence in future projections of climate change .

4.1.2 Spatial soil carbon patterns

Modelled soil carbon was found to be poorly spatially correlated with the empirical data amongst models in both CMIP5 and CMIP6 (Fig. a). An improvement in CMIP6 ESMs was seen in the spatial patterns across the mid-latitudes, which were generally overestimated in CMIP5. However, significant underestimations of modelled soil carbon in the northern latitudes still remain, which have a significant impact on model predictions of global total soil carbon stocks (Fig. a). This systematic underestimation was previously reported in the literature as a limitation of the CMIP5 models, where found northern latitude soil carbon to be less consistent with the empirical data than on a global scale. This limitation remains amongst models in the CMIP6 generation, where it was found that the majority of CMIP6 models underestimate northern latitude soil carbon stocks regardless of whether the global soil carbon stocks are underestimated.

However, an exception to this northern latitude underestimation is seen within CMIP6 in the models CESM2 and NorESM2-LM. These ESMs include the land surface model (LSM) CLM5 , which is the first LSM to include the representation of vertically resolved soil carbon in its CMIP simulations. This representation enables the inclusion of separate carbon pools at varying depths in the soil, which aims to more consistently simulate soil carbon with the real world . This is of particular importance in the northern latitudes, where carbon stocks are expected to exist at much greater depths than the 1 m considered in this study . This can be seen in Table , where increased magnitudes of soil carbon stocks are shown when increased depths are considered using the empirical datasets. A more thorough evaluation of soil carbon in both CESM2 and NorESM2-LM is suggested for future research, with a particular focus on this improved northern latitude soil carbon stock simulation; however, this evaluation of deeper soil carbon stocks (below 1 m) is beyond the scope of this study.

Accurately simulating soil carbon in the northern latitude regions is of particular importance as it is a major part of the total global soil carbon pool . Additionally, much of the carbon stored in these soils is held within permafrost, which is known to be particularly sensitive to climate change. Permafrost thaw under climate change has the potential to release significant amounts of carbon into the atmosphere over a short period of time with increased warming , representing a positive feedback within the climate system. Permafrost dynamics are generally poorly represented in ESMs, where found CMIP6 ESMs to have a similar representation compared with CMIP5. Underestimating soil carbon in the northern latitudes may result in underestimation of the impact of this feedback in future climate change projections. Future improvements are needed to improve the simulation of soil carbon stocks globally, but particularly in the northern latitudes.

4.2 Drivers of soil carbon change

To allow for a more in-depth understanding of the inconsistencies found between modelled and empirical soil carbon, the simulation of above- and below-ground controls of soil carbon was also evaluated. Simulations of contemporary soil carbon can be disaggregated into the effects of litterfall, which is well approximated by plant NPP, and effective soil carbon turnover time (${\mathit{\tau }}_{\mathrm{s}}$), which is affected by both the temperature and moisture of the soil . If models are to reliably simulate soil carbon in a way that is consistent with empirical data, the spatial drivers of soil carbon, NPP and ${\mathit{\tau }}_{\mathrm{s}}$, must also be simulated consistently with empirical data. Isolating the effects of NPP and ${\mathit{\tau }}_{\mathrm{s}}$ on soil carbon helps us to break down the simulation of soil carbon to help understand the limitations and inconsistencies seen amongst the models.

4.2.1 NPP

An improved simulation of NPP is suggested in the ESMs included from CMIP6 compared with the ESMs from CMIP5. This conclusion is suggested by an increased number of models in our CMIP6 ensemble having global total NPP values consistent with empirical data (Table ), the overestimation of tropical NPP amongst CMIP5 models being seen to be reduced amongst the CMIP6 models (Fig. b), and more agreement being seen within CMIP6 relative to CMIP5 in the simulation of mid-latitude and northern latitude NPP (Fig. b). Modelled NPP was found to be robustly more consistent with the empirical data in our CMIP6 ensemble compared with the CMIP5 ensemble in all statistical evaluation metrics. Since CMIP5, multiple models have seen an addition of a dynamic nitrogen cycle , where the models with nitrogen cycles are highlighted in Fig.  by the shaded bars. The results suggest an improvement in the simulation of NPP with the addition of dynamic nitrogen in models. However, CMIP6 models that do not represent a nitrogen cycle also mostly see improvements in the simulation of NPP, suggesting NPP is more constrained by observations in the newest generation of models. CanESM5 is the only ESM within CMIP6 included here to not see an overall improvement in the simulation of NPP, where NPP is found to be overestimated compared with the benchmark dataset. It is likely that the inclusion of a nitrogen cycle in this model would limit this overestimated NPP and improve consistency with the observations .

Despite this apparent improved simulation of NPP in CMIP6, the spatial correlation between modelled soil carbon and NPP was found to be inconsistent with the equivalent empirically derived relationship. This result was previously shown for the CMIP5 models and has been more recently shown for the CMIP6 models , both agreeing with the results found here. The majority of CMIP6 models were found to have positive $C_{\mathrm{s}}$–NPP spatial correlations, as opposed to a negligible spatial correlation found in the observations (Fig. a). Despite NPP driving the spatial pattern of soil carbon stocks due to carbon input from vegetation, a positive correlation was not expected in the real world due to regions with high soil carbon not correlating with regions of high NPP. For example, in the observationally derived data, soil carbon stocks are greatest in the northern latitudes due to long turnover times in these regions, whereas NPP is lower due to cold temperatures in these regions limiting vegetation growth. The three CMIP6 models which did not significantly overestimate this correlation (CESM2, NorESM2-LM, and BCC-CSM2-MR) are three of the models with the most empirically consistent proportion of soil carbon stocks in the northern latitudes. Conversely, the tropical regions see high NPP values, but warmer temperatures result in faster turnover times and lower soil carbon stocks. NPP is expected to increase in the future under climate change , which means an overly positive correlation in models could result in a subsequent increase in modelled projections of soil carbon stocks. An overestimation of future soil carbon storage could result in an overestimation of the future carbon sink and an inaccurate global carbon budget .

4.2.2 Soil carbon turnover time

The systematic improvements suggested from the evaluation of NPP simulation within our CMIP6 ensemble are not suggested for the simulation of ${\mathit{\tau }}_{\mathrm{s}}$, where the simulation of ${\mathit{\tau }}_{\mathrm{s}}$ appears to remain inconsistent with the empirical data in CMIP6 from CMIP5. Improvements are suggested within CMIP6 relative to CMIP5, such as more agreement within the ensemble at the mid-latitudes and in the tropical regions; however, less agreement is seen in the northern latitudes (Fig. c). Northern latitude ${\mathit{\tau }}_{\mathrm{s}}$ is generally underestimated in models, which corresponds to the underestimation of soil carbon seen in these regions. This has been previously identified in ESMs, where it was found that the underestimation of global ${\mathit{\tau }}_{\mathrm{s}}$ amongst the CMIP5 models is primarily due to low values in the northern latitudes . The reduced agreement in CMIP6 is due to long ${\mathit{\tau }}_{\mathrm{s}}$ values existing in the northern latitudes of CESM2 and NorESM2-LM alongside the general ensemble underestimations (Fig. ). The increased northern latitude ${\mathit{\tau }}_{\mathrm{s}}$ values in CESM2 and NorESM2-LM are likely to be due to the representation of vertically resolved soil carbon pools, which allows for differential ${\mathit{\tau }}_{\mathrm{s}}$ values for pools at varying depths. Despite these individual improvements since CMIP5, large discrepancies exist within the CMIP6 ensemble between modelled and empirical ${\mathit{\tau }}_{\mathrm{s}}$.

To simulate ${\mathit{\tau }}_{\mathrm{s}}$ consistently with observations, the relationship of ${\mathit{\tau }}_{\mathrm{s}}$ with both temperature ($T$) and moisture ($\mathit{\theta }$) must also be simulated in a way that is consistent with observations. Generally, the ${\mathit{\tau }}_{\mathrm{s}}$–$T$ relationship is consistently simulated; however, there is variation in the modelled temperature sensitivity of ${\mathit{\tau }}_{\mathrm{s}}$ across the ensemble. The ${\mathit{\tau }}_{\mathrm{s}}$–$\mathit{\theta }$ relationship is less consistently represented, where the majority of CMIP6 models do not match the empirically derived relationship. Despite a positive dependence of soil respiration on soil moisture in the empirical data, many of the CMIP6 models display a contradictory positive ${\mathit{\tau }}_{\mathrm{s}}$–$\mathit{\theta }$ correlation (Fig. ). Many of the models use functions that increase respiration with soil moisture (see Sect. ), so the increase in ${\mathit{\tau }}_{\mathrm{s}}$ with increasing soil moisture indicated by positive ${\mathit{\tau }}_{\mathrm{s}}$–$\mathit{\theta }$ correlations in the models is unexpected. We note that this effect occurs most strongly in the models with a very strongly negative ${\mathit{\tau }}_{\mathrm{s}}$–$T$ relationship (Fig. c), so it could in fact be an artefact of a negative correlation between temperature and soil moisture. In this context it is also important to consider what soil moisture in LSMs represents . The aim within the models is to act as the lower boundary condition for atmospheric models, and therefore their soil parameters may historically have been tuned to give appropriate evaporation rates and not necessarily to represent the soil moisture itself in an accurate way, so it may be more relevant to consider the large-scale emergent patterns of ${\mathit{\tau }}_{\mathrm{s}}$ than the direct relationships between soil moisture and respiration. It is noted that the empirical relationship shows ${\mathit{\tau }}_{\mathrm{s}}$ reducing with higher soil moisture, which suggests that the observations are picking up more on longer turnover times in dry areas rather than in saturated areas such as peatlands. This may be due to having only surface soil moisture information, whereas peatlands, while saturated at depth, typically have a water table $\sim$ 10 cm below the surface and can be very dry at the surface . Thus, while models do not include the necessary processes for peat formation , this is unlikely to be the cause of the discrepancy since it would lead to even more of a positive ${\mathit{\tau }}_{\mathrm{s}}$–$\mathit{\theta }$ correlation in the models.

Different processes control soil carbon formation in different ecosystems, including stabilization by clay particles, transformation by microbes, and nitrogen and phosphorous availability . In the present study, the largest discrepancies in both soil carbon and turnover times are seen in permafrost and peatland areas (see Figs.  and ). For example, the western Siberian peatland complex stands out in the majority of the panels in these figures as an area of high model error. This is partly because the soil carbon turnover times and quantities of soil carbon are largest in these regions but also partly because of the specific controlling processes in these ecosystems. A key part of soil carbon development in permafrost regions is the fact that organic material can be preserved in frozen soil, including via cryoturbation and yedoma deposits, which have not yet been thoroughly represented in models . There are a variety of other factors, such as plants storing significantly more of their carbon below ground instead of above ground in cold climates and recalcitrant vegetation such as mosses, which are not represented in most ESMs . Peatland formation is controlled primarily by waterlogging, which reduces the oxygen available for decomposition, but there are a huge number of additional physical and biogeochemical feedbacks that take place . These kinds of small-scale processes and inhomogeneities are difficult to resolve in global models with $\sim$ 100 km${}^{\mathrm{2}}$ grid cells, and this should be weighed up against their relative impact on global carbon budgets when considering including these processes in ESMs. However, it is suggested that the large-scale discrepancies such as in the permafrost and large peatland areas can and should be resolved in future model versions.

Our results suggest that much of the uncertainty associated with modelled soil carbon stocks can be attributed to the simulation of below-ground processes. The apparent improved consistency of NPP with empirical data suggests that considerable efforts have been made to achieve an improved representation of above-ground processes in CMIP6 ESMs since the release of the CMIP5 ensemble. However, the same improvements are not apparent in the simulation of ${\mathit{\tau }}_{\mathrm{s}}$ as systemic limitations remain in the new generation of ESMs considered in this study, suggesting that the same progress in the model development of below-ground processes has not been achieved between CMIP5 and CMIP6. Moreover, focus on above-ground processes without consideration of below-ground processes can result in inconsistencies of soil carbon stocks. For example, the inclusion of a nitrogen cycle has been shown to lead to a reduction in soil carbon in the model (see Fig. 6 in ), so tuning of the baseline turnover rates is required to keep soil carbon stocks consistent with observed values.

The required improvement of soil carbon pool turnover rates has previously been identified for the CMIP5 ensemble , and more recently, found that the difference in turnover times amongst the CMIP6 models is responsible for approximately 88 % of the variation seen in global soil carbon stocks amongst the models and stated that constraining key parameters which control soil carbon turnover processes is a key area for future model development. A key development seen in CMIP6 since CMIP5 is the representation of vertically resolved soil carbon. Models which simulate non-vertically resolved soil carbon typically turn over all the carbon based on the temperature near the soil surface. This could lead to reduced quantities of soil carbon and an underestimation of northern latitude soil carbon stocks, due to near-surface soil being warmer than the deeper soil and as turnover is known to respond exponentially to temperature . Overall, further improvements in the representation of soil carbon turnover time, with a particular focus on the northern latitudes, are identified as a key area for future model development.

5 Conclusions

The ability of Earth system models (ESMs) to simulate present-day soil carbon is vital to help predict reliable global carbon budget estimates, which are required for Paris Agreement targets. In this study, CMIP6 ESMs have been evaluated against empirical datasets to assess their ability to represent soil carbon and related controls: net primary productivity (NPP) and the effective soil carbon turnover time (${\mathit{\tau }}_{\mathrm{s}}=C_{\mathrm{s}}$ / $R_{\mathrm{h}}$). The evaluation is completed by comparison with the previous generation of CMIP5 ESMs to assess where improvements have been made and to identify priorities for future model development. Below the key conclusions from this study are listed.

• The spatial patterns of soil carbon in CMIP6 models appear to be more in agreement with each other than they were in CMIP5 and are more consistent with observations in the mid-latitudes, although caveats around the uncertainty in observations and the ensemble design make this conclusion uncertain. However, soil carbon is still heavily underestimated in high northern latitudes (with the exception of the two CMIP6 models that represent deep soil carbon).

• Overall, we are not able to identify significant improvements in the simulation of the observed spatial pattern of soil carbon across the globe from the CMIP5 to CMIP6 generation.

• There is good evidence that spatial patterns of contemporary NPP are better simulated in CMIP6 than in CMIP5 generation models when compared with satellite-derived estimates.

• However, spatial patterns of ${\mathit{\tau }}_{\mathrm{s}}$ continue to be poorly represented in CMIP6 models, in comparison with estimates derived from observational datasets of soil carbon and heterotrophic respiration ($R_{\mathrm{h}}$).

• Importantly, soil carbon simulations in both the CMIP5 and CMIP6 ESM generations seem to be spuriously highly correlated with NPP, which may make soil carbon in these models over-responsive to future projected changes in NPP.

Taken together, these conclusions point to a need for a much greater emphasis on improving the representation of below-ground soil processes in the next generation (CMIP7) of ESMs.

Appendix A

Figure A1

Maps comparing empirical datasets of soil carbon ($C_{\mathrm{s}}$). The benchmark dataset is a map plot showing $C_{\mathrm{s}}$ approximated to a depth of 1 m by combining the Harmonized World Soils Database (HWSD) and Northern Circumpolar Soil Carbon Database (NCSCD) , where the NCSCD was used where overlap occurs. Additional map plots are shown for empirical $C_{\mathrm{s}}$ estimated by the World Inventory of Soil property Estimates (WISE30sec) , the named “S2017” from , the Global Soil Dataset for use in Earth system models (GSDE) , and the Global Gridded Surfaces of Selected Soil Characteristics (IGBP-DIS) .

[Figure omitted. See PDF]

Figure A2

Maps of empirical carbon flux datasets. Net primary production (NPP) is approximated using the MODIS NPP dataset , and heterotrophic respiration ($R_{\mathrm{h}}$) is approximated using the CARDAMOM $R_{\mathrm{h}}$ dataset .

[Figure omitted. See PDF]

Figure A3

Maps of the difference in soil carbon ($C_{\mathrm{s}}$) between the historical simulation (1950–2000) for the CMIP5 models and the benchmark dataset.

[Figure omitted. See PDF]

Figure A4

Maps of the difference in NPP between the historical simulation (1995–2005) for the CMIP5 models and the benchmark dataset.

[Figure omitted. See PDF]

Figure A5

Maps of the difference in ${\mathit{\tau }}_{\mathrm{s}}$ between the historical simulation for the CMIP5 models and the benchmark datasets, where ${\mathit{\tau }}_{\mathrm{s}}$ is defined as the ratio of $C_{\mathrm{s}}$ (1950–2000) to $R_{\mathrm{h}}$ (1995–2005).

[Figure omitted. See PDF]

Figure A6

Maps of soil carbon ($C_{\mathrm{s}}$) in the historical simulation (1950–2000) for the CMIP6 models.

[Figure omitted. See PDF]

Figure A7

Maps of soil carbon ($C_{\mathrm{s}}$) in the historical simulation (1950–2000) for the CMIP5 models.

[Figure omitted. See PDF]

Figure A8

Maps of NPP in the historical simulation (1995–2005) for the CMIP6 models.

[Figure omitted. See PDF]

Figure A9

Maps of NPP in the historical simulation (1995–2005) for the CMIP5 models.

[Figure omitted. See PDF]

Figure A10

Maps of soil carbon turnover times (${\mathit{\tau }}_{\mathrm{s}}$) in the historical simulation for the CMIP6 models, where ${\mathit{\tau }}_{\mathrm{s}}$ is defined as the ratio of $C_{\mathrm{s}}$ (1950–2000) to $R_{\mathrm{h}}$ (1995–2005).

[Figure omitted. See PDF]

Figure A11

Maps of soil carbon turnover times (${\mathit{\tau }}_{\mathrm{s}}$) in the historical simulation for the CMIP5 models, where ${\mathit{\tau }}_{\mathrm{s}}$ is defined as the ratio of $C_{\mathrm{s}}$ (1950–2000) to $R_{\mathrm{h}}$ (1995–2005).

[Figure omitted. See PDF]

Table A1

Table presenting global carbon fluxes (PgC yr${}^{-\mathrm{1}}$), NPP, and $R_{\mathrm{h}}$ for the 11 CMIP6 models included in this study and the empirical benchmark datasets.

Earth system model	NPP (PgC yr${}^{-\mathrm{1}}$)	$R_{\mathrm{h}}$ (PgC yr${}^{-\mathrm{1}}$)
ACCESS-ESM1.5	45.6	45.1
BCC-CSM2-MR	51.2	48.9
CanESM5	75.5	75.0
CESM2	43.9	38.3
CNRM-ESM2-1	45.6	40.3
GFDL-ESM4	52.6	43.7
IPSL-CM6A-LR	46.4	39.9
MIROC-ES2L	59.1	52.7
MPI-ESM1.2-LR	58.9	53.4
NorESM2-LM	43.5	38.2
UKESM1-0-LL	60.8	57.5
Ensemble mean	53.0 $\pm$ 9.39	48.4 $\pm$ 10.5
Benchmark datasets	56.6 $\pm$ 14.3	51.7 $\pm$ 21.8

Table A2

Table presenting global carbon fluxes (PgC yr${}^{-\mathrm{1}}$), NPP, and $R_{\mathrm{h}}$ for the 10 CMIP5 models included in this study and the empirical benchmark datasets.

Earth system model	NPP (PgC yr${}^{-\mathrm{1}}$)	$R_{\mathrm{h}}$ (PgC yr${}^{-\mathrm{1}}$)
BNU-ESM	44.3	42.5
CCSM4	42.9	41.4
CanESM2	59.0	58.8
GFDL-ESM2G	74.4	62.7
GISS-E2-R	31.0	39.5
HadGEM2-ES	69.1	67.0
IPSL-CM5A-LR	76.6	62.4
MIROC-ESM	47.1	41.2
MPI-ESM-LR	73.5	59.9
NorESM1-M	45.0	41.3
Ensemble mean	56.3 $\pm$ 15.4	52.8 $\pm$ 10.7
Benchmark datasets	56.6 $\pm$ 14.3	51.7 $\pm$ 21.8

Data availability

The datasets analysed during this study are available online: CMIP5 model output (https://esgf-node.llnl.gov/search/cmip5/, ESGF, 2022a), CMIP6 model output (https://esgf-node.llnl.gov/search/cmip6/, ESGF, 2022b), the Harmonized World Soils Database (HWSD) and Northern Circumpolar Soil Carbon Database (NCSCD) (https://github.com/rebeccamayvarney/CMIP_soilcarbon_evaluation, ; ), the World Inventory of Soil property Estimates (WISE30sec) (https://www.isric.org/explore/wise-databases, ), Sanderman et al. (2017) soil carbon estimates (1 and 2 m) (10.7910/DVN/QQQM8V, HARVARD, 2022), the Global Soil Dataset for use in Earth system models (GSDE) (http://globalchange.bnu.edu.cn/research/soilw, Land-Atmosphere Interaction Research Group at Sun Yat-sen University, 2022), Global Gridded Surfaces of Selected Soil Characteristics (IGBP-DIS) (https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=569, ORNL-DAAC, 2022), MODIS net primary production (https://modis.gsfc.nasa.gov/data/dataprod/mod17.php, NASA, 2022), CARDAMOM heterotrophic respiration (https://datashare.ed.ac.uk/handle/10283/875, CARDAMOM, 2022), the Copernicus Climate Change Service (C3S) soil moisture gridded dataset (10.24381/cds.d7782f18, ), and the WFDEI meteorological forcing data (https://rda.ucar.edu/datasets/ds314.2/, NCAR, 2022).

The supplement related to this article is available online at: https://doi.org/10.5194/bg-19-4671-2022-supplement.

Author contributions

RMV, SEC, and PMC outlined the evaluation and drafted the manuscript, and RMV completed the analysis and produced the figures. EJB provided the empirical datasets and gave helpful advice for addressing the empirical uncertainty. All the co-authors provided guidance on the study at various times and suggested edits to the draft manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We thank the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modelling groups for producing and making their model output available to enable studies such as this evaluation study. We also thank the providers of the empirical datasets, which enabled us to complete the evaluation in this study. The Taylor diagrams presented in this study were produced using the code available at https://gist.github.com/ycopin/3342888 (last access: 14 January 2022).

Financial support

This research has been supported by the European Research Council, Emergent Constraints on Climate–Land feedbacks in the Earth System project (ECCLES; grant no. 742472), and the European Research Council, Climate–Carbon Interactions in the Current Century (4C) project (grant no. 821003) (Rebecca M. Varney and Peter M. Cox). Sarah E. Chadburn was supported by a Natural Environment Research Council independent research fellowship (grant no. NE/R015791/1). Eleanor J. Burke was supported by the Joint UK BEIS/Defra Met Office Hadley Centre Climate Programme (grant no. GA01101).

Review statement

This paper was edited by Alexey V. Eliseev and reviewed by two anonymous referees.