Soil organic carbon stocks in European croplands and grasslands: How much have we lost in the past decade?

The EU Soil Strategy 2030 aims to increase soil organic carbon (SOC) in agricultural land to enhance soil health and support biodiversity as well as to offset greenhouse gas emissions through soil carbon sequestration. Therefore, the quantification of current SOC stocks and the spatial identification of the main drivers of SOC changes is paramount in the preparation of agricultural policies aimed at enhancing the resilience of agricultural systems in the EU. In this context, changes of SOC stocks (Δ SOCs) for the EU + UK between 2009 and 2018 were estimated by fitting a quantile generalized additive model (qGAM) on data obtained from the revisited points of the Land Use/Land Cover Area Frame Survey (LUCAS) performed in 2009, 2015 and 2018. The analysis of the partial effects derived from the fitted qGAM model shows that land use and land use change observed in the 2009, 2015 and 2018 LUCAS campaigns (i.e. continuous grassland [GGG] or cropland [CCC], conversion grassland to cropland (GGC or GCC) and vice versa [CGG or CCG]) was one of the main drivers of SOC changes. The CCC was the factor that contributed to the lowest negative change on Δ SOC with an estimated partial effect of −0.04 ± 0.01 g C kg−1 year−1, while the GGG the highest positive change with an estimated partial effect of 0.49 ± 0.02 g C kg−1 year−1. This confirms the C sequestration potential of converting cropland to grassland. However, it is important to consider that local soil and environmental conditions may either diminish or enhance the grassland's positive effect on soil C storage. In the EU + UK, the estimated current (2018) topsoil (0–20 cm) SOC stock in agricultural land below 1000 m a.s.l was 9.3 Gt, with a Δ SOC of −0.75% in the period 2009–2018. The highest estimated SOC losses were concentrated in central‐northern countries, while marginal losses were observed in the southeast.


| INTRODUC TI ON
Globally, soils store approximately 2135 Gt of soil organic carbon (SOC) in the first metre (Wang et al., 2022) making them the largest terrestrial reservoir of carbon (Terrer et al., 2021;Todd-Brown et al., 2014).However, factors such as land use and management practices coupled with rising temperatures due to climate change, may potentially trigger a transition of soils from carbon storage to becoming a significant source of atmospheric carbon dioxide (CO 2 ) (Bond-Lamberty & Thomson, 2010;Lal, 2004;Rodrigues et al., 2021).Given the magnitude of this terrestrial carbon pool, even minor shifts in soil carbon stocks, on the scale of a few percentage points, could substantially influence atmospheric CO 2 concentrations and the C balance at global scale (Paustian et al., 2016(Paustian et al., , 2019;;Smith et al., 2020).
In this context, the European Commission strengthens the contribution of the land use, land use change and forestry (LULUCF) sector to the European Union's (EU) increased climate ambition and recognizes the need to reverse the current declining trend of carbon removals (McGrath et al., 2023).The EU Soil Strategy 2030, 'Reaping the benefits of healthy soils for people, food, nature, and climate' (COM(2021) 699 final), recognizes that targeted and sustained sustainable soil management practices can contribute to achieve climate neutrality by increasing the amount of C stored in agricultural soils.
The amended LULUCF legislation proposes an EU target for net removals of greenhouse gases of at least 310 million tonnes of CO 2 equivalent per year in the sector by 2030, distributed among member states as binding targets.Therefore, quantifying current SOC changes and identifying their drivers is paramount in developing targeted strategies and policies to increase SOC stock and enhance the resilience of agricultural systems.Indeed, adopting practices that improve SOC stock has been shown to concurrently enhance the physical attributes of the soils (Watts & Dexter, 1998), reduce susceptibility to erosion and significantly exceed the crop yield performance and stability exhibited by conventionally managed agricultural systems, particularly in periods of drought stress (Lal, 2004;Oldfield et al., 2019;Pan et al., 2009;Paustian et al., 2019;Zhang et al., 2016).However, an incomplete understanding of how land use and edaphic drivers influence SOC changes at different scales (Bispo et al., 2017;Stockmann et al., 2015) adds complexity to the development of EU policies.A large body of research based on mid-to-long-term field experiments has greatly contributed to increase knowledge on SOC dynamics and the development of SOC sequestration rates, particularly following land use change (Richter et al., 2007).
In line with the proposed Intergovernmental Panel on Climate Change (IPCC) Tier 2 approach, CO 2 emission factors (where available) are currently used to monitor and report changes in SOC stocks at the country level.Although those factors, often derived from mid-to-long-term experiments better reflect local conditions, there is a potential issue.Over time, they might be used to estimate linear changes in SOC, even when the actual SOC temporal dynamics display a nonlinear pattern (Lugato, Bampa, et al., 2014), especially involving land use changes.Furthermore, the existing geographical spread of mid-to-long-term experiments and the limited representativeness of soil management practices prevent a comprehensive coverage of all regions in Europe or globally (Bellassen et al., 2022;Smith et al., 2020).Therefore, the conclusions drawn from these point estimates may not accurately reflect larger spatial scales (Mishra et al., 2012;VandenBygaart et al., 2004).
At the highest level of the IPCC Tier approach, Tier 3, efforts are made to overcome the limitations of the lower tiers.This approach utilizes measured-based methods, such as a high spatial resolution soil monitoring system, and validated process-based models to accurately assess SOC changes.However, the Tier 3 approach, despite reducing the uncertainty in estimating SOC changes, is currently only adopted in a limited number of countries due to the high technical expertise requirements and economic burden associated with deploying a high spatio-temporal resolution soil monitoring system (Paustian et al., 2019;Smith et al., 2020).At European scale, spatial estimates of SOC changes have been assessed by employing process-based models (Lugato, Bampa, et al., 2014;Lugato, Panagos, et al., 2014).The high spatial resolution modelled SOC stock estimates were found to be in agreement with aggregated SOC stock at the regional administrative level.However, these estimates were validated with a single point in time and do not allow for the assessment of the model's sensitivity to small detectable changes that generally occur in a short-term time frame (5-10 years).
In the context of inventory-driven Tier 3 method, spatially stratified temporal observational data can be employed to assess changes in SOC stocks over time, particularly following the implementation of management practices in the short to mid time frame (Stockmann et al., 2013).While such data already exist and, in conjunction with data-driven modelling approaches, have been instrumental in assessing SOC changes at national level [e.g. in Argentina (Heuvelink et al., 2021), Hungary (Szatmári et al., 2019), England and Wales (Bellamy et al., 2005), Finland (Heikkinen et al., 2022), Denmark (Taghizadeh-Toosi et al., 2014) and Sweden (Andrén et al., 2008)

| MATERIAL S AND ME THODS
Changes in SOC content, hereafter defined as Δ SOCc (g C kg −1 year −1 ) across EU and UK during the period 2009-2018, were assessed by fitting a quantile generalized additive model (qGAM) (Fasiolo et al., 2021) using the mgcViz R package (Fasiolo et al., 2020) on Δ SOCc calculated from SOC content data obtained from the revisited points of LUCAS topsoil (0-20 cm) surveys of 2009 and 2018.
The trained model was then used to predict Δ SOCc at spatial level using gridded predictors at 500 m resolution.The predicted gridded Δ SOCc were then converted in SOC stocks (Δ SOCs) at 0-20 cm depth.

| Modelling approach and data source
The qGAM was employed to estimate the location-specific median of Δ SOCc across EU and UK.GAM models are semi-parametric regression models that have the ability to capture the nonlinear relationship between response and explanatory variables, via so-called smooth effects.
Standard GAMs require selecting a parametric model (e.g.Gaussian) for the conditional distribution of the response variable.This choice is important because choosing the wrong distribution can lead to invalid p-values for the significance of the effects and to inaccurate estimates of, for example, the conditional mean of the response to be fitted to the data, which in turn requires assumptions about the distribution of the conditional error.This can lead to uncertainties when selecting an appropriate model.In contrast, the qGAMs are non-parametric models that are able to fit the conditional quantile of interest, without making any parametric assumption.This means that it does not require prior assumptions about the conditional distribution of the conditional response variable.This makes qGAMs more robust under heteroscedasticity, which occurs when the variance of the response variable varies with the explanatory variables.Furthermore, qGAMs are fitted by minimizing a penalized pinball loss, which makes them less sensitive to outliers than regression models fitted by penalized least-squares regression.This is important in the context of this study, where the data on Δ SOCc may contain outliers due to natural variation in SOC content, as well as potential errors in data collection or measurement.
The response variable Δ SOCc, expressed as the arithmetic differences between SOC content of 2018 and 2009 was calculated from the revisited points of LUCAS soil surveys collected from specific agricultural lands (i.e.cropland and grassland).The LUCAS survey is a project to monitor land use and land cover changes across the EU.A soil module was performed in 2009, 2015 and 2018 across EU Member States.The LUCAS programme is the largest comprehensive and harmonized source of topsoil data across the EU.For each survey around 22,000 locations spanning various land covers were sampled and analysed for their physical and chemical properties following ISO standard procedures.The sample analyses were performed by a single laboratory, contributing to data comparability by avoiding uncertainties due to analysis based on different methods or different calibrations in multiple laboratories.Further details regarding LUCAS soil module and soil related chemical analysis can be found in Orgiazzi et al. (2018).Of the roughly 22,000 initial locations sampled across all land uses in the 2009/12 soil survey, 5722-specifically from cropland and grassland-were revisited in the two subsequent soil surveys (2015 and 2018).These 5722 locations underwent filtering to exclude sites with organic-rich soils (SOC > 160 g C kg −1 ), those with over 5% CaCO 3 , and those missing either SOC values or particle size analysis, leaving 5307 locations.
Additionally, a meticulous quality control process was implemented to pinpoint changes between 2009 and 2018 unrelated to agricultural practices.This involved inspecting the cover photos from each LUCAS survey (d 'Andrimont et al., 2022).The objective was to identify changes, such as cropland conversions to service roads or discrepancies in sampling centroids across the three surveys (2009/12, 2015 and 2018).From this refined dataset of 4482 observations, the Δ SOCc values were derived.
Since the first LUCAS collection for Romania (RO) and Bulgaria (BG) was completed in the 2012, the total Δ SOCc was divided by the number of years between the first Y t 1 and the last soil survey Y t 2 to normalize the differences in terms of time interval across EU + UK (i.e. 6 years for BG and RO and 9 years for the rest of the countries) using Equation (1).
The effect included in the model for Δ SOCc were: (1) The annual precipitation, precipitation seasonality and SOCc of 2009 were log transformed, while the soil clay content was square root transformed to achieve better dispersion of the observed variables (Wood, 2006).Apart from the median (q0.5), the first (q0.25) and the third (q0.75) quantiles were also modelled to assess the variation of Δ SOCc.The model formula was as follows: The contribution of each variable included in the model's predictions was assessed by visually inspecting partial dependence plots (PD plots).PD plots serve as an effective tool to illustrate the marginal effect of a single predictor variable on the predicted outcome, while holding all other variables constant (Hastie et al., 2009).This technique isolates and visualizes the relationship between a specific variable and the response, providing a clearer understanding of how changes in a predictor correspond to changes in the response variable.The performance of the final model was assessed using tenfold cross-validation.The dataset was randomly divided into 10 equally sized folds, with nine of them being used to train the qGAM model and one to assess the quality of Δ SOCc predictions.This process was repeated 10 times, with each iteration setting aside a different fold, to ensure a robust evaluation of the model's performance.The model's predictions were compared against the observations, and assessed via the root mean squared error, index of agreement (d, Equation 3) and the coefficient of determination (r 2 ).Statistical analyses and graphical presentations (ggplot2; Wickham, 2016) were performed in the R statistical environment (R Core Team, 2021).

| Upscaling model predictions
The trained model was used to upscale Δ SOCc predictions across the EU + UK, using spatially explicit covariates at a resolution of 500 m.The initial SOCc in 2009 (Figure S1) and soil clay content (Figure S2) were obtained from the European Soil Data Centre (Ballabio et al., 2016;Panagos et al., 2022).Long-term MAT (°C, Figure S3), inter-annual temperature variation (Figure S4), annual precipitation (P, mm, Figure S5) and inter-annual precipitation variation were obtained from the high-resolution WorldClim datasets (Fick & Hijmans, 2017).The land cover/use data were obtained from the Corine Land Cover dataset (https:// land.coper nicus.eu/ paneurop ean/ corin e-land-cover ), using the closest years to the LUCAS

| Land use change scenarios
To assess the potential impact of land use change on cumulative Δ SOCs, two different land use change modelled scenarios were examined: converting cropland to grassland (C_to_G) and vice versa (G_to_C).The Δ SOCs estimates obtained from both scenarios were ranked on a scale from 0 to 100, representing the lowest negative Δ SOCs and the highest positive Δ SOCs within each country respectively.For each scenario (C_to_G and G_to_C), cumulative Δ SOCs was calculated for each gradual conversion of land use, using 5% increments, ranging from the lowest to the highest Δ SOCs rank for C_to_G and from the highest to the lowest Δ SOCs rank for G_to_C.
The starting point of the gradual conversion differed for each scenario: for C_to_G, it was set to the lowest rank, in order to maximize the positive effect of grassland on SOC gains in high-loss hotspots, and to identify the minimum cropland area that needs to be converted into grassland to achieve a net balance of gains and losses across the study area.For G_to_C, the starting point was set to the highest rank, in order to minimize the negative effect of cropland expansion on SOC stocks.

| Modelling Δ SOC using qGAM
The revisited LUCAS points in Figure 1 showed Δ SOCc values ranging from −0.83 to 0.98 g C kg −1 year −1 (q0.05 and q0.95 respectively), with a median of 0.004 g C kg −1 year −1 .Notably, continuous grassland (GGG) and the transition from cropland to grassland (CGG) displayed the highest median Δ SOCc values of 0.3 g C kg −1 year −1 and 0.28 g C kg −1 year −1 respectively.Conversely, continuous cropland (CCC) and the transition from grassland to cropland (GCC) exhibited the lowest median Δ SOCc values of 0.001 g C kg −1 year −1 and −0.01 g C kg −1 year −1 respectively.
The qGAM accurately predicted Δ SOCc with a mean absolute error of 0.3 g C kg −1 year −1 (Figure 2), demonstrating its reliability.
The high d value of 0.70 and r 2 value of .40further confirmed the good correlation between predicted and observed values.The partial effect plots of the qGAM covariates (Figure 3) revealed that continuous grassland (GGG) had a significant positive influence on Δ SOCc, contributing an estimated 0.49 ± 0.02 g C kg −1 year −1 (p < .001, Figure 3a).The transition from cropland to grassland also had a positive effect, although only significant for CGG with a contribution of 0.32 ± 0.03 g C kg −1 year −1 (p < .001).Conversely, continuous cropland (CCC) showed a significant negative influence on SOC changes, with an estimated contribution of −0.04 ± 0.01 g C kg −1 year −1 (p < .001).
The combined effect of initial SOC content (SOC 2009 ) and soil clay content is shown in Figure 3b.Higher soil clay content (2) at lower SOC 2009 content had a positive effect on Δ SOCc, while increasing SOC 2009 content had a negative influence on Δ SOCc.
Additionally, for the MAT covariate (Figure 3c), a negative influence up to −0.25 g C kg −1 year −1 was observed for temperatures above 12°C with low inter-annual variability.In contrast, high annual mean precipitation (P, Figure 3d) at any level of inter-annual variability (P_CV) had a positive influence, providing a gain of up to 1 g C kg −1 year −1 .

| Spatial prediction of Δ SOCs across Europe
Total

| Land use change effect on SOC changes across Europe
The ranking of Δ SOCs within each European country for the conversion of grassland to cropland (G_to_C), cropland to grassland (C_to_G) and the estimates of total changes in SOC stock for each gradual conversion are presented in Figure 5.The gradual conversion of grassland to cropland starting from the top 5% best performing grasslands (i.e.rank >95) always provided a net SOCs losses (Figure 5a,b).The relative SOC losses in G_to_C conversion scenario went from −0.8% if the top 5% grasslands would be converted (2 Mha) to −4.8% if the total current area of grassland (45 Mha) would be converted to cropland in a 9 year period.In contrast, the gradual conversion of cropland to grassland starting from the worst 5% performing cropland area (i.e.rank <5) provided a net increase in SOC stocks (Figure 5c,d).With the conversion of the worst 7.5% croplands within each country for a total of 8.6 Mha a net 0 Δ SOCs across Europe could be achieved in a 9-year period.With the further conversion of cropland to grassland an increase up to 12.1% of SOC stock was estimated.

| Drivers of SOC changes across European agricultural soils
Changes in SOC content by land use class between 2009 and 2018 (Figure 1) confirm the beneficial gain in SOC content obtained from continuous grassland or the conversion of cropland to grassland.Continuous grassland (GGG, Figure 3) or the conversion of cropland to grassland (CGG) could contribute to an increase in SOC content by 0.48 ± 0.01 and 0.33 ± 0.04 g C kg −1 year −1 , (~1.2 and 0.8 t C ha −1 year −1 respectively) across European pedoclimatic conditions.This aligns with previous modelling studies (Lugato, Bampa, et al., 2014;Poeplau et al., 2011;Schneider et al., 2021) and meta-analysis (Conant et al., 2017;Deng et al., 2016), which found that grassland can increase SOC stocks between 0.1 and 1 t C ha −1 year −1 .A recent ensemble modelling of land surface-atmosphere CO 2 flux across Europe (McGrath et al., 2023) revealed that grasslands act as net sinks, while croplands are net sources of CO 2 .Although there is agreement in trends between the present study and the ensemble models results, the magnitudes of these roles differ significantly.
In grassland, the majority of the C inputs are allocated belowground via root degradation and C-rich roots exudates (Jackson et al., 2017).These belowground C inputs are more likely to be incorporated in the mineral-associated organic matter (MAOM), a stable form or soil organic matter with a residence time in the soil that spans decades to centuries, therefore contributing to long-term carbon sequestration in soil (Bai & Cotrufo, 2022).In contrast, qGAM partial effect analysis highlighted a net negative influence on SOC content for continuous cropland (CCC) up to −0.03 ± 0.005 g C kg −1 year −1 .
Continuous cropland, especially monoculture, or the conversion of grassland to cropland (e.g.GCC), even under conservation agriculture practices, have been identified as major detrimental factors for SOC storage in agricultural systems (DuPont et al., 2010;Wiesmeier et al., 2019).
Besides the lower soil C inputs of croplands compared to grasslands (Janzen et al., 2022), SOC losses in cropland have been attributed to a faster SOC turnover due to increased aeration (that promotes soil organic matter oxidation) and the disruption of soil aggregates, both promoted by tillage events (Herzfeld et al., 2021).The latter is not a common practice in grassland.However, at a European scale, SOC losses have also been estimated for long-term continuous grassland, such as in Ireland (Figure 4), where SOC content was estimate to be above 150 t C ha −1 or >50 g C kg −1 (Figure 5).The qGAM partial effect plots showed a negative effect on SOC changes up to −2 g C kg −1 year −1 at high initial SOC content (SOC 2009 ) (Figure 3b), which could potentially be explained by saturation of stable organic matter (MAOM) and subsequent accumulation of more labile C in the particulate organic matter fraction more sensitive to perturbation (Cotrufo et al., 2019).
However, the amount of SOC that can be stored and protected in the soil profile is not only constrained by the available mineral surface area (Chenu et al., 2019;Derrien et al., 2023;Lavallee et al., 2020;Wiesmeier et al., 2019) but also by other soil edaphic and environmental factors.In this study, it was found that the combination of low initial SOC content and high soil clay content positively influenced Δ SOCc (Figure 3).Conversely, in regions characterized by low soil clay content, high SOC losses have been estimated, including in countries with favourable climatic conditions for SOC sequestration (Figure 3).Favourable climatic conditions for SOC accumulation are generally identified as humid and relatively low temperatures (Wiesmeier et al., 2019).In this study, long-term MAT <13°C and annual precipitation above 700 mm have a positive influence on SOC between 2009 and 2018 (Figure 3).Contrary, warm and semi-arid climatic conditions (i.e.Mediterranean climate) promoted SOC losses (Figure 3c).The patterns of these effects align with findings from previous studies conducted across Europe (Rial et al., 2017;Schneider et al., 2021).While MAT is an important factor in determining SOC changes, it alone does not fully explain the spatio-temporal trend of losses across Europe.Instead, the combination of MAT and temperature seasonality has a greater impact on SOC changes.It is not surprising that high and steady MAT has a negative effect on SOC stocks.However, low MAT (<10°C) combined with high levels of temperature seasonality, were also found to have a negative effect on SOC.These findings are consistent with a recent assessment of SOC trends in Finland (Heikkinen et al., 2022), where increased temperature during the summer period was identified as one of the drivers of SOC losses.Indeed, for Nordic countries (Denmark, Finland, Sweden and Estonia) even the q0.75 model (Figure 4) estimated SOC losses.These estimates are in agreement with similar studies conducted at national level (Bellamy et al., 2005;Lettens et al., 2005;Meersmans et al., 2011;Taghizadeh-Toosi et al., 2014).This estimated decline in SOC, is concerning given that these countries have the highest SOC stocks, and therefore large potential for losses.
Projected climate scenarios suggest that by 2050, MAT will increase by 1-2°C and there will be extended periods of drought.
These conditions could lead to a monotonic increase in SOC losses into the atmosphere, particularly at high latitudes, exacerbating the negative feedback on climate change.
It is important to note that soil systems have a finite capacity for storing SOC.For instance, observed increases in SOC under specific land management scenarios, such as the conversion of cropland to When only a single set of data in time is available, data-driven approaches to estimate SOC changes at continental scale tend to derive a static map of SOC employing geostatistical models and then use the derived model covariates coefficients to project SOC predictions in time (e.g.Yigini & Panagos, 2016).This approach assumes that covariates, fitted to predict the spatial distribution of the dependent variable, are also fit to predict the temporal trend.
However, this is not always the case, especially when SOC changes are relatively small in a relatively short time frame (e.g.3-5 years).
For instance, estimates of the spatial distribution of the dependent variable (in this case, large background SOC stocks) have an associated uncertainty.If the associated uncertainty is larger than the naturally occurring slow C changes in short time frame, it becomes clear that final estimates of the changes are affected by a large uncertainty.Indeed, when estimating perceivable changes in SOC stocks at the field, region or continental scale, an added complexity arises from detecting relative small SOC changes that are driven by climate within the considered time frame (Smith et al., 2020;Stanley et al., 2023).This study, by leveraging on the temporal component For the current study, soil bulk density, also corrected for the rock fragment content, was estimated using an empirical pedotransfer function, as the bulk density in LUCAS was not measured in 2009.Frequently, national and continental soil inventories estimate soil bulk density using established pedotransfer functions (Poeplau et al., 2017).The widely recognized pedotransfer function from Hollis et al. (2012) was used and its predictive capacity was further validated.Despite the large variability within the measured soil bulk density, the values predicted by this function aligned strongly with those measured in the 2018 LUCAS soil survey (Figure S7).However, it is important to note that soil bulk density measurements are subject to random errors that can reach up to 40% (Zhou et al., 2023).Furthermore, the LUCAS sampling depth is confined to the 0-20 cm depth range, which does not account for dilution effects in SOCc attributed to tillage events that can reach depths of up to 0-30 cm in croplands.Despite these limitations, the LUCAS soil dataset remains the only harmonized soil survey at a European scale.Future refined LUCAS surveys are expected to address these limitations, improving the accuracy and reliability of soil data at the European level.
Multi-model ensemble approach proved to be an effective way to reduce the uncertainty in climate and crop modelling (Asseng et al., 2013;Martre et al., 2015).The benefit of multi-model ensemble approach has also been proven for SOC estimates at spatial scale.to the way each model simulates dynamic processes, sets parameters and uses input variables (Asseng et al., 2013;Riggers et al., 2019;White et al., 2011).Ensemble data-driven approaches have also shown promising results in estimating SOC (Biney et al., 2022;Szatmári et al., 2019) as well as in predicting soil mineral-associated organic carbon fractions (Xiao et al., 2022).Although modelling approaches have demonstrated their ability to effectively reproduce observational data, questions remain regarding their transferability and scalability where no or scarce data are available at the continental scale.The reliability of models and their parameterization and verification require comprehensive, harmonized and temporally resolved datasets that to date, mainly exit in temperate zones.
To this end, comprehensive information from under-represented zones can be obtained from detailed spatial soil surveys coupled with long-term field experiments.This information could enable the development of accurate and robust models that can effectively reproduce observational data and inform decision-making in the context of sustainable land management from continental to regional scale.
However, the transferability of such modelling frameworks from regional to field scale is questioned due to the inverse relationship between spatial scale and the uncertainty in SOC estimate (Oldfield et al., 2022).In fact, Ogle et al. (2010)

| Strategic land use change effect on SOC changes
Grassland is widely recognized as the agricultural practice that promotes the highest SOC storage.In fact, converting current cropland to grassland could increase SOC by more than 10% within a 9-year period across Europe.In contrast, converting grassland to cropland, which encompasses 47 Mha, could lead to a net loss of up to 4.8% of current SOC stocks.Despite the scenario of converting 116 Mha of cropland to grassland is rather utopic, as food security targets through crop production must be maintained in the EU, the incremental conversion approach can provide useful indication under different socio-economic scenarios.For instance, the proposed strategic approach of only converting the top 7.5% highrisk SOC losses cropland within each Member State, could balance SOC losses and gains across the EU ('Net Zero Emissions').This approach would enable the preservation of cropland necessary for food security while promoting SOC storage in areas with the greatest potential for success.For effective climate change mitigation and emissions offsetting, it is necessary to maintain grassland converted from cropland over time (beyond the time frame considered in this study), as this is crucial for SOC stocks to remain at their maximum levels following the establishment of a new equilibrium (Moinet et al., 2023).However, the applicability of this measure needs to take into account its direct economic impact on farms business, as farmers generate higher profits from using arable land for crop production (Gocht et al., 2016).Income from private soil carbon sequestration initiatives is unlikely to bridge this financial gap as they cannot ensure permanence beyond the length of the contract (Paul et al., 2023).Over time, it is expected that the rates of SOC sequestration will gradually decrease until steady state is achieved (Minasny et al., 2017), which can lead to reduced carbon credits for farmers who participate in carbon crediting schemes.Consequently, the income generated from SOC sequestration schemes is likely to decrease over time.Furthermore, SOC sequestration rates are constrained by various environmental factors, such as temperature and rainfall (Figure 3c,d), as well as site-specific factors such as soil clay content (Figure 3b), which cannot be altered by management practices.Therefore, crediting farmers solely based on the absolute amount of carbon sequestered, apart from the uncertainties associated with measuring, reporting and verification (previous section), will lead to a disparity in opportunities across the continent and hinder the homogenous uptake of the practice at a continental scale.
But, if the pricing of sequestered CO 2 is not merely dependent on the absolute amount, but also includes additional environmental co-benefits and considers site-specific environmental constraints, the crediting system could be steered towards greater balance and impartiality.Therefore, finding ways to incentivize long-term maintenance of permanent grasslands across the whole continent, while also addressing the economic challenges faced by farmers, will be crucial for the success of this emissions mitigation strategy.
The proposed spatial mitigation strategy has limitations in that it only examined the impact of pedoclimatic conditions and land use changes on SOC stocks from a general and simplistic perspective, as it did not consider the effect of other specific agricultural soil management.Within this context, it is acknowledged that some scenario representations, especially those involving extreme land use changes, are theoretical constructs meant to highlight potential impacts rather than propose direct policy actions.However, the proposed strategic approach, which considers the inherent SOC sequestration constraints of specific regions, could help policymakers identify areas where the implementation of mitigation practices would be most effective.Empowering policymakers with such spatial information could facilitate the delineation of targeted spatial ac- ] they have not yet been employed to assess SOC changes at the European continental level.A comprehensive and harmonized pan-European analysis is key to understanding the combined effects of land use and climate that transcend national borders.Given the previously unexplored nature of a pan-European assessment, this research primarily aims to leverage spatially stratified temporal observational data to investigate changes in SOC stocks of croplands and grasslands at the continental scale across Europe.Achieving this broader view is particularly critical as the insights derived will directly inform and assist policymakers and stakeholders in evaluating the objectives of the European Green Deal.For the purposes, in this study, a data-driven spatial SOC changes modelling approach was employed, utilizing SOC temporal change data obtained from the Land Use/Land Cover Area Frame Survey (LUCAS) conducted across EU member states in 2009, 2015 and 2018.This novel approach provided for the first time estimate of SOC changes at the European continental level based on observational stratified spatio-temporal data.This methodology offers the potential to further elucidate how climate, land use and management and edaphic factors influence SOC changes across the EU and the United Kingdom (UK).
(i) Linear effect of land use information for each revisited LUCAS point in 2009, 2015 and 2018 for cropland (C) and grassland (G).These data provided useful information to assess the effect of land use and land use change, that is, continued grassland (GGG) or cropland (CCC), grassland to cropland (GGC or GCC) and vice versa (CGG or CCG) on Δ SOCc.(ii) Smooth, nonlinear interaction of SOC content (g C kg −1 ) and soil clay content (%) at 0-20 cm depth in the 2009/12, obtained from the 2009/12 LUCAS soil data.(iii) Smooth, nonlinear interaction of annual long-term mean precipitation (P, mm) with precipitation seasonality (coefficient of variation, P_CV) and of annual long-term temperature (mean annual temperature [MAT], °C) with temperature seasonality (standard deviation, T_SD) extracted from the WorldClim climatic datasets (Fick & Hijmans, 2017).The WorldClim datasets have average monthly climate data for minimum, mean and maximum temperatures and for precipitation for 1970-2000 at a resolution of 1 km.
sampling surveys and resampled to a resolution of 500 m.The predicted values were then converted in Δ SOC stocks (Δ SOCs) at 0-20 cm depth.The soil bulk density used to convert SOC concertation to stocks was predicted using an empirically derived pedotransfer function, adapted fromHollis et al. (2012) and corrected for coarse fragments content(Poeplau et al., 2017) similar to Lugato, Panagos, et al. (2014) and Schneider et al. (2021).Comparison between predicted adjusted soil bulk density and measured soil bulk density obtained from LUCAS soil survey 2018 is shown in Figure S7.
Figure4c) quantiles for grassland and cropland for the period

F
I G U R E 2 Density scatter plot of predicted against observed Δ SOCc (g C kg −1 year −1 ) at 0-20 cm for 2009-2018.Solid line is the 1:1 line.d, index of agreement; r 2 , coefficient of determination; RMSE, root mean squared error; SOC, soil organic carbon.

F
Partial dependence plots for predictive variables included in quantile generalized additive model.(a) Land use 2009-2015-2018 (C, cropland; G, grassland).(b) Interaction of initial SOC content (SOC 2009 ) and soil clay content.(c) Interaction of annual long-term temperature (MAT, °C) with temperature seasonality (standard deviation, T_SD).(d) Interaction of annual long-term mean precipitation (P, mm) with precipitation seasonality (coefficient of variation, P_CV).SOC, soil organic carbon.MAT, mean annual temperature.

F I G U R E 4
Spatial prediction of Δ SOCs (0-20 cm) using qGAM during the 2009-2018 period for (a) q0.5, (b) q0.25 and (c) q0.75.(d) Latitudinal variations of mean Δ SOCs (0-20 cm) for qGAMs and DayCent.Map lines delineate study areas and do not necessarily depict accepted national boundaries.qGAM, quantile generalized additive model; SOC, soil organic carbon.grassland, may plateau as the soil system approaches this equilibrium.Notably, SOC changes exhibit a sinusoidal function between new and prior steady states.Recent literature, including Gocke et al. (2023), suggests ambiguity about when this equilibrium state is attained, and even the 20-year time frame suggested by the IPCC (2006) may be inadequate for capturing these complex dynamics.In the context of short-to-medium-term studies/surveys, like LUCAS surveys from 2009 to 2018, calculating the delta from repeated data points serves as an insightful methodology, capturing both steadystate and non-steady-state conditions of SOC.Addressing the multi-dimensional complexities that influence SOC, this study strategically combined static and dynamic drivers (land cover/use change) capturing cumulative effects observed over the 9-year period related to complex anthropogenic and environmental driver interactions.F I G U R E 5 Ranking of SOCs changes from the worst (0) to the best (100) for each country during the conversion of grassland to cropland [G_to_C (a)] and cropland to grassland [C_to_G (c)], along with the corresponding total SOC stock changes [G_to_C (b), C_to_G (d)] resulting from gradual land use conversions in 5% increments.The ranking ranges from the lowest to the highest Δ SOCs for C_to_G and from the highest to the lowest Δ SOCs for G_to_C.Map lines delineate study areas and do not necessarily depict accepted national boundaries.SOC, soil organic carbon.4.2 | Uncertainty in data-driven Δ SOC modelling across Europe The qGAM model was trained and validated (Figure 2) with spatiotemporally extended data collected on specific revisited locations across Europe between 2009 and 2018.The revised soil sampling regime of the LUCAS soil module, for the first time, reduced the uncertainty on data source for model training at continental scale.
of LUCAS soil survey, directly modelled SOC changes.Modelling directly SOC changes reduced the possible error associated with covariates influence, and reduced the uncertainty in estimating small SOC changes.The estimated changes in SOC (Δ SOCc) over time attributed to initial SOC (SOC 2009 ) were relatively small, ranging from −2.5 to 1.3 g C kg −1 year −1 (Figure3b).In contrast, the background of SOC content (SOC 2009 ) can vary significantly, ranging from 1.5 to over 50 g C kg −1 .Given the small magnitude of SOC change relative to the much larger background, even a minor error in estimating the background of SOC could significantly impact the accuracy of estimating the change over time.Overall, the Δ SOCs spatial estimates obtained with qGAM model were in agreement with Δ SOCs estimates obtained with the tested mechanistic model DayCent.Large discrepancy between models was mainly observed at the extreme latitudes within the mean ΔSOCs distribution (Figure4d).The observed discrepancies between model predictions are attributed to the training dataset variance, distribution and size for data-driven approaches, while for mechanistic models to the inherent characteristics of the model itself.Data-driven approaches may face potential limitations due to biases present in the training dataset, such as unevenly distributed sampling locations, and the absence of temporal measurements for certain variables, such as soil bulk density.
Riggers et al. (2019) combined and evaluated predictions of six different model concepts across Germany and found that the ensemble model predictions outperform the single model predictions.Multimodel ensemble approach reduces the structural uncertainty linked tion plans to achieve compliance with the latest LULUCF regulation, which requires each Member State to offset accounted emissions from land use with equivalent removals in the sector.