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 of soil health and climate regulation by soils. The Rock-Eval® thermal analysis technique is increasingly used as it represents a powerful method for SOM characterization by providing insights on bulk SOM chemistry and thermal stability. In this study, we applied this technique 15 on a large soil sample set from the first campaign (2000–2009) of the French monitoring network of soil quality: RMQS. Based on our analyses on ca. 2000 composite surface (0–30 cm) samples taken all over mainland France, we observed a significant impact of land cover on both SOM thermal stability and elemental stoichiometry. Cropland soils had a lower mean value of hydrogen index (a proxy for SOM H/C ratio) and a higher thermal stability than grasslands and forests. Regarding the oxygen index (a proxy for SOM O/C ratio), we observed significant differences in values for croplands, 20 grasslands and forests. Positive correlations between the temperature parameters on the one hand and the clay content and pH on the other hand highlight the protective effect of clay on organic matter and 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 repartition of data on the whole French territory allow to build a national interpretative 25 referential for these indicators in surface soils.


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 and key nutrients to plants and soil fauna. SOC also regulates the water cycle through controlling soil structure (Rawls et al., 2003). 30 From a climatic perspective, soils can act both as a source or a sink of carbon (Amundson, 2001;Eglin et al., 2010).
3 analysis is a particularly fast and simple tool to use, and is therefore well suited to the study of large soil sample sets, such as the ones collected in the context of national or continental soil monitoring networks.
However, so far, the different existing soil monitoring networks worldwide 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;70 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; GIS Sol; https://www.gissol.fr/le-gis/programmes/rmqs-34). The RMQS network has been designed for the long-term monitoring of the soil quality of the whole French territory by collecting information and sampling soils on, every 15 years in average, a set of 2170 sites at the locations of a regular, square grid thus forming a 75 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 7 major land cover types (croplands, grasslands, forests, vineyards & orchards, wastelands, poorly human-disturbed environments, and gardens).
In this study, we aimed in the first place to verify that the Rock-Eval® method was suited 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 80 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 at studying the relationships between land-cover, climate and soil properties and the SOC-related indicators derived from Rock-Eval® thermal analysis. 85 2 Material and methods 2.1 Topsoil sampling and processing A full description of the RMQS and of the soil sampling process of its first sampling campaign is available in Jolivet et al., 2006. Briefly, the soil is monitored at the locations of a regular, square grid with a resolution of 16 km. A sampling site was settled when possible at the center of the cell; otherwise, an alternative site was taken within a 1 km radius from the center of 90 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 then mixed to provide a composite sample. Subsoil samples were also taken, but 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 days on average. The samples were then quartered according to NF ISO 11464 to obtain a sub-sample of ca. 650 g. They were then crushed by hand to 95 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 100 mortar to reach the particle size requirements of Rock-Eval® thermal analysis of soils (below ca. 250 µm).

Physical and chemical soil analyses
The 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 in this study the following ones (Jolivet et al., 2006): particle-size measurements without decarbonation in g.kg -1 of sample (Robinson pipette 105 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; coarse sand: 200-2000 µm); pH in water (NF ISO 10390, dilution with ⅕); total carbonate content in g.kg -1 of sample (volumetric method NF EN ISO 10693) to estimate the total inorganic carbon Cinorg = Total carbonate ✕ 0.12 in g.kg -1 of sample; total carbon content (g.kg -1 of sample) determined by elemental analysis using dry combustion on non-decarbonated soil; organic carbon content derived from the elemental 110 analysis (TOCea; g.kg -1 of sample) calculated as Total carbon -Cinorg (NF ISO 10694); total nitrogen in g.kg -1 of sample (dry combustion NF ISO 13878); CEC in cmol+.kg -1 of sample (cobaltihexammine chloride extraction NF X31-130); free iron oxides in g/100 g measured with the Tamm method in the dark and Mehra--Jackson method (INRA standard/NF ISO 22036).

Rock-Eval® thermal analysis 115
2.3.1 Thermal analysis process Rock-Eval® thermal analyses were carried out at the UMR 7193 ISTeP (Sorbonne Université, Paris, France) on the 2037 recovered samples 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, France). For each analysis, the sample was placed in a special high-temperature resistant 120 stainless steel pod allowing the transport gas to pass through. It first underwent a pyrolysis step under an inert N2 atmosphere. After a three-minute 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 CO2 (CO2_PYR Rock-Eval® thermogram) were detected by an infrared detector. The second step is an oxidation (laboratory air atmosphere with CO2 and 125 H2O previously removed, i.e., in presence of oxygen): the sample experienced a one-minute isotherm at 300 °C, then was raised to 850 °C following a 20 °C/min ramp, and finally remained on a five-minutes isotherm at 850 °C. The evolution of Eval® thermograms). The five resulting five thermograms were processed using the Geoworks software (Vinci Technologies, Geoworks V1.6R2), except for the parameters R-index and, I-index (defined in 2.3.2) which were computed 130 using homemade Python scripts according to the formula proposed by Sebag et al. (2016).
Our Rock-Eval® thermal analyses campaign included duplicated soil analyses (one every eight samples), 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 a possible drift in the analysis. The Rock-Eval® thermal analysis of a soil sample measures its total organic carbon (TOCre6) and total inorganic carbonC (MinC) contents that sum to total carbon content (see Behar et 135 al., 2001 for a detailed description). The organic carbon yield of Rock-Eval® thermal analysis was defined as TOCre6/TOCea, while 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 analyzed by the thermal analysis method corresponds to SOM measured 140 using the reference elemental analysis method. We therefore 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.

Rock-Eval® parameters
Many usual Rock-Eval® parameters were calculated from the thermograms (Table A). First there are parameters related to carbon quantities : the total organic carbon (TOCre6; g.kg -1 sample); the total inorganic carbon (MinC; g.kg -1 sample); the 145 amount of pyrolyzable organic carbon (PC; g.kg -1 sample); the ratio of pyrolyzable organic carbon over total organic carbon (PC/TOCre6; no unit); the carbon released during the first pyrolysis isotherm (PseudoS1; g.kg -1 sample); the carbon released as hydrocarbons during the pyrolysis except during the first isotherm (S2; g.kg -1 sample); the ratio of carbon released as hydrocarbons during the 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 150 different intervals of integration depending on the thermogram. The upper limits of the integration ranges were selected to exclude CO and CO2 signals derived from carbonates. The temperature parameter T50_HC_PYR (respectively T70_HC_PYR and T90_HC_PYR; °C) is defined as the temperature at which 50 % (respectively 70 % and 90 %) of the hydrocarbon effluents hasve been emitted during the pyrolysis ramp (the initial isotherm is excluded; the integration ends at 650 °C). Similarly, T30_CO2_PYR (respectively T50_CO2_PYR, T70_CO2_PYR and T90_CO2_PYR; °C) is the 155 temperature at which 30 % (respectively 50 %, 70 % and 90 %) of the CO2 hasve 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 hasve been emitted during the pyrolysis ramp (the beginning isotherm is excluded; the integration ends at 560 °C).
T50_CO2_OX (respectively T70_CO2_OX and T90_CO2_OX; °C) is the temperature at which 50 % (respectively 70 % and 90 %) of the CO2 hasve been emitted during the oxidation phase (the integration ends at 611 °C); T50_CO_OX 160 (respectively T70_CO_OX; °C) is the temperature at which 50 % (respectively 70 %) of the CO have 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 (related to the thermolabile organic carbon released as hydrocarbon effluents, Sebag et al., 2016;no unit) and the R-index (the proportion of thermostable organic carbon released as hydrocarbon effluents after 400 °C, Sebag et al., 2016;no unit). Finally, we calculated the three following Rock-Eval® parameters, related to the SOM 165 stoichiometry. The Hydrogen Index, HI, is the ratio of emitted hydrocarbons to TOCre6 (unit g HC.kg -1 TOCre6); it is calculated following Eq. (1): where S2 is the hydrocarbons signal during pyrolysis (Behar et al., 2001). The Oxygen Index, OIre6, is the ratio of organic oxygen to TOCre6 (unit g O2.kg -1 TOCre6); it is calculated following Eq. (2): 170 where S3 and S3CO are respectively the organic CO2 and organic CO signals during pyrolysis (Behar et al., 2001;Cécillon et al., 2018). The ratio of hydrogen amount to 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 on the following parameters in more detail: T50_HC_PYR; 175 T90_HC_PYR; T50_CO2_PYR; T50_CO2_OX; I-index; R-index; HI and OIre6. The results obtained for some other Rock-Eval® parameters are presented as Supplementary Information. The temperature parameters T90_HC_PYR, T50_CO2_PYR and T50_CO2_OX were selected as they are derived from the 3 different thermograms contributing the most to the Rock-Eval® signals, and as they are well correlated to the proportion of centennially stable SOC in temperate soils (Cécillon et al., 2021) and were used in some previous studies (e.g. Barré et al., 2016;Poeplau et al., 2019). The parameters T50_HC_PYR, 180 I-index and R-index were selected as they were 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 as they are usual Rock-Eval® parameters and they give insights on the elemental stoichiometry of SOM.

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

Statistical analysis
We calculated linear regressions without 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 190 amount of the samples. We chose to use no intercepts as the analysis of several empty pods only showed very weak signal (TOCre6 < 0.2 gC.kg -1 ).
All the samples collected from the systematic sampling grid, regardless of their land cover, were analyzed using Rock-Eval thermal analysis. This includes 847 croplands, 571 forests, 496 grasslands, 57 vineyards, 16 wastelands, 46 poorly humandisturbed sites and 4 gardens. Considering the very small number of samples for wastelands and gardens compared to the 195 whole set, we decided not to include them in the following statistical treatments regarding the land covers. The number of poorly human-disturbed samples can be considered sufficient for statistical treatment, however they represent a very heterogeneous set of samples (10 miscellaneous subclasses such as peatlands, alpine grasslands, water edge vegetation, heath, dry siliceous meadows, etc.). We did not consider relevant to analyze them as a whole.
To assess the effect of land cover on the Rock-Eval® parameters, we performed pairwise comparisons of medians by non-200 parametric Kruskal-Wallis tests (P <0.05) followed by Wilcoxon tests, with P < 0.05 for each pair. The correction of pvalues in the framework of the multiple comparisons was done with the Holm--Bonferroni method. Correlations between parameters were calculated using the Spearman method. We conducted a principal component analysis (using the R library FactoMineR) using all 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 205 iron oxyhydroxides contents and C/N ratio (Fig. D). The data processing and statistical analysis were carried out with the R software (V4.1.2; R Core Team 2021. Packages integrated to R: base, datasets, graphics, grDevices, methods, stats, utils. Packages added: corrplot, car, ggplot2, ggpubr, factoextra, plot3D, rstatix, sf, tmap). 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 210 3.1 Carbon yields of Rock-Eval® thermal analysis Figure 1a presents the TOCre6 plotted against TOCea. We observed a high correlation despite a few points far from the regression line (R²=0.96, n = 2037) 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 regarding the Rock-Eval® organic carbon yields (yields ranging from 0.7 and 1.3) left out 145 samples. Another sample was left out because of its TOC content: with a value 215 of 0.57 g. kg -1 , this sample contains 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 consists in samples from 785 croplands, 526 forests, 481 grasslands, 42 vineyards, 14 wastelands, 40 poorly human-disturbed sites and 3 gardens. A principal components analysis (PCA) conducted on all the topsoil samples showed no cluster for the samples with poor organic carbon yields ("rejected") compared to the samples with good yields ("accepted") ( Fig. D). However, there was a significant difference 220 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 to 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 to accepted samples.
The remaining sample selection logically showed a better agreement between TOCre6 and TOCea, yet with on average lower TOCre6 values compared with TOCea (Rock-Eval® organic carbon yield of 0.87, R²=0.99, n = 1891; Fig. 1b). The 225 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).  respectively the maximum between the minimum value or the 1 st quartile minus 1.5 times the interquartile range (max[min;Q1-1.5*(Q3-Q1)]) and the minimum between the maximum or the 3 rd quartile plus 1.5 times the interquartile range (min[max;Q3+1.5*(Q3-Q1)]). Different letters indicate significant differences in the distribution of the values for the land uses with the Kruskal-Wallis test (P < 0.05) and multiple Wilcoxon 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 250 stability was significantly higher in croplands and vineyards & orchards compared to forests and grasslands. Topsoil organic carbon was slightly but significantly less thermally stable in forests than in grasslands (Fig. 2a, 2b, 2c). Notably, three other Rock-Eval® parameters related to SOC thermal stability in the HC_PYR thermogram (T50_HC_PYR, I-index and R-index) showed a different response to land cover (Fig. 2d, 2e, 2f). The T50_HC_PYR and R-index indicated no significant difference in thermal stability in forests and croplands. The I-index indicated a value significantly lower in forests than in 255 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 Supplementary Information (Table B). The HI (respectively OIre6 and C/N) mean value is 214 g HC.kg -1 TOCre6 (respectively 177 g O2.kg -1 TOCre6 and 12.05). 260 We observed significantly higher average values for the HI in both grasslands and forests compared to croplands and vineyards & orchards (Fig. 3, Fig. Cb). In contrast, grasslands and forests showed smaller values of OIre6 compared to croplands and vineyards & orchards (Fig. 3, Fig. Cc).
In addition, Figure 3 highlights that the distribution of the C/N ratio on the Rock-Eval® pseudo van Krevelen diagrams (HI=f(OIre6)) depends on land cover. We observed a slight trend of the C/N ratio with the hydrogen and oxygen indices: the 265 C/N ratio was higher for high HI and low OIre6. This trend was more pronounced for croplands and forests. Samples are limited to those with Rock-Eval® organic carbon yields ranging from 0.7 to 1.3.

Correlations between Rock-Eval® indicators of SOM thermal stability and elemental stoichiometry and 270
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, 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, however with smaller correlation coefficient 275 values. They strongly and positively (correlation coefficient > 0.3) correlated with the water pH, the carbonate content and the cation exchange capacity, while the relationships with iron oxyhydroxides content were much lower. The three selected temperature parameters were all significantly positively correlated to Mean Annual Temperature (MAT) and negatively correlated to Mean Annual Precipitation (MAP), although the correlations were weak. Table 1: Spearman correlation coefficients of the Rock-Eval® temperature and stoichiometric parameters with pedoclimatic 280 variables (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  for the RMQS topsoil (0-30 cm) samples. 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 indicate the p-value: 0 | *** | 0.001 | ** | 0.01 | * | 0.05 | • | 0.1 | X.

285
Regarding the indicators of SOM stoichiometry, HI and C/N correlated negatively with the clay and silt contents, contrary to OIre6 which correlated positively. HI and C/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 oxyhydroxides content measured by the Mehra--Jackson method. As for the thermal parameters, correlations with the climatic variables were on average smaller. 290 Additionally, the correlation coefficient of TOCre6 with HI was 0.35 (respectively -0.34 with OIre6, 0.37 with C/N, -0.26 with T90_HC_PYR, -0.21 with T50_CO2_PYR, -0.05 with T50_CO2_OX).
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 295 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, north-east, Corsica, south-east and in the 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 from the rejected samples: in particular, the soils in the Landes contain on average more sand, which is characteristic -as stated above -of the rejected samples. Visually, we noticed an autocorrelation of the values, HI and 300 T50_CO2_PYR presenting on average opposite trends (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 as in the Paris Basin and in the south-west and south-east part of the country. Brittany, Normandy and the Landes are somewhat exceptions to this rule as they show high HI values and a relatively low SOC thermal stability. 305 Figure 5 shows the land cover at each sampling site.  which concern a few dozens of samples, could have different origins such as error on sample labeling, aliquoting, 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 analyzed in Rock-Eval® were stored for about fifteen years. We can therefore expect slightly better yields when elementary and Rock-Eval® analysis are performed with less time between both, and on the exact same 325 powders. This is what we plan for the samples of the second RMQS sampling campaign. The very different values between TOCea and TOCre6 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 330 likely to differ from the one used to determine the TOCre6, than when 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 results of TOCea and TOCre6.
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 organic carbon yield ranging from 0.7 to 1.3. This subjective 335 threshold is a quality threshold to ensure that the samples analyzed using Rock-Eval® were the same as the samples analyzed using elemental analysis on which rely all studies conducted on the first campaign of the RMQS. This selection only marginally improved the average organic carbon yield (0.87; Fig. 1b) and organic carbon was still underestimated by Rock-Eval®. Conversely, 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 340 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 is not always verified (Baudin et al., 2015; see e.g. Behar et al. (2001) for a definition of the Rock-Eval® peaks). Of note, as MinC and TOCre6 are very well correlated to Cinorg and TOCea (R² > 0.98), it 345 should therefore be possible to draw a correction formula to assess TOCea and Cinorg using Rock-Eval® with high accuracy. This would allow determining simultaneously in less than 1 h organic C and inorganic C with no risk of error due to erroneous decarbonation.

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. 350 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 compared to forest or grassland soils (Poeplau and Don, 2013). Indeed, fresh organic carbon inputs to soil are usually higher in forest and grassland compared to croplands where human exportation of biomass is higher (Murty et al., 2002). As a result, SOC fractions with lower mean residence time in soils and lower thermal stability can be more abundant in forests and grasslands compared to croplands. For instance, 355 several studies reported that carbon in particulate organic matter (a relatively more labile form of SOC) contributes more to total SOC in forest and grassland compared to cropland (e.g. Guo and Gifford, 2002;Poeplau et al., 2011;Poeplau and Don, 2013;Lugato et al. 2021). Moreover, agricultural practices may also speed up SOC mineralization 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-360 enhancing practices, croplands contain on average less biogeochemically labile SOC than forests and grasslands.
Thermal stability, as assessed using T90_HC_PYR, T50_CO2_PYR and T50_CO2_OX, was the highest in vineyards, orchards and croplands compared to 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). 365 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, R-index and I-index (Fig. 2). Cécillon et al. (2021) reported for soils with highly contrasted biogeochemical stability that the relationship between thermal stability and biogeochemical stability was weaker for T50_HC_PYR, R-index and I-index. Our results showed that this relationship even disappears when considering data sets with more heterogeneous topsoil samples. The use of the Rock-Eval® temperature parameters 370 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 to clay content and negatively correlated to 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 much 375 reduced number of study sites. Soil clay fractions interact with microbial compounds which results in the formation of organo-mineral complexes which SOC has a high biogeochemical stability (e.g. Lehmann and Kleber, 2015). We can therefore hypothesize that clay-rich soils are also richer in biogeochemically stable carbon. The positive correlation between clay content and SOC thermal stability, and the good correlations between 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 380 stabilities. Iron oxides are mineral compounds that are also supposed to protect SOC from decomposition. To 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 oxides contents is relatively small in our set of topsoils.
T90_HC_PYR, T50_CO2_PYR and T50_CO2_OX were all positively correlated to pH. Such a correlation between 385 T50_CO2_OX and pH was already 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 acidic bogs. We can therefore hypothesize that acidity slows down SOM mineralization which can favor the accumulation of labile SOC components. As these labile SOC fractions would appear as thermally unstable, it would explain the positive relationship between pH and Rock-Eval® indicators of SOC thermal stability. 390 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 expect SOC labile fractions to be more rapidly processed at higher temperature. It would be in line with the observed positive correlations between MAT and the three selected thermal stability indicators. The relatively weak (Spearman rho 395 value below 0.2) correlations can be due to the fact that MAT will also play on carbon inputs to the soil. Indeed, if higher SOC mineralization was balanced by increased biomass inputs it would mess up dampen the relation between MAT and 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 influences net primary production and therefore the soil carbon 400 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 was obtained on an 8km ✕ 8km grid and do not have the same precision as 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. Massif Central, Alps, Pyrenees) 405 where forest, grassland, and low MAT dominate, and presented by Martin et al. (2011) as with relatively high SOC contents, 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 with warmer MAT, dominated by vineyards and croplands, 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 410 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 repeated application of exogenous organic matter.

Elemental stoichiometry of soil organic matter in French topsoils
Higher values of HI and lower values of OIre6 were observed in forests and grasslands compared to croplands and vineyards. This trend was observed in previous studies (Disnar et al., 2003;Saenger et al., 2013;Sebag et al., 2016). It also 415 confirms that HI and OIre6 can be good proxies of 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) showed a high variability of the C/N ratio between land cover classes: the C/N ratio was higher in forest topsoils than in grasslands, as well as in croplands and vineyards. This is classically explained by 420 the fact that SOC is on average less processed in forests and grasslands compared to croplands and vineyards (Cotrufo et al., 2019), as well as by the higher C/N ratio of the biomass inputs to soil in forests. Indeed, the biotransformation of organic matter tends to lower its C/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 C/N (Fig. 3). similar to those observed for SOM thermal stability. HI and OIre6 are respectively negatively and positively correlated to pH (Table 1) 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 the presence of clays can promote the protection of microbially processed H-depleted SOM. Similarly to what was observed for SOM thermal stability, relationships between elemental stoichiometry 430 and climate variables are weak, probably because climate plays on 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) illustrated the effect of land cover, climate and pedological variables on SOM elemental stoichiometry. Regions dominated by grassland and forest (Fig. 4d5) (Fig. 4) also illustrated the relationships previously observed between these Rock-Eval® parameters (Barré et al., 2016;Cécillon et al., 2021).

Conclusion
This study is an unprecedented effort to make widespread thermal analysis measurements on a national soil quality 440 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 and 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 445 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 of SOC biogeochemical stability, while others parameters could not (T50_HC_PYR, R-index and I-index). 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 PARTYSOC 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 of SOC 450 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 it was observed for the AMG model of SOC dynamics in several French long term agronomic experiments (Kanari et al., 2022).