Elemental stoichiometry and Rock-Eval® thermal stability of organic matter in French topsoils

The quality and quantity of soil organic matter (SOM) are key elements that impact soil health and climate regulation by soils. The Rock-Eval^® thermal analysis technique is becoming more commonly used, as it represents a powerful method for SOM characterization by providing insights into bulk SOM chemistry and thermal stability. In this study, we applied this technique on a large soil sample set from the first campaign (2000–2009) of the French Soil Quality Monitoring Network (RMQS – Réseau de mesures de la qualité des sols). Based on our analyses of ca. 2000 composite surface (0–30 cm) samples collected across mainland France, we observed a significant impact of land cover on both the SOM thermal stability and elemental stoichiometry. Cropland soils had a lower mean hydrogen index value (a proxy for the SOM $\mathrm{H}/\mathrm{C}$ ratio) and a higher thermal stability than grasslands and forests. Regarding the oxygen index (a proxy for the SOM $\mathrm{O}/\mathrm{C}$ ratio), we observed significant differences among the values for croplands, grasslands, and forests. Positive correlations of the temperature parameters with the clay content and pH highlight the protective effect of clay on organic matter as well as the impact of pH on microorganisms' mineralization activity. Surprisingly, we found weak effects of climatic parameters on the thermal stability and stoichiometry of SOM. Our data suggest that topsoil SOM is on average more oxidized and biogeochemically stable in croplands. More generally, the high number and even distribution of data across the whole French territory allow one to build a national interpretative reference for these indicators in surface soils.

1 Introduction

The fate of soil organic carbon (SOC) is crucial from both soil health and climatic perspectives. In terms of soil health, SOC plays an important functional role. Its decomposition by microorganisms provides energy to the whole soil food web as well as key nutrients to plants and soil fauna. SOC also regulates the water cycle by controlling soil structure (Rawls et al., 2003). From a climatic perspective, soils can act as a source or a sink of carbon (Amundson, 2001; Eglin et al., 2010). Maintaining or increasing SOC stocks has become a key policy issue for the coming decades (Rumpel et al., 2018) that raises a number of important scientific challenges regarding our knowledge of SOC dynamics (Dignac et al., 2017).

The evolution of SOC stocks depends on the balance between soil carbon inputs (mostly by plants) and outputs (mostly by microbial decomposition). The persistence of SOC determines soil carbon outputs; thus, estimating the biogeochemical stability of SOC to microbial decomposition (i.e. the difficulty for microorganisms to mineralize SOC) is of paramount importance to infer SOC dynamics (Schmidt et al., 2011; Lehmann and Kleber, 2015). Indeed, a better knowledge of SOC persistence would allow us to refine our estimates of the soil carbon inputs needed to maintain or enhance SOC stocks. However, estimating the biogeochemical stability of SOC is a challenging task because its turnover time encompasses a broad spectrum (ranging from days or weeks to centuries; Balesdent and Guillet, 1982), resulting from a series of interacting SOC stabilization mechanisms. Indeed, SOC can be protected from microbial decomposition due to its chemical nature (e.g. pyrogenic SOC), its interactions with soil mineral surfaces, or its spatial inaccessibility for microbes (Baldock and Skjemstad, 2000; Von Lützow et al., 2006).

Several routine techniques have been proposed to separate fractions that are labile, intermediate, or stable at various timescales (von Lützow et al., 2007; Bispo et al., 2017; Chenu et al., 2015). However, none of these techniques manage to precisely isolate homogeneous fractions with the same biogeochemical stability (von Lützow et al., 2007; Poeplau et al., 2018; Cécillon et al., 2021). Common methods include biological respiration measurements performed during laboratory incubations of soils (e.g. Collins et al., 2000) and various physical (particle size or density) and/or chemical (aqueous or organic extraction) SOC fractionation methods (von Lützow et al., 2007).

Thermal analysis methods have been used for several decades to study the characteristics of soil organic matter (SOM). Many different methods, such as thermogravimetry, differential scanning calorimetry, and evolved gas analysis, exist that measure different variables (Plante et al., 2009). A multitude of variations in temperature ramps, compositions of reaction atmosphere, and measured parameters are encountered within each class of methods. Some thermal analysis methods provide indicators that are related to SOM biogeochemical stability: the more biogeochemically stable SOM is, the more thermally stable and depleted in energy and hydrogen it is (Barré et al., 2016; Sanderman and Grandy, 2020). Among thermal analysis methods, the use of Rock-Eval^® thermal analysis is becoming increasingly common to derive thermal indicators related to SOC biogeochemical stability (Gregorich et al., 2015; Saenger et al., 2015; Cécillon et al., 2018, 2021; Poeplau et al., 2019; Chassé et al., 2021).

The Rock-Eval^® method was developed in the 1970s. Initially intended for the characterization of petroleum source rocks and sediments in order to estimate their potential for hydrocarbon extraction (Espitalié et al., 1977), this method was then adapted to the study of SOM (Disnar et al., 2003). This technique allows for the measurement of the organic and inorganic carbon content of a soil sample as well as the measurement of numerous indicators of the thermal stability and elemental stoichiometry of SOM. Espitalié et al. (1977) showed that the Rock-Eval^® hydrogen index (HI) and oxygen index (OIre6) are a good proxies for the respective $\mathrm{H}/\mathrm{C}$ and $\mathrm{O}/\mathrm{C}$ ratios of organic matter in kerogens and, later, soils. Many temperature parameters can also be calculated to provide insights into the thermal stability of SOM (Gregorich et al., 2015; Sebag et al., 2016; Cécillon et al., 2018). With a rate of one sample per hour at a reasonable price (below USD 50 per sample), Rock-Eval^® thermal analysis is a particularly fast and simple tool to use; therefore, it is well suited to the study of large soil sample sets, such as those collected in the context of national or continental soil monitoring networks.

However, to date, the different existing global soil monitoring networks have not used thermal analysis methods to infer SOC biogeochemical stability. Some of them have focused on SOC physical fractionation schemes, in combination with infrared spectroscopy or environmental variables (e.g. Vos et al., 2018; Viscarra-Rossel et al., 2019; Lugato et al., 2021; Sanderman et al., 2021). Here, we used Rock-Eval^® thermal analysis to investigate the thermal stability and elemental stoichiometry of topsoil samples of the first campaign of the French Soil Quality Monitoring Network (RMQS – Réseau de mesures de la qualité des sols; GisSol; https://www.gissol.fr/le-gis/programmes/rmqs-34, last access: 5 June 2022). The RMQS network has been designed for the long-term monitoring of the soil quality of the whole French territory; the network collects information and samples soils (every 15 years on average) at a set of 2170 sites located on a regular, square grid, forming a systematic sample (Jolivet et al., 2006; English version to be available online). The first campaign took place between 2000 and 2009 in mainland France, covering seven major land cover types (croplands, grasslands, forests, vineyards and orchards, wastelands, environments with little human disturbance, and gardens).

In this study, we first aimed to verify that the Rock-Eval^® method was suitable to characterize SOM on archived soil samples at the scale of a monitoring network. For this purpose, we checked if the organic and inorganic carbon yields of the Rock-Eval^® thermal analysis for soil samples, calculated by comparing Rock-Eval^® estimates to reference methods, were acceptable. Second, we computed several common Rock-Eval^®-based indicators in order to perform an unprecedented country-wide evaluation of the thermal stability and elemental composition of the SOM. Third, thanks to the numerous environmental data available at each RMQS site, we aimed to study the relationships between land cover, climate, and soil properties and the SOC-related indicators derived from Rock-Eval^® thermal analysis.

2 Material and methods

2.1 Topsoil sampling and processing

A full description of the RMQS network and the soil sampling process of its first sampling campaign is available in Jolivet et al. (2006). Briefly, soil is monitored at locations across the French territory on a regular, square grid with a resolution of 16 km. Sampling sites were selected, when possible, at the centre of the cell; otherwise, an alternative site was selected within a 1 km radius of the centre of the cell. This resulted in a total of 2170 RMQS sites in mainland France. At each selected site, 25 topsoil samples (0–30 cm or tilled layer depths) were taken with a spiral soil auger from a 20 m × 20 m sampling area and then mixed to provide a composite sample. Subsoil samples were also taken, but they were not considered in the present study.

The composite samples (5 to 10 kg of bulk soil) were air-dried at 30 ^∘C in trays for 8 to 10 d on average. The samples were then quartered according to NF ISO 11464 to obtain a subsample of ca. 650 g. They were then crushed by hand to break aggregates while preserving calcareous and/or ferro-manganic nodules and sieved at 2 mm. The remains of the composite samples were stored in water-tight plastic buckets. An aliquot of each air-dried and sieved composite sample was then finely ground using a Cyclotec 1093 (FOSS).

Of the 2170 archived aliquots of finely ground topsoil samples from the first RMQS sampling campaign in mainland France, 2037 were recovered and used for this study. When necessary, the samples were manually ground again using an agate mortar to reach the particle size requirements for Rock-Eval^® thermal analysis of soils (below ca. 250 µm).

2.2 Physical and chemical soil analyses

Physical and chemical soil analyses were carried out on the 2 mm sieved composite samples at the Laboratoire d'Analyse des Sols (INRAE, Arras, France). Among the large set of soil properties measured, we selected the following in this study (Jolivet et al., 2006): particle size measurements without decarbonation, in grams per kilogram of sample (Robinson pipette and underwater sieving, method validated in relation to standard NF X31-107), leading to five fractions (clay: ≤ 2 µm; fine silt: 2–20 µm; coarse silt: 20–50 µm; fine sand: 50–200 µm; and coarse sand: 200–2000 µm); pH in a suspension of soil diluted with water (dilution 1 : 5, NF ISO 10390); total carbonate content, in grams per kilogram of sample (volumetric method, NF EN ISO 10693), to estimate the total inorganic carbon, calculated as Cinorg = Total carbonate × 0.12; total carbon content, in grams per kilogram of sample, determined by elemental analysis using dry combustion on non-decarbonated soil; organic carbon content derived from the elemental analysis (TOCea), in grams per kilogram of sample, calculated as Total carbon – Cinorg (NF ISO 10694); total nitrogen, in grams per kilogram of sample (dry combustion, NF ISO 13878); cation exchange capacity (CEC), in centimoles per kilogram of sample (cobaltihexammine chloride extraction, NF X31-130); and free iron oxides, in grams per 100 g, measured with the Tamm method in the dark and Mehra–Jackson method (INRA standard/NF ISO 22036).

2.3 Rock-Eval^® thermal analysis

2.3.1 Thermal analysis process

Rock-Eval^® thermal analyses on the 2037 recovered samples were carried out at the ISTeP – UMR 7193 (Sorbonne Université, Paris, France) according to the routine classically used for soil samples (Disnar et al., 2003; Baudin et al., 2015). Approximately 60 mg of each finely ground topsoil sample was used for the Rock-Eval^® thermal analysis on a RE6 turbo device (Vinci Technologies, 2021). For each analysis, the sample was placed in a special high-temperature-resistant stainless-steel pod, allowing the transport gas to pass through. It first underwent a pyrolysis step under an inert N₂ atmosphere. After a 3 min isotherm at 200 ^∘C, the sample was heated to 650 ^∘C following a temperature ramp of 30 ^∘C min⁻¹. The flame ionization detector (FID) monitored the gaseous emissions of carbon from hydrocarbon compounds (HC_PYR Rock-Eval^® thermogram), while CO (CO_PYR Rock-Eval^® thermogram) and CO₂ (CO2_PYR Rock-Eval^® thermogram) were detected by an infrared detector. The second step is an oxidation (laboratory air atmosphere with CO₂ and H₂O previously removed, i.e. in the presence of oxygen): the sample experienced a 1 min isotherm at 300 ^∘C, was raised to 850 ^∘C following a 20 ^∘C min⁻¹ ramp, and finally remained at a 5 min isotherm at 850 ^∘C. The evolution of CO and CO₂ was again monitored using the infrared detector during the oxidation phase (CO_OX and CO2_OX Rock-Eval^® thermograms). The five resulting thermograms were processed using the Geoworks software (Geoworks V1.6R2, Vinci Technologies, 2021), except for the R-index and I-index parameters (defined in Sect. 2.3.2) which were computed using homemade Python scripts according to the formula proposed by Sebag et al. (2016).

Our Rock-Eval^® thermal analyses campaign included duplicate soil analyses (one every eight samples), which were performed in order to check the reproducibility of the analyses, along with standard analyses (one every nine samples) to check the calibration of the device and identify possible drift in the analysis. The Rock-Eval^® thermal analysis of a soil sample measures its total organic carbon (TOCre6) and total inorganic carbon (MinC) contents that sum to the total carbon content (see Behar et al., 2001, for a detailed description). The organic carbon yield of Rock-Eval^® thermal analysis was defined as TOCre6 $/$ TOCea, its inorganic carbon yield was defined as MinC $/$ Cinorg, and its total carbon yield was defined as (TOCre6 + MinC) $/$ (TOCea + Cinorg). We used the organic carbon yield of Rock-Eval^® thermal analysis to select soil samples among duplicates: only the one with the best yield was conserved. When assessing SOM thermal stability and elemental stoichiometry, it is essential to ensure that SOM analysed by the thermal analysis method corresponds to SOM measured using the reference elemental analysis method. Therefore, we proposed a quality criterion for Rock-Eval^® thermal analysis based on its organic carbon yield, with an arbitrary acceptable range of yields from 0.7 to 1.3.

2.3.2 Rock-Eval^® parameters

Many usual Rock-Eval^® parameters were calculated from the thermograms (Table A1 in the Appendix). First, there are parameters related to carbon quantities: the total organic carbon (TOCre6; grams per kilogram of sample); the total inorganic carbon (MinC; grams per kilogram of sample); the amount of pyrolyzable organic carbon (PC; grams per kilogram of sample); the ratio of pyrolyzable organic carbon over total organic carbon (PCM / TOCre6; no unit); the carbon released during the first pyrolysis isotherm (PseudoS1; grams per kilogram of sample); the carbon released as hydrocarbons during pyrolysis except during the first isotherm (S2; grams per kilogram of sample); and the ratio of carbon released as hydrocarbons during pyrolysis except during the first isotherm over the pyrolyzable organic carbon (S2 / PC; no unit). Second, there are temperature parameters related to the SOC thermal stability. Their calculation was performed over different intervals of integration depending on the thermogram. The upper limits of the integration ranges were selected to exclude CO and CO₂ signals derived from carbonates. The temperature parameters T50_HC_PYR, T70_HC_PYR, and T90_HC_PYR (^∘C) are defined as the respective temperatures at which 50 %, 70 %, and 90 % of the hydrocarbon effluents have been emitted during the pyrolysis ramp (the initial isotherm is excluded; the integration ends at 650 ^∘C). Similarly, T30_CO2_PYR, T50_CO2_PYR, T70_CO2_PYR, and T90_CO2_PYR (^∘C) are the temperatures at which 30 %, 50 %, 70 %, and 90 % of the CO₂ has been emitted during the pyrolysis ramp (the beginning isotherm is excluded; the integration ends at 560 ^∘C); T50_CO_PYR (^∘C) is the temperature at which 50 % of the CO has been emitted during the pyrolysis ramp (the beginning isotherm is excluded; the integration ends at 560 ^∘C). T50_CO2_OX, T70_CO2_OX, and T90_CO2_OX (^∘C) are the respective temperatures at which 50 %, 70 %, and 90 % of the CO₂ has been emitted during the oxidation phase (the integration ends at 611 ^∘C); T50_CO_OX and T70_CO_OX (^∘C) are the respective temperatures at which 50 % and 70 % of the CO has been emitted during the oxidation phase (the integration ends at 850 ^∘C). We also calculated two other parameters previously used in assessing the thermal stability of SOC: the I index, which is related to the thermolabile organic carbon released as hydrocarbon effluents (no unit; Sebag et al., 2016) and the R index, which is the proportion of thermostable organic carbon released as hydrocarbon effluents after 400 ^∘C (no unit; Sebag et al., 2016). Finally, we calculated three Rock-Eval^® parameters, related to the SOM stoichiometry, as outlined in the following. The hydrogen index (HI) is the ratio of emitted hydrocarbons to TOCre6 (in grams of hydrocarbons per kilogram of TOCre6); it is calculated as follows:

$$\begin{array}{}\text{(1)}& \mathrm{HI}={\displaystyle \frac{\mathrm{S}\mathrm{2}\times \mathrm{100}}{\mathrm{TOCre}\mathrm{6}}},\end{array}$$

where S2 is the hydrocarbon signal during pyrolysis (Behar et al., 2001). The oxygen index (OIre6) is the ratio of organic oxygen to TOCre6 (in grams of O₂ per kilogram of TOCre6); it is calculated as follows:

$$\begin{array}{}\text{(2)}& \mathrm{OIre}\mathrm{6}={\displaystyle \frac{\mathrm{16}}{\mathrm{28}}}\times {\displaystyle \frac{\mathrm{S}\mathrm{3}\mathrm{CO}\times \mathrm{100}}{\mathrm{TOCre}\mathrm{6}}}+{\displaystyle \frac{\mathrm{32}}{\mathrm{44}}}\times {\displaystyle \frac{\mathrm{S}\mathrm{3}\times \mathrm{100}}{\mathrm{TOCre}\mathrm{6}}},\end{array}$$

where S3 and S3CO are the organic CO₂ and organic CO signals during pyrolysis, respectively (Behar et al., 2001; Cécillon et al., 2018). The ratio of the hydrogen amount to the oxygen amount is HI $/$ OIre6 (no unit).

As presented above, the treatment of the five thermograms can result in the production of a multitude of Rock-Eval^® parameters. We have decided to present the results of the following parameters in more detail: T50_HC_PYR, T90_HC_PYR, T50_CO2_PYR, T50_CO2_OX, the I index, the R index, HI, and OIre6. The results obtained for some other Rock-Eval^® parameters are presented in Table A2 in the Appendix. The temperature parameters T90_HC_PYR, T50_CO2_PYR, and T50_CO2_OX were selected for two reasons: (1) they are derived from the three different thermograms contributing the most to the Rock-Eval^® signals; (2) they are well correlated with the proportion of centennially stable SOC in temperate soils (Cécillon et al., 2021) and have been used in some previous studies (e.g. Barré et al., 2016; Poeplau et al., 2019). The T50_HC_PYR, I-index, and R-index parameters were selected because they have been used in several previous studies (e.g. Gregorich et al., 2015; Sebag et al., 2016; Matteodo et al., 2018; Soucémarianadin et al., 2018). HI and OIre6 were selected because they are usual Rock-Eval^® parameters and provide insights into the elemental stoichiometry of SOM.

2.4 Climate data

Climate data were extracted from the French SAFRAN database (https://publitheque.meteo.fr/okapi/accueil/okapiWebPubli/index.jsp, last access: 31 March 2022). The daily data were averaged over the 1969–1999 period in order to compute the mean annual temperature (MAT) and mean annual precipitation (MAP) for each site.

2.5 Statistical analysis

We calculated linear regressions without an intercept using the measurements of the organic, inorganic, and total carbon yield of Rock-Eval^® thermal analysis to verify the ability of the Rock-Eval^® thermal analysis to accurately measure the carbon amount of the samples. We chose to use no intercepts because the analysis of several empty pods only showed a very weak signal (TOCre6 < 0.2 g C kg⁻¹).

All of the samples collected from the systematic sampling grid, regardless of their land cover, were analysed using Rock-Eval thermal analysis. This included 847 croplands, 571 forests, 496 grasslands, 57 vineyards, 16 wastelands, 46 sites with little human disturbance, and 4 gardens. Considering the very small number of samples for wastelands and gardens compared with the whole set, we decided not to include them in the following statistical treatments regarding land covers. The number of samples comprising environments with little human disturbance can be considered sufficient for statistical treatment; however, these samples represent a very heterogeneous set (10 miscellaneous subclasses, such as peatlands, alpine grasslands, water edge vegetation, heath, and dry siliceous meadows). Thus, we did not consider it relevant to analyse them as a whole.

To assess the effect of land cover on the Rock-Eval^® parameters, we performed pairwise comparisons of medians using non-parametric Kruskal–Wallis tests (p<0.05) followed by Wilcoxon tests, with p<0.05 for each pair. The correction of p values within the framework of the multiple comparisons was done using the Holm–Bonferroni method. Correlations between parameters were calculated using the Spearman method. We conducted a principal component analysis (with the FactoMineR R library; Lê and Husson, 2008) using all of the observations and 11 pedoclimatic parameters: clay, total silt and total sand contents, pH in water, residual water content, carbonate content, mean annual temperature, mean annual precipitation, Tamm and Mehra–Jackson iron oxyhydroxide contents, and the $\mathrm{C}/\mathrm{N}$ ratio (Fig. A2 in the Appendix). The data processing and statistical analysis were carried out using R software (V4.1.2; R Core Team 2021): the base, datasets, graphics, grDevices, methods, stats, and utils packages were already integrated to R, and the corrplot (Wei and Simko, 2021), car (Fox and Weisberg, 2019), ggplot2 (Wickham, 2016), ggpubr (Kassambara, 2020), factoextra (Kassambara and Mundt, 2020), plot3D (Soetaert, 2021), rstatix (Kassambara, 2021), sf (Pebesma, 2018), and tmap (Tennekes, 2018) packages were added. The point maps of the Rock-Eval^® hydrogen index and T50_CO2_PYR values were obtained using the tmap and sf R packages.

3 Results

3.1 Carbon yields of Rock-Eval^® thermal analysis

Figure 1a presents TOCre6 plotted against TOCea. We observed a high correlation (R² = 0.96, n=2037), despite a few points far from the regression line, and an average carbon yield, corresponding to the slope of the regression, equal to 86 %. Limiting the Rock-Eval^® dataset to samples passing our quality check with respect to the Rock-Eval^® organic carbon yields (yields ranging from 0.7 and 1.3) led to the omission of 145 samples. Another sample was left out because of its TOC content: with a TOC value of 0.57 g kg⁻¹, this sample contained too little organic carbon for the data from the Rock-Eval^® thermal analysis to be routinely exploitable (Khedim et al., 2021). The remaining sample set consisted of samples from 785 croplands, 526 forests, 481 grasslands, 42 vineyards, 14 wastelands, 40 sites with little human disturbance, and 3 gardens. A principal components analysis (PCA) conducted on all of the topsoil samples showed no cluster for the samples with poor organic carbon yields (“rejected”) compared with the samples with good yields (“accepted”) (Fig. A2 in the Appendix). However, there was a significant difference between the medians of the two groups for many pedoclimatic parameters. In particular, the total sand content was on average 76 % higher in the rejected samples compared with accepted samples (101 % higher for the coarse sand and 35 % higher for the fine sand), and the carbonate content was also 67 % higher in the rejected samples compared with accepted samples.

The remaining sample selection logically showed a better agreement between TOCre6 and TOCea, although with lower TOCre6 values on average compared with TOCea (Rock-Eval^® organic carbon yield of 0.87, R² = 0.99, n=1891; Fig. 1b). The inorganic carbon content for this sample selection was slightly overestimated by the Rock-Eval^® thermal analysis (Rock-Eval^® inorganic carbon yield of 1.07, R² = 0.98, n=1891). Finally, the total carbon content measured with Rock-Eval^® thermal analysis for this sample selection is consistent with the total carbon measured using the elemental analysis (Rock-Eval^® total carbon yield of 0.96, R² = 0.99, n=1891).

Figure 1 Carbon yields of Rock-Eval^® thermal analysis. Panel (a) presents the organic carbon yield of TOCre6 as a function of TOCea for all analysed RMQS topsoil (0–30 cm) samples. Panel (b) shows the organic carbon yield of TOCre6 as a function of TOCea for the RMQS topsoil samples limited to those with Rock-Eval^® organic carbon yields ranging from 0.7 to 1.3. Panel (c) presents the inorganic carbon yield of MinC as a function of Cinorg for the RMQS topsoil samples limited to those with Rock-Eval^® organic carbon yields ranging from 0.7 to 1.3. Panel (d) shows the total carbon yield of TOCre6 + MinC as a function of TOCea + Cinorg for the RMQS topsoil samples limited to those with Rock-Eval^® organic carbon yields ranging from 0.7 to 1.3.

3.2 Soil organic matter thermal stability in French topsoils and its relationships with land cover

The summary statistics of many different Rock-Eval^® temperature parameters for the 1891 RMQS topsoil samples with satisfactory Rock-Eval^® organic carbon yields are compiled in Table A2 in the Appendix.

Figure 2 shows the box plots for the six selected parameters (T50_HC_PYR, T90_HC_PYR, T50_CO2_PYR, T50_CO2_OX, the I index, and the R index), focusing on the four major land cover types.

Figure 2 Effect of land cover on topsoil organic carbon thermal stability for the RMQS topsoil (0–30 cm) samples under the four major land covers in France: croplands, forests, grasslands, and vineyards and orchards. The individual panels show the (a) T90_HC_PYR distribution, (b) T50_CO2_PYR distribution, (c) T50_CO2_OX distribution, (d) T50_HC_PYR distribution, (e) I-index distribution, and (f) R-index distribution. Samples are limited to those with Rock-Eval^® organic carbon yields ranging from 0.7 to 1.3. For the box plots, the solid black midline in each box is the median, the lower and upper edges are the respective first (Q₁) and third (Q₃) quartiles, and the lower and upper whiskers are the maximum between the minimum value or the first quartile minus 1.5 times the interquartile range (max [min; $Q_{\mathrm{1}-\mathrm{1.5}}\times$ (Q₃–Q₁)]) and the minimum between the maximum or the third quartile plus 1.5 times the interquartile range (min [max; $Q_{\mathrm{3}+\mathrm{1.5}}\times$ (Q₃–Q₁)]), respectively. Different letters indicate significant differences in the distribution of the values for the land uses according to a Kruskal–Wallis test (p<0.05) and a pairwise Wilcoxon rank sum test (p<0.05). The box width is proportional to the square root of n.

We observed similar results for the temperature parameters T90_HC_PYR, T50_CO2_PYR, and T50_CO2_OX: thermal stability was significantly higher in croplands and in vineyards and orchards compared with forests and grasslands. Topsoil organic carbon was slightly but significantly less thermally stable in forests than in grasslands (Fig. 2a, b, c). Notably, three other Rock-Eval^® parameters related to SOC thermal stability in the HC_PYR thermogram (T50_HC_PYR, the I index, and the R index) showed a different response to land cover (Fig. 2d, e, f). T50_HC_PYR and the R index indicated no significant difference in thermal stability in forests and croplands. The I index indicated a value that was significantly lower in forests than in croplands.

3.3 Elemental stoichiometry of soil organic matter in French topsoils and its relationships with land cover

The summary statistics of different elemental stoichiometry parameters for the 1891 RMQS topsoil samples with satisfactory Rock-Eval^® organic carbon yields are compiled in Table A2 in the Appendix. The HI, OIre6, and $\mathrm{C}/\mathrm{N}$ mean values are 214 g of hydrocarbons per kilogram of TOCre6, 177 g of O₂ per kilogram of TOCre6, and 12.05, respectively.

We observed significantly higher average HI values in both grasslands and forests compared with croplands and with vineyards and orchards (Figs. 3, A1b). In contrast, grasslands and forests showed smaller OIre6 values compared with croplands and with vineyards and orchards (Figs. 3, A1c).

In addition, Fig. 3 highlights that the distribution of the $\mathrm{C}/\mathrm{N}$ ratio on the Rock-Eval^® pseudo van Krevelen diagrams (HI =f(OIre6)) depends on land cover. We observed a slight trend in the $\mathrm{C}/\mathrm{N}$ ratio with the hydrogen and oxygen indices: the $\mathrm{C}/\mathrm{N}$ ratio was higher for high HI and low OIre6. This trend was more pronounced for croplands and forests.

Figure 3 Rock-Eval^® pseudo van Krevelen diagrams (hydrogen index = f(oxygen index)) for the RMQS topsoil (0–30 cm) samples in (a) croplands, (b) grasslands, (c) forests, and (d) vineyards and orchards. Colours indicate the values of the $\mathrm{C}/\mathrm{N}$ ratio. Samples are limited to those with Rock-Eval^® organic carbon yields ranging from 0.7 to 1.3.

3.4 Correlations between Rock-Eval^® indicators of SOM thermal stability and elemental stoichiometry and pedoclimate

Table 1 presents the Spearman correlation coefficient values of the Rock-Eval^® temperature and stoichiometric parameters with the selected pedoclimatic variables. The three selected temperature parameters (T90_HC_PYR, T50_CO2_PYR, and T50_CO2_OX) correlated significantly and positively with the clay content and negatively with the sand content. T90_HC_PYR and T50_CO2_PYR also correlated positively with silt content, although with smaller correlation coefficient values. They strongly and positively (correlation coefficient >0.3) correlated with the water pH, the carbonate content, and the cation exchange capacity, whereas the relationships with the iron oxyhydroxide content were much lower. The three selected temperature parameters were all significantly positively correlated with the mean annual temperature (MAT), and they were negatively correlated with the mean annual precipitation (MAP), although the correlations were weak.

Regarding the indicators of SOM stoichiometry, HI and $\mathrm{C}/\mathrm{N}$ correlated negatively with the clay and silt contents, whereas OIre6 correlated positively with these parameters. HI and $\mathrm{C}/\mathrm{N}$ also correlated negatively with the pH, the cation exchange capacity, and (to a lesser extent) with the carbonate content. They showed a slight negative correlation with the iron oxyhydroxide content measured using the Mehra–Jackson method. As for the thermal parameters, correlations with the climatic variables were smaller on average.

Additionally, the correlation coefficients of TOCre6 with HI, OIre6, $\mathrm{C}/\mathrm{N}$, T90_HC_PYR, T50_CO2_PYR, and T50_CO2_OX were 0.35, −0.34, 0.37, −0.26, −0.21, and −0.05, respectively.

Table 1 Spearman correlation coefficients of the Rock-Eval^® temperature and stoichiometric parameters with the following pedoclimatic variables for the RMQS topsoil (0–30 cm) samples: TOCre6, particle size distribution, pH in water, carbonate content, cation exchange capacity (CEC), iron oxyhydroxides, mean annual temperature (MAT) averaged over 1969–1999, and mean annual precipitation (MAP) averaged over 1969–1999. The analysis was limited to samples with Rock-Eval^® organic carbon yields ranging from 0.7 to 1.3. Absolute values ≥ 0.3 are in bold. The asterisks and superscript letters indicate the p value: ${}^{***}$ between 0 and 0.001, ${}^{**}$ between 0.001 and 0.01, ^* between 0.01 and 0.05, ^a between 0.05 and 0.1, and ^b >0.1.

3.5 Distribution of some Rock-Eval^® indicators of SOM thermal stability and elemental stoichiometry over the French mainland territory

Figure 4 shows the point maps of the HI and T50_CO2_PYR values over the French mainland territory. The missing topsoil samples (133 not included in the initial sample set and 146 rejected due to poor C yields) are distributed over the whole territory with some clusters in the north of the French Alps, the north-east, Corsica, the south-east, and in Landes. The first three clusters come from the 133 samples not included in the initial set. The Landes and south-east clusters are from both the absent samples and the rejected samples: in particular, the soils in Landes contain more sand on average, which is characteristic – as stated above – of the rejected samples. Visually, we noticed an autocorrelation of the values, with HI and T50_CO2_PYR presenting opposite trends on average (the Spearman correlation coefficient between HI and T50_CO2_PYR is −0.69). Mountainous regions (notably the French Alps, the Pyrenees, and the Massif Central) exhibit higher HI values and lower SOC thermal stability. Conversely, plain areas usually presented higher SOC thermal stability and lower HI values, such as in the Paris Basin and in the south-western and south-eastern parts of the country. Brittany, Normandy, and Landes are somewhat exceptions to this rule, as they show high HI values and a relatively low SOC thermal stability. Figure 5 shows the land cover at each sampling site.

Figure 4 Point maps of two Rock-Eval^® parameters – (a) hydrogen index values and (b) T50_CO2_PYR values – on the French mainland territory for the RMQS topsoil (0–30 cm) samples; (c) a map of the land cover at each sampling site (the numbers in square brackets correspond to the number of sites for each land cover in our final dataset; n=1891), modified from Jolivet (2011); and (d) a map of the main regions used for the interpretation.

4 Discussion

4.1 Carbon yields of Rock-Eval^® thermal analysis

Our average, the organic carbon yield (0.86; Fig. 1a) was in line with previous studies. Indeed, Disnar et al. (2003) (0.91) and Cécillon et al. (2018, 2021) (organic carbon yield from 0.90 to 0.96 depending on the sites) obtained slightly higher yields, whereas Saenger et al. (2013) reported lower yields (0.77). However, some samples presented high discrepancies between their TOCea and TOCre6 values. Samples with a TOCre6 value strongly differing from its corresponding TOCea value were systematically reanalysed using Rock-Eval^®, which confirmed their first TOCre6 measurement. The outliers with respect to the organic carbon yield were, thus, not related to a problem in their Rock-Eval^® measurement. These very different values, which concern a few dozen samples, could have different origins, such as error in sample labelling, the division into aliquots, grinding, or storage conditions. Indeed, for the same sample, the powders used for the elemental analysis and the Rock-Eval^® thermal analysis did not come from the same aliquot. In addition, the elemental analyses were performed shortly after sampling, whereas the samples analysed in Rock-Eval^® were stored for about 15 years. Therefore, we can expect slightly better yields when elementary and Rock-Eval^® analysis are performed with less time between both as well as when they are performed on the exact same powders. This is what we plan for the samples of the second RMQS sampling campaign. The very different TOCea and TOCre6 values could also be due, for some samples, to a mismeasurement of the total carbonate content, leading to a miscalculation of the inorganic and organic carbon contents. This hypothesis could be plausible, as the median value of the carbonate content was significantly higher in the rejected samples. The last hypothesis originates from the high content of sand in the rejected samples: sandy samples are more heterogeneous; thus, the material used to determine the TOCea is more likely to differ from that used to determine the TOCre6, compared with when the sand content is lower. Moreover, the physical state of organic matter in sandy soils can be different from other soils. Disnar et al. (2003) encountered “pellets” of SOM in sandy soils, which can strongly influence the TOCea and TOCre6 results.

The samples presenting a high discrepancy between TOCea and TOCre6 were not considered further in the analysis. As stated above, we restricted our study to the samples with an organic carbon yield ranging from 0.7 to 1.3. This subjective threshold is a quality threshold to ensure that the samples analysed using Rock-Eval^® were the same as the samples analysed using elemental analysis, on which all studies conducted on the first campaign of RMQS rely. This selection only marginally improved the average organic carbon yield (0.87; Fig. 1b), and organic carbon was still underestimated by Rock-Eval^®. Conversely, the inorganic carbon yield was slightly overestimated (1.07; Fig. 1c). As a result, the yield of total carbon (organic + inorganic carbon) was close to 1.00 (Fig. 1d). This suggests that almost all sample carbon is detected by the Rock-Eval^® machine in the five thermograms but that a small part of the organic carbon is erroneously attributed to inorganic carbon. This may be due to a slight misplacement of the boundary between organic and inorganic carbon, probably in the S3 and S3CO signals. Also, the S3^′CO signal is attributed half to organic carbon and half to inorganic carbon due to potential Boudouard reactions, which are not always verified (Baudin et al., 2015; see also e.g. Behar et al., 2001, for a definition of the Rock-Eval^® peaks). It is of note that, as MinC and TOCre6 are very well correlated with Cinorg and TOCea (R²>0.98), it should be possible to draw a correction formula to assess TOCea and Cinorg using Rock-Eval^® with high accuracy. This would allow for the simultaneous determination of organic C and inorganic C in less than 1 h with no risk of error due to erroneous decarbonation.

4.2 Thermal stability of soil organic carbon in French topsoils

We have observed that the thermal stability defined according to different Rock-Eval^® parameters varies in French topsoils. We can investigate whether these variations are consistent with our knowledge of SOC biogeochemical stability. SOC biogeochemical stability is on average higher in croplands and vineyards than in forest or grassland soils (Poeplau and Don, 2013). Indeed, fresh organic carbon inputs to soil are usually higher in forest and grassland compared with croplands, where human exportation of biomass is higher (Murty et al., 2002). As a result, SOC fractions with a lower mean residence time in soils and a lower thermal stability can be more abundant in forests and grasslands than in croplands. For instance, several studies have reported that carbon in particulate organic matter (a relatively more labile form of SOC) contributes more to total SOC in forest and grassland compared with croplands (e.g. Guo and Gifford, 2002; Poeplau et al., 2011; Poeplau and Don, 2013; Lugato et al. have 2021). Moreover, agricultural practices may also speed up SOC mineralization, thereby further limiting the accumulation of labile SOC fractions. For instance, Balesdent et al. (1990) observed that the tillage practices lead to a significantly higher mineralization than no tillage. Combining the effects of lower carbon inputs and mineralization-enhancing practices, croplands contain less biogeochemically labile SOC on average than forests and grasslands.

Thermal stability, as assessed using T90_HC_PYR, T50_ CO2_ PYR and T50_CO2_OX, was the highest in vineyards and orchards and in croplands compared with forest and grassland soils (Fig. 2). These results suggest that, overall, SOC thermal stability, as assessed using these Rock-Eval^® parameters, is related to SOC biogeochemical stability. This is in good agreement with previous results obtained on smaller datasets (Barré et al., 2016; Poeplau et al., 2019; Cécillon et al., 2021). On the contrary, there was no consistent relationship between thermal stability and expected biogeochemical stability when the thermal stability was measured using T50_HC_PYR, the R index, and the I index (Fig. 2). For soils with highly contrasted biogeochemical stability, Cécillon et al. (2021) reported that the relationship between thermal stability and biogeochemical stability was weaker for T50_HC_PYR, the R index, and the I index. Our results showed that this relationship even disappears when considering datasets with more heterogeneous topsoil samples. The use of the Rock-Eval^® temperature parameters T90_HC_PYR, T50_CO2_PYR, and T50_CO2_OX should, therefore, be preferred when seeking to measure thermal stability indicators directly related to biogeochemical stability.

T90_HC_PYR, T50_CO2_PYR, and T50_CO2_OX were all strongly and positively correlated with the clay content and negatively correlated with the sand content (Table 1). In a previous study, Soucémarianadin et al. (2018) did not observe any correlation between T50_CO2_OX and clay or sand content; however, their study was conducted on forest soils only and on a greatly reduced number of study sites. Soil clay fractions interact with microbial compounds, resulting in the formation of organo-mineral complexes in which SOC has a high biogeochemical stability (e.g. Lehmann and Kleber, 2015). Therefore, we can hypothesize that clay-rich soils are also richer in biogeochemically stable carbon. The positive correlation between clay content and SOC thermal stability as well as the good correlations between the CEC, which depends on the first order of the clay content, and SOC thermal stability would then be another illustration of the link between SOC thermal and biogeochemical stabilities. Iron oxides are mineral compounds that are also supposed to protect SOC from decomposition. In this respect, the inconsistent (Mehra–Jackson iron) or even negative correlations (Tamm iron) between T90_HC_PYR, T50_CO2_PYR, and T50_CO2_OX and iron oxides were not expected. These weak correlations could be attributed to the fact that the range of iron oxide contents is relatively small in our set of topsoils.

T90_HC_PYR, T50_CO2_PYR, and T50_CO2_OX were all positively correlated with pH. Such a correlation between T50_CO2_OX and pH has already been observed by Soucémarianadin et al. (2018) for a set of French forest soils. Acidity may protect SOM from degradation by microorganisms (Clivot et al., 2021), by reducing their activity, which is actually observed in low-pH bogs. Therefore, we can hypothesize that acidity slows down SOM mineralization which can favour the accumulation of labile SOC components. As these labile SOC fractions would appear thermally unstable, it would explain the positive relationship between pH and Rock-Eval^® indicators of SOC thermal stability.

T90_HC_PYR, T50_CO2_PYR, and T50_CO2_OX showed weak but significant positive correlations with MAT averaged over 1969–1999 (Table 1). Such a correlation has also been observed in Soucémarianadin et al. (2018) for French forest soils. As soil microbial activity and, thus, SOC mineralization increase with temperature (Rey and Jarvis, 2006), we can expect the SOC labile fractions to be more rapidly processed at higher temperature. This would be in line with the observed positive correlations between MAT and the three selected thermal stability indicators. The relatively weak (Spearman ρ value below 0.2) correlations could be due to the fact that MAT also plays a role in carbon inputs to the soil. Indeed, if higher SOC mineralization was balanced by increased biomass inputs, it would dampen the relation between MAT and the SOC biogeochemical stability. In a similar way, the weak negative correlation between MAP and thermal stability may be explained by the complex effect of MAP on SOC biogeochemical stability: increased soil moisture stimulates SOC processing up to a certain point (Moyano et al., 2013) and also influences net primary production and, therefore, soil carbon inputs. In any case, the relationships between SOC and MAP or MAT are hard to disentangle (Chen et al., 2019). Another explanation for the weak values is that the climatic data were obtained on an 8 km × 8 km grid and do not have the same precision that they would have had if a weather station had been deployed at each site. This probably adds noise to the correlation.

The point map representing SOC thermal stability over mainland France (Fig. 4b) illustrates the relationships between SOC thermal stability, land cover, climate, and pedological variables. Mountainous regions (e.g. the Massif Central, Alps, and Pyrenees) dominated by forest and grassland with a low MAT and relatively high SOC contents (the latter according to Martin et al., 2011) had a lower SOC thermal stability. Plains dominated by croplands with intensive agricultural practices and with relatively low SOC contents, such as the Paris Basin, showed high SOC thermal stability. The southern part of France, which has a warmer MAT, dominant vineyard and cropland land cover, and relatively low SOC contents, also presented high SOC thermal stability. The lower SOC thermal stabilities observed in Brittany and Normandy (which are agricultural regions) could be explained by the higher proportion of livestock. Therefore, in addition to the presence of grasslands in these regions, the cultivated soils in Brittany and Normandy are more likely to receive the repeated application of exogenous organic matter.

4.3 Elemental stoichiometry of soil organic matter in French topsoils

Higher HI and lower OIre6 values were observed in forests and grasslands compared with croplands and vineyards. This trend has also been observed in previous studies (Disnar et al., 2003; Saenger et al., 2013; Sebag et al., 2016). This confirms that HI and OIre6 can be good proxies for SOC biogeochemical stability. Indeed, as previously observed, biogeochemically stable SOC is more oxidized and H depleted (Barré et al., 2016; Poeplau et al., 2019; Cécillon et al., 2021).

The pseudo van Krevelen diagrams (Fig. 3) show a high variability in the $\mathrm{C}/\mathrm{N}$ ratio between land cover classes: the $\mathrm{C}/\mathrm{N}$ ratio is higher in forest topsoils than in grasslands, croplands, or vineyards. This is classically explained by the fact that SOC is on average less processed in forests and grasslands compared with croplands and vineyards (Cotrufo et al., 2019) as well as by the higher $\mathrm{C}/\mathrm{N}$ ratio of the biomass inputs to soil in forests. Indeed, the biotransformation of organic matter tends to lower its $\mathrm{C}/\mathrm{N}$ ratio and oxidize it (Cleveland and Liptzin, 2007). This is in good agreement with the observed trends of decreasing HI and increasing OIre6 with decreasing $\mathrm{C}/\mathrm{N}$ (Fig. 3).

SOM elemental stoichiometry presented correlation patterns with land cover, climate, and pedological variables that were similar to those observed for SOM thermal stability. As shown in Table 1, HI and OIre6 are negatively and positively correlated with pH, respectively, as previously observed by Soucémarianadin et al. (2018) in French forest soils. This would be in line with acidity slowing down the mineralization of H-enriched labile SOC fractions (Clivot et al., 2019). The negative correlation between clay content and HI could be explained by the fact that the presence of clays can promote the protection of microbially processed H-depleted SOM. Similar to what was observed for SOM thermal stability, relationships between elemental stoichiometry and climate variables are weak, probably because climate plays a role in both soil carbon inputs and outputs in opposite ways (climate conditions enhancing SOC mineralization usually also enhance fresh SOM inputs).

The point map of HI in mainland France (Fig. 4a) illustrates the effect of land cover, climate, and pedological variables on SOM elemental stoichiometry. Regions dominated by grassland and forest (Fig. 4d), such as mountainous regions, the Landes forest, or the forest-dominated eastern part of France, are characterized by a relatively H-enriched SOM. Conversely, regions with a high MAT and dominant cropland, vineyard, and orchard land covers are characterized by a relatively H-depleted SOM. Both point maps of thermal stability and HI (Fig. 4) also illustrate the relationships previously observed between these Rock-Eval^® parameters (Barré et al., 2016; Cécillon et al., 2021).

5 Conclusions

This study is an unprecedented effort to carry out widespread thermal analysis measurements on a national soil quality monitoring network. It demonstrated that Rock-Eval^® may be used as a rapid and cost-effective method to assess the thermal stability and elemental stoichiometry of SOM on national soil monitoring networks. The very satisfying organic and inorganic carbon yields could make Rock-Eval^® thermal analysis a very suitable tool for research work in carbonate soils or even for routine soil analysis if commercial laboratories take advantage of the method. Our results highlighted the influence of land cover and pedoclimatic variables on SOM thermal stability and elemental stoichiometry. They suggested that some Rock-Eval^® temperature parameters describing SOC thermal stability (T90_HC_PYR, T50_CO2_PYR, and T50_ CO2_OX) could be used as reliable proxies for SOC biogeochemical stability, whereas other parameters (T50_HC_PYR, the R index, and the I index) could not. Our study also opened wide perspectives for future research. In the short term, these Rock-Eval^® results on French topsoils can be used as input to the PARTY_SOC machine learning model (Cécillon et al., 2021) to infer the size of the centennially stable SOC fraction. They can also be compared to other proxies for SOC biogeochemical stability, such as SOM physical fractionation results. In the medium term, it will be interesting to test whether this analytical information can be used to improve the accuracy of SOC stock evolution simulations at the scale of a national soil monitoring network, as was observed for the Andriulo–Mary–Guérif (AMG) model of SOC dynamics in several French long-term agronomic experiments (Kanari et al., 2022).

Appendix A

Figure A1 Effect of land cover on topsoil organic carbon stoichiometry for the RMQS topsoil (0–30 cm) samples under the four major land covers in France: croplands, forests, grasslands, and vineyards and orchards. The individual panels show the (a) $\mathrm{C}/\mathrm{N}$ distribution, (b) HI distribution, and (c) OIre6 distribution. Samples are limited to those with Rock-Eval^® organic carbon yields ranging from 0.7 to 1.3. For the box plots, the solid black midline of each box is the median, the lower and upper edges are the respective first (Q₁) and third (Q₃) quartiles, and the lower and upper whiskers are the maximum between the minimum value or the first quartile minus 1.5 times the interquartile range (max [min; $Q_{\mathrm{1}-\mathrm{1.5}}\times$ (Q₃–Q₁)]) and the minimum between the maximum or the third quartile plus 1.5 times the interquartile range (min [max; $Q_{\mathrm{3}+\mathrm{1.5}}\times$ (Q₃–Q₁)]), respectively. Different letters indicate significant differences in the distribution of the values for the land uses according to a Kruskal–Wallis test (p<0.05) and a pairwise Wilcoxon rank sum test (p<0.05). The box width is proportional to the square root of n.

Figure A2 Score of the 2037 samples on axes 1 and 2 of the principal component analysis on 11 pedoclimatic parameters: clay, total silt and total sand contents, pH in water, water content, carbonate content, mean annual temperature, mean annual precipitation, Tamm and Mehra–Jackson iron oxyhydroxide contents, and the $\mathrm{C}/\mathrm{N}$ ratio for the RMQS topsoil (0–30 cm) samples. Samples with an organic carbon yield between 0.7 and 1.3 are plotted in light blue, whereas samples with an organic carbon yield < 0.7 or > 1.3 are plotted in dark blue.

Table A1 Description of the Rock-Eval^® parameters and their calculation.

Table A2 Minimum, maximum, mean, first quartile, third quartile, median, and standard deviation values of the Rock-Eval^® parameters for the RMQS topsoil (0–30 cm) samples limited to those with Rock-Eval^® organic carbon yields ranging from 0.7 to 1.3 (n=1891, “total”) and by land cover.