Grassland contribution to soil organic carbon stock under climate change scenarios in Basque Country (Spain)

In this study, we estimated the contribution of managed grasslands to the “4 per 1000” initiative in Basque Country under two climate scenarios (RCP4.5 and RCP8.5) adopted by the IPCC in its Fifth Assessment Report. For this purpose, the RothC model was calibrated and validated with a historic database of grassland soil organic carbon (SOC) (1983 to 2019). The results at field-scale show a rate of increase of 1.26 t C ha−1 year−1 at a depth of 0–30 cm over a 36-year simulation. The model was run at the regional scale in short- (2020–2040), medium- (2041–2070), and long-term (2071–2090) future climate scenarios. For all the simulations, agricultural practices and available data for grassland systems were considered. RothC model projections showed how SOC stock responses varied depending on initial SOC and climate subregions, with higher values for a lower initial SOC stock under the highest precipitation regime subregion. A 4 per 1000 storage rate could be achieved in grassland soils with an initial SOC < 80 t C ha−1. The overall trends showed that future climate change will lead to a decrease in the SOC stock in grasslands with a higher initial SOC if appropriate practices are not implemented to maintain the SOC stock.


Introduction
Soil organic carbon (SOC) constitutes a part of the carbon (C) cycle, incorporating atmospheric C and acting as a sink of C (FAO 2017). The "4 per 1000-Soils for Food Security and Climate" initiative (hereinafter referred to as the "Initiative") was launched with the objective of increasing global SOC stocks by 0.4% per year over 20 years as compensation for global greenhouse gas emissions (GHGs) due to anthropogenic sources (Minasny et al. 2017). Previous studies have indicated that adopting recommended management practices to maintain the SOC concentration could potentially increase the global stock by 1.2-3.1 Pg C year −1 over 25 to 50 years (Lal 2013).
Grasslands are one of the most important and widespread terrestrial ecosystems, covering approximately 40% of the global land surface and containing the largest share (39%) of terrestrial soil C stocks (White et al. 2000). SOC comprises 47% of the mitigation potential for agriculture and grasslands (Bossio et al. 2020). Grassland systems in Basque Country are an important resource for livestock production. They represent 45% of the total cultivated area (Departamento de Desarrollo Económico 2020), thus offering considerable potential as C sinks through the connection between livestock and grasslands for climate change mitigation. In this region, sustainable grassland intensification aims to increase plant yields destined for forage production while maintaining soil function as a C sink. The remaining challenge is to reliably quantify rates of change in soil C across intensively managed grasslands.
Climate change is a key driver of change in soil C dynamics (Gottschalk et al. 2012). Increasing mean air temperatures are likely to accelerate SOC decomposition and loss of SOC in the future if soil moisture is not a limiting factor. However, the magnitude and the duration of this effect on soil acting as a C sink in intensive livestock systems are uncertain (Fornara et al. 2020).
Due to the complexity of the soil-plant-atmosphere system, it is helpful to study SOC dynamics through simulation models combined with local measured data. SOC turnover models have been used to simulate SOC dynamics under various land uses, management practices, and climatic conditions (Smith et al. 1997b;Smith 2002). Different soil-process models have been developed and applied to understand soil functions (e.g., DayCent-Parton et al. (1994); DNDC- Li et al. (1992); RothC- Coleman and Jenkinson (1996)). The RothC model is a process-based model to estimate SOC dynamics that has been widely used in a variety of ecosystems (Coleman et al. 1997;Smith et al. 2008;Xu et al. 2011) and climate regions (Falloon and Smith 2006;Powlson et al. 2011;Meyer et al. 2018). It has also been used to predict future SOC dynamics under the effect of climate change (Gottschalk et al. 2012;Chenu et al. 2019;Afzali et al. 2019). SOC modelling under climate change conditions has also been previously studied in Spain at a regional scale (Nieto et al. 2010;Segura et al. 2016;Jebari et al. 2018), but there are no modelling studies specific to the C stock in grasslands in this region. In addition, most of the modelling efforts have focused on forests and arable lands, while grasslands have received less attention (Soussana et al. 2006;Wiesmeier et al. 2019) due to the complexity of the process-based modelling of grassland systems compared to annual cropping systems.
The objectives of this study were to (1) calibrate the RothC model with long-term data on intensively managed grasslands for livestock production at the field scale and (2) evaluate the long-term (2020-2090) contribution of the Basque Country grassland soils to the Initiative under two climate change scenarios (RCP4.5 and RCP8.5) adopted by the IPCC in its Fifth Assessment Report (IPCC 2013).

Study area
Basque Country is situated on the north coast of Spain and has a highly diverse climate. On the northern slope, the dominant climate is Atlantic (Jones et al. 2020) and has led to the development of pastures and forests for livestock and timber production, while on the southern slope, a Mediterranean transitional climate is suited for crop production and dominated by extensive crops.
Two studies were performed, one at a field scale using long-term historical soil data to validate the functionality of the RothC model. Information obtained from this work was used to analyse the effect of climate change on SOC dynamics at the regional scale.

Field-scale approach
Data from 1983 to 2019 from grassland fields located in Basque Country (42°59′0.04"N, 2°36′27.8"W) were used to run the RothC model (Fig. 1). The regional climate is typically Atlantic (Jones et al. 2020), with the area of study being a hybrid zone transitioning to a Mediterranean climate, with a mean annual (1971-2016) temperature of 11.4 °C and a mean annual rainfall of 1410 mm. The site has a mean elevation of 560 m above sea level and an area of 22,000 hectares, and the soils are mainly Acrisols. The clay content of the soils ranges from 18.2 to 34.5%, and the bulk density ranges from 0.94 to 1.41 g cm −3 (0-30 cm) (Supplementary file. Section 1, Tables 3 and 4).
The sampled soils were part of a monitoring programme established to determine the soil fertilization recommendation at 31 sites established in the area, of which 21 were selected for the study due to having the same management practices and conditions. Datasets covered a period of 36 years and included full daily meteorological records and at least two measurements of soil organic matter content.
Soil sampling was carried out following standard methods as required for the fertilization assessment (details of data sampling are given in Supplementary file. Section 1). For the last sampling event, replicated soil samples were taken at a depth of 0-30 cm in 3 geolocated random plots on each site with a Monterra™ Garmin. Soil texture was determined by laser diffraction with a Mastersizer 2000, organic matter (%) was measured according to Walkley and Black (1934), and soil bulk density was measured using the core method (Grossman and Reinsch 2018). The most recently measured (2018-2019) SOC stock ranged from 77.27 to 122.14 t C ha −1 . These data agree with those reported for Basque Country by Calvo de Anta et al. (2020), with grassland soils (0-30 cm) recorded at 103 ± 61 t C ha −1 .
Site locations were chosen to provide representative land use and management practices for forage production. Plant composition was predominantly composed of two botanical groups: Lolium spp. and Trifolium spp. Based on surveys with farmers, we established, for modelling purposes, that grasslands were fertilized with slurry twice a year (5.3 t C ha −1 ), in April and September, and the forage yield averaged 9.5 t DM ha −1 .

Regional-scale approach
A regional-scale approach was undertaken to study the effect of climate change on the SOC in grasslands of Basque Country. Basque Country has a surface area of 7234 km 2 , of which 24% of land is occupied by grasslands. The soil C stock map of Basque Country showed a range from 20.2 to 179.4 t C ha −1 for grassland soil, where most of the soils with a lower SOC stock were found in the southern part of the region. This map shows soil organic C and texture from the first 30 cm of soil (details are given in Supplementary file. Section 2).

Climate data sources and projections
Historic climate data (monthly precipitation and air temperature) from 1971 to 2019 for the field-scale study area were obtained from the Basque Agency of Meteorology. 1 To begin simulations, climate data from at least 12 previous years were averaged. From 1983 onwards, annual climate data were used. The monthly open-pan evaporation was calculated from potential evapotranspiration (Hargreaves-Samani equation) using the ETo calculator, developed by the Land and Water Division of FAO (Allen et al. 1998) and dividing it by 0.75, as suggested by Coleman and Jenkinson (1996).
To obtain a climate classification of Basque Country for the regional simulations, an isocluster unsupervised classification was conducted introducing average daily temperature, average annual precipitation, and annual ETo for 1971-2016 as input bands and using 3 class differentiations and a minimum class size of 20 cells. The three resulting climate subregions are shown in Fig. 1. Subregion 1 is located in the southern part of Basque Country and is characterized by the lowest temperature and rainfall values. Subregions 2 and 3 are located on the coast, and they do not differ in temperature, while subregion 3 is the most humid (Table 1).
The climate projections used in this study were obtained from KLIMATEK 2016. Data are available online 2 and include historical data and projections for three periods. On the one hand, historical data on daily temperature and precipitation for the period 1971-2015 were obtained by interpolating historical data from 180 meteorological stations distributed throughout Basque Country. On the other hand are climate projections based on the EURO-CORDEX Fig. 1 Study area at the fieldscale approach (hatched area) and subregions selected based climate parameters (yellow area, subregion 1; green area, subregion 2; blue area, subregion 3) initiative, which provides regional climate projections for Europe up to 0.11° resolution (EUR-11, ~ 12.5 km) as part of the global CORDEX framework (the World Climate Research Programme Coordinated Regional Downscaling Experiment). 3 In Basque Country, thirteen simulations have been used for RCP4.5 and for RCP8.5, combining six different RCMs and five different GCMs (more detail in Supplementary file, Sect. 3, Table 5). These thirteen EURO-CORDEX simulations were used by combining the delta method (Zahn and von Storch 2010) with historical data to achieve climate projections at a resolution of ~ 1 × 1 km 2 .

RothC model description
The Rothamsted Carbon Model (RothC-26.3) was used to simulate SOC dynamics (Coleman and Jenkinson 1996). RothC-26.3 is a process-oriented model for the turnover of SOC in nonwaterlogged topsoil that runs at a monthly time step.
As described by Coleman et al. (1997), the RothC model separates SOC into four active fractions: decomposable plant material (DPM), resistant plant material (RPM), microbial biomass (BIO), humified organic matter (HUM), and a small amount of inert organic matter (IOM) that is resistant to decay. Each fraction, except for IOM, undergoes decomposition by first-order kinetics at its own unique rate, which is determined by using modifiers for soil moisture, temperature, and plant cover, which are the main input parameters for running the model. The initial IOM pool was set using the equation proposed by Falloon et al. (1998).
The model can be run in 'inverse' mode when inputs are calculated from known changes in soil organic matter or in 'forward' mode using known inputs to calculate the changes in soil organic matter. A more detailed description of the model is given by Coleman and Jenkinson (1996).

Field-scale approach
The model was first calibrated for fields selected from the study area ( Fig. 1) considering land management and historical SOC data as model inputs. The model performance was then validated using data available from the rest of the grassland sites.
RothC was first run to equilibrium in inverse mode (i.e., the C input was adjusted so that the modelled SOC matched the initial measured SOC of 27.13-81.59 t C ha −1 ). We assumed that the soils were in equilibrium with regard to SOC. To establish this initial SOC, we used inverse mode simulation, determining that grasslands need a monthly plant C input to the soil of 0.11-0.33 t C ha −1 . Once the starting C content had been established, the model was used to predict changes in the organic carbon in 19 plots for model validation.
The DPM/RPM ratio or litter quality factor was set to the default value for grasslands (1.44). A value of 0.25 was used for the first year of the series for those fields converted to grassland from forest.
The average measured grassland yield in the area of study (9.5 t DM ha −1 ) was used to calculate the monthly plant input, which was distributed across the different months depending on the soil cover. The contribution of C to the roots was calculated based on the shoot-to-root ratio described by Bolinder et al. (2002) for grasslands, yielding a proportion of 1.5 considering a shoot C content in grasslands of 42%. After setting these equilibrium conditions, the model was run in forward mode  to adjust the inputs for plants and farmyard manure (FYM). The FYM inputs were 66.54% of the total C incorporated into the soil annually.

Regional scale
A total of 16,207 grassland plots (857.44 km 2 ) were selected to carry out regional simulations. For the land management file, we used the same plant C inputs per month for subregions 2 and 3 (0.25 t C ha −1 ) since they were associated with Atlantic production of grasslands in Basque Country with an Table 1 Climatic parameters at the field-scale and regional-scale in the current and final periods (2071)(2072)(2073)(2074)(2075)(2076)(2077)(2078)(2079)(2080)(2081)(2082)(2083)(2084)(2085)(2086)(2087)(2088)(2089)(2090)  average of 10.85 t DM ha −1 . For subregion 1 (under the influence of a Mediterranean climate), the grassland yields were lower (4.93 t DM ha −1 ), and the plant C inputs were set as 0.12 t C ha −1 month −1 (Oyanarte and Rodríguez 1993). The C input from the FYM application for the 3 different climatic zones was set the same as that in the farm-scale study. We simulated SOC changes for the period 2020-2090 under two climate change scenarios (RCP4.5 and RCP8.5). Three different periods were simulated to analyse the short-(2020-2040), medium-(2041-2070), and long-term (2071-2090) effects of the different climate projections. Soils were classified into 5 ranges to detect SOC changes: < 50, 50-79, 80-109, 110-139, and > 140 t C ha −1 . We automatized the total number of simulations of grassland sites for the period of 2020 to 2090 using the RothC Microsoft Excel version provided by Meyer et al. (2018). For the regional-scale study, soil, climate, and land use data were processed in GIS format (ArcMap 10.6). We assumed constant plant and fertilizer C inputs without considering the effect of future climate. We further tested the sensitivity of the model to C input from plants and FYM. To reflect the variability in C inputs from plants within a subregion, simulations were performed considering the range of grassland yield production data from previous regional studies to estimate plant residues (Oyanarte and Rodríguez 1993). Although grassland yield production was not available for each spatial unit, we included the variability resulting from using the maximum, median, and minimum C inputs from plants. Spatial data on organic amendment application were not available; thus, we assumed a 20% variation in carbon inputs derived from animal manure over the annual mean value. A total of 16,207 units × 2 scenarios (RCP4.5 and RCP8.5) × 2 parameters (C input from plants and from FYM) simulations were performed.

Statistical analysis
We compared model output in relation to a set of data that had not been calibrated. Simulated and observed data were statistically analysed to determine the degree of fit. Methods for evaluating the accuracy of a simulation are discussed in detail elsewhere (Smith et al. 1997a). The quantitative calculations described below provide information on a distinct aspect of the accuracy of the simulation. The coefficient of determination (R 2 ), root mean square error (RMSE), standard deviation (SD), Nash-Sutcliffe efficiency model (EF), and average relative error (ARE%) were calculated. The R 2 , RMSE, and EF are statistical tests to determine the level of model fit for the observed and simulated data. The RMSE value ranges from 0 to a positive value, which should be on the same order of magnitude as the standard deviation (Willmott and Matsuura 2005); the closer it is to 0, the closer the simulation is to reality. The maximum value of EF is 1, indicating a perfect fit (Nash and Sutcliffe 1970). Finally, ARE% was calculated to quantify model overestimation (positive) or underestimation (negative) (Yang et al. 2014).
where Pi is the simulated value, Oi is the observed value, O is the average of the observed values, P is the average of the simulated values, and n is the number of measured values.
A statistical analysis using the least significant difference (LSD) test was performed for regional-scale climate change scenarios.

Field-scale model performance
Once the model was set at equilibrium, the simulated SOC was similar to the measured value, with a range of 27.13 to 81.59 t C ha −1 at 30 cm depth (Supplementary file, Table 3). The model estimated SOC values of 73.86 to 112.58 t C ha −1 in 2019. The field-specific temporal dynamics of C accumulation are given in Fig. 5 (Supplementary file). Good agreement between the modelled and observed SOCs was observed. The RothC model could generally simulate SOC changes at most grassland sites, as reflected in the statistical indexes calculated (Fig. 2). A linear relationship was found between the measured SOC and simulated values (R 2 = 0.71). Positive EF values are indicative of satisfactory model simulations. RMSE values of model performance were smaller than the standard deviation of observations, demonstrating that the model can simulate the measured values. The RothC model marginally overestimated the SOC (0.002% of ARE).
Our simulations, which considered land management and changes in climate parameters, suggested that SOC stocks continuously increased from 1983 up to 2019. Simulations show an average C stock change rate of 1.26 t C ha −1 year −1 for grassland soils during the period of 1983-2019, which represented a total accumulation of 44.17 t C ha −1 . (1) Regional-scale soil organic carbon dynamics in response to climate projections Table 2 shows the simulated change in SOC from 2020 to 2090 for different ranges of an initial SOC recorded in the region for the 3 climatic subregions identified. Grassland soils of Basque Country would have a total potential C accumulation of 1.62 Tg under RCP4.5 and 1.37 Tg under RCP8.5, and the difference was statistically significant (p < 0.001). The overall trend indicates faster rates of increased SOC stock in soils with a lower initial SOC (< 50 and 50-79 t C ha −1 ), with values higher than 0.56 t C ha −1 year −1 for subregion 1 (Table 2). In contrast, a negative impact on SOC was observed in the soils of subregions 2 and 3 with a high initial SOC stock (> 110 t C ha −1 ), with a maximum average SOC reduction of 0.68 t C ha −1 year −1 . The specific patterns of SOC change under the RCP4.5 and RCP8.5 scenarios are shown below. Figure 3 also shows that there is a range of variability in the cumulative SOC. The ranges vary depending on FYM or plant production inputs. The medians obtained under both scenarios show differences in the accumulation of SOC of 21.11 t C ha −1 during 2020-2090. This variability was 4.11% lower in the RCP8.5 scenario than that in the RCP4.5 scenario. The worst and best scenarios were found in response to plant production inputs, with accumulations of only 1.03 t C ha −1 for low plant input in RCP8.5 (second quartile) and 43.59 t C ha −1 for high plant input in RCP4.5 (third quartile) over 70 years. This result shows the higher range of variation in C plant inputs under RCP4.5.

Characteristics of soil organic carbon change under the RCP4.5 and RCP8.5 scenarios
The results showed a similar trend of SOC variation in the short term (2020-2040) under both the RCP4.5 and RCP8.5 scenarios. From that time, it can be observed that the SOC stock increases at a lower rate until a stable trend is reached in 2070. Under RCP4.5, 0.25 Tg more is estimated than under RCP8.5.
According to the simulations under RCP4.5, in general, grassland soils with a current SOC < 110 t C ha −1 would increase their SOC over the period 2020-2090 (Table 2). The rate of SOC stock would increase more under soils with lower SOC stocks, especially for grassland soil with a current SOC < 50 t C ha −1 , where the increase in SOC would reach, on average, 48.40 t C ha −1 for the period 2020-2090. Similarly, the largest SOC reduction (13.56 t C ha −1 on average) was found in soils with the highest initial SOC stocks (> 110 t C ha −1 ).
Among climatic subregions, coastal sites (subregions 2 and 3) would show higher losses of soil C than the southern part of the region (Table 2). In fact, SOC losses are minimal in subregion 1 under RCP4.5 and only in soils with high concentrations of carbon in 2020. However, SOC losses in subregion 3 occurred at a lower initial SOC (80-109 t C ha −1 ), reaching a total SOC reduction of 46.88 t C ha −1 in the soils with the highest initial SOC. Figure 4 shows the simulated percentage of annual change in SOC compared with 2020 values under both climatic scenarios (RCP4.5 and RCP8.5). Overall, simulations show an increase in SOC over the period 2020-2090 under RCP4.5, which included 89.06% of the total grassland soils studied, but the rate of SOC increase would slow over time. In 2090, most soils will show an increase between 0 and 0.8%. The areas of higher SOC change are estimated in subregion 1, where the soils start with a lower SOC stock than those for subregions 2 and 3. In subregion 1, most soils show an SOC increase from 0.39 to > 0.94%. Under RCP8.5, an SOC increase would take place at a rate lower than 0.4% in comparison to RCP4.5. Most of the SOC reductions are predicted in subregions 2 and 3, where losses of SOC would begin sooner than in subregion 1.

Contribution to the 4 per 1000 objective
If we express the values from 1983 to 2019 in terms of the Initiative, it can be said that grassland soils in the area of study incorporate an average of 2.7%, more than the Initiative objective.
Overall, our results suggest that the Initiative would be feasible as long as SOC stocks are maintained or even increased under both scenarios in soils with an initial SOC stock of < 80 t C ha −1 . For most topsoil data, the increase in SOC stock varied within the subregion and by initial SOC stock. Thus, during the whole period, soils with a low initial SOC (< 50 t C ha −1 ) presented the highest SOC variation per year, with an average increase of 1-1.65%. Soils with an initial SOC of up to 50-80 t C ha −1 increased 0.4% per year in subregions 1 and 2, while those in subregion 3 presented an increase of 0.23% per year under both climate change scenarios.
Therefore, at least a 0.4% increase in SOC per year would be reached from 2020 to 2090 in grasslands located in climatic subregions 1 and 2 with initial SOC stocks up to 80 t C ha −1 under both scenarios (Table 2). For subregion 3, this goal would be met for grasslands containing up to 50 t C ha −1 .
In the short term (2020-2040), the objective of increasing SOC by 0.4% would be met for 83.51% of the area occupied by grasslands under both climate scenarios (Fig. 4). These percentages would be reduced to 61.48% and 47.60% and to 57.14% and 38.78% under the RCP4.5 and RCP8.5 scenarios in the medium and long term, respectively. Future data from soil monitoring programmes could validate these estimates.

Model performance
The first phase of this study was to calibrate and validate the RothC model with a historical dataset of grassland soil parameters. The site was selected as representative of the climatic subregions defined for the Basque Country region. The fields are located in an area that exhibits climate variability in the region and represents grassland management for forage production. Acceptable correlation parameters were found between most observed and simulated SOC values, suggesting that the model was able to capture temporal dynamics in SOC. This result could be attributed to the fact that model initialization was adjusted for each field, which is a favourable aspect to reduce the error, as reported by Smith et al. (1997a). The model fit the data with an R 2 of 0.71, which is comparable to similar studies in the literature (Cerri et al. 2007;Nakamura et al. 2011;Afzali et al. 2019).
The increase of 1.26 t C ha −1 year −1 in the 0-30 cm layer during 36 years of simulation agrees with results published by Conant et al. (2001), who reported a range of 0.11 to 3.04 t C ha −1 year −1 in grassland under different management practices. They observed that C accumulation rates were highest during the first 40 years after the application of treatments aimed at grassland improvement. Fornara et al. (2016) reported a soil C accumulation rate of 0.88 t C ha −1 year −1 in intensively managed grassland fields in Ireland during a 45-year-long grassland experiment, reflecting high application rates of cattle slurries that led to significant increases in soil C stocks.
The increase in SOC recorded can be attributed to repeated practices of C input for forage production, which is consequently modelled in our simulation setup. Nevertheless, this value should be interpreted with caution, as there are other factors affecting soil C dynamics, such as soil type, deep soil horizons, and microorganism interactions (Chenu et al. 2019).
In our study area, several good practices had taken place for years, such as organic fertilization, diversification of grass species, and use of legumes. In temperate grasslands, it has been reported that the intensity of management, number of cuttings, and nitrogen (N) fertilization strongly affect SOC storage (Soussana et al. 2006). Grassland fertilization has been used for a long time in our study area to increase forage production. Additionally, legume species can increase soil N, resulting in superior soil fertility and increasing aboveground and belowground production.
Animal manure is assumed to be more easily decomposed than plant material (Coleman and Jenkinson 1996). In contrast to inorganic fertilization, a portion of increased soil C is attributable directly to the addition of manure C to the soil, and the effect of this variable depends on the quantity and quality of the organic material applied (Bhogal et al. 2018). In this sense, characterizing the decomposability of Values on the X axis represent C inputs for plant and FYM used in the simulations. A value of 0.115-0.253 t C ha −1 was used as the mean plant production of the subregions. Lower and higher values were set at 0.14-0.07 and 0.21-0.37 t C ha −1 respectively. For FYM, 5.3 t C ha −1 and a variation of 20% were used to cover a reasonable range of organic fertilizer applications the organic amendments applied would help improve model performance.
In addition, it is well recognized that clay offers physicochemical protection to SOC, allowing long residence times of organic matter (Hassink 1997;Six et al. 2002). Soils sampled at the field scale show a clay composition of approximately 25%, which would contribute to preserving the SOC that can be stored in association with this fraction.
Furthermore, even though our study site is located at the intersection of different climate subregions, weather data tended to be cooler than those at coastal sites, which may have contributed to maintaining and increasing SOC stocks, as observed in other studies (Minasny et al. 2017), together with continuous organic fertilizer application for forage production. Climate constitutes an important factor for SOC storage, and strong relationships between SOC and

Regional estimation
It is standard practice to upscale SOC models from the plot to regional level with the aim of predicting changes in SOC. Regional analysis involves subdividing the study area into several spatial units (i.e., grid cells) for which sets of driving variables (e.g., climate, soil characteristics, and land management) are derived and then supplied to the model. In this study, these grids were developed using a regional climate model at 1 × 1 km resolution. After establishing the baseline SOC level for 2020, a 70-year simulation was run for grassland soils classified depending on their SOC stock. Correlations found between the observed and simulated SOC can be used to estimate the SOC change in response to climate change (Falloon and Smith 2006;Ludwig et al. 2007). At the regional scale, most grassland soils fall within the range of 50-109 t C ha −1 , which coincides with the value given by Soussana et al. (2006) for land under permanent grassland in France with average stocks of nearly 70 t C ha −1 . Our results show that soils with the lowest initial SOC stock would present the highest SOC increase under all subregions and climate scenarios, reaching an annual increase ranging from 1 to 1.65%.
This value considerably decreased for soils containing 50-79 t C ha −1 , with an annual increase from 0.23 to 0.94%. According to Minasny et al. (2017), who analysed the accumulation potential from 20 regions around the world, there is a tendency of a higher C accumulation potential in soils with a low initial SOC stock. Conant et al. (2001) also established that large changes were slightly more likely for soils with low initial C and slightly less likely for soils with high C. As expected, low-C soils may have a higher capacity for C storage. Wan et al. (2011) also predicted the highest decrease in net SOC, possibly because of the relatively higher original SOC stock than that in other areas. In our study, future estimates in relation to the historical period showed a 2.95 times slower accumulation rate under the RCP4.5 and RCP8.5 scenarios for soils in the 50-79 t C ha −1 initial SOC range.
Climate change was expected to have an impact on SOC in the long term, while changes in the short term would more likely be driven by land management practices and land use change (Pedroli and Meiner 2017). Our results show a higher response to SOC change when high C plant inputs were simulated under RCP4.5. Our explanation is that differences in yield production among subregions in response to climate change would be lower under scenario RCP8.5 than those under RCP4.5, thus leading to lower variability among plant C input parameters, which would reduce SOC change variability. Soils in the initial SOC range of 110-139 t C ha −1 showed a slight increase under both scenarios in subregions 1 and 2 (RCP4.5), while those in subregions 3 and 2 (RCP8.5) presented even negative SOC change. In fact, we found the lowest SOC stock change in soils located in subregion 3, with no difference between climate scenarios. These results agree with those reported for Irish grassland fields by Byrne and Kiely (2009). These authors found that the SOC stocks decrease with increasing rainfall. In our study, the rainfall records for subregion 3 differed significantly compared to those for subregions 1 and 2. The rate modifier for soil moisture depends on precipitation, clay content, and evaporation rate (Coleman and Jenkinson 1996). In addition, subregion 1, with current cooler average temperatures and lower predicted temperature increases in RCP4.5 and RCP8.5, showed higher SOC increases than subregions 2 and 3, even for soils with high organic matter due to lower SOM mineralization rates. In this sense, previous studies have found that SOC stocks are highest under cool conditions even at the subregional scale (Wiesmeier et al. 2019).
Under temperate climate, SOC stock change takes place slowly, and periods of years and even decades can be required to allow for quantification of rates of change (Poulton et al. 2018). Our data from regional simulations show that there is substantial variability in the rates of C storage over time, but it appears that increased soil C may persist for grassland soil for 70 years. Previous long-term studies have reported that changes in SOC may take as long as 60 years for degraded soils (Potter et al. 1999). Nevertheless, it is difficult to determine how long soil C will continue to increase with improved management due to the interaction of other factors, such as climate and soil type. Nevertheless, it is well established that SOC does not increase indefinitely but moves towards equilibrium (Smith 2004). In our study, equilibrium will be attained sooner for region 3 soils under scenario RCP8.5 with organic matter content in the range of 80-110 t C ha −1 . As indicated by Stockmann et al. (2013), increases in SOC following a change in management are initially rapid and then slow, reaching a new quasi-equilibrium at some point in the future. Grassland sites with an initial SOC ≥ 80 t C ha −1 had been subjected to improved management for years; thus, the change in SOC stock was lower or even negative for soils with a high initial SOC, with soil C loss up to 0.42% year −1 under RCP8.5. Regarding clay content, soils with an initial SOC higher than 80 t C ha −1 would apparently not be affected by clay content, as SOC was high enough to saturate the clay particles. These observations are in agreement with those of Hassink (1997), who also found that the amount of organic compounds that can be bound to clay particles is limited.
Thus, C loss would be attributed to the combination of initial SOC stock, clay, and climate change scenarios or climate characteristics at the subregional scale. In this sense, Bellamy et al. (2005) reported C losses in soils with a high initial SOC up to 2% per year. They found that these C losses took place irrespective of land use, suggesting a link to climate change. The effects of temperature will interact in a complex way with changes in soil moisture brought about by changing rainfall and evapotranspiration patterns. In wet anoxic soils, increased temperature will lead to increased evaporation, which in turn will accelerate the rate of decomposition. The increasing rainfall and soil water content could have promoted higher decomposition rates (Sommers et al. 1981), counteracting the effect of higher temperatures. In our study, the climate subregion had a greater influence than the climate change scenario due to the larger differences in rainfall found between subregions.
Predictions under the RCP4.5 scenario presented either a significantly higher increase or lower reduction in SOC than under the RCP8.5 scenario in subregions 1 and 2. Given that the RCP4.5 scenario was associated with lower changes in rainfall and temperature, modelled changes in SOC using the RCP4.5 scenario seemed to favour SOC stock increases or maintenance depending on the initial SOC stock. This observation is in agreement with the fact that SOC stocks generally increase as the mean annual temperature decreases provided there is no limitation in soil water content. Kirschbaum (2000) concluded that global warming is likely to reduce SOC by stimulating rates of decomposition. In fact, whether SOC stocks increase or decrease under climate change depends upon which process dominates at a given location. Increasing temperatures can accelerate decomposition and promote the loss of soil C but can also be slowed by lower soil moisture.

Contribution to the 4 per 1000 objective
This work showed that grassland management carried out in the study area for 36 years led to SOC increases of 27 per 1000 at 30 cm depth. Poulton et al. (2018) found that for permanent pastures under organic fertilization application every 4 years, which is considerably less frequent than ours, reported an annual SOC increase from − 0.1 to 7.4 per 1000 during the first 20 years.
We evaluated to what extent grassland soils in the region would comply with the 70-year simulation with the Initiative for three regional climatic subregions under the same usual grassland management to assess the effect of climate change. The rates of SOC increments averaged 0.42% for all grassland soils evaluated. According to Minasny et al. (2017), although C accumulation rates varied between countries and climatic conditions, there was a trend according to the type of management, and SOC stock rates were 0.5 t C ha −1 in the case of organic amendments. At the end of the simulation period, in 44-54% of cases, SOC increases occurred by more than 0.4% per year. SOC increases occurred mostly at a rate greater than 1.48% per year at 0-30 cm depth during the first 20 years of the simulation period under both scenarios. Small differences were found between scenarios RCP4.5 and RCP8.5 within each period, with increases exceeding the 4 per 1000 per year objective continuing in some areas after 70 years. Simulations show that SOC changes would take place slowly, as also observed by Poulton et al. (2018) in longterm agricultural field experiments. Poulton et al. (2018) reported a rapid increase in SOC during the first 20 years following the application of organic fertilizer to soil with initially poor in organic matter (< 30 t C ha −1 ). The rate of SOC increase after applying FYM for many years was predicted to be higher (> 0.85% per year) in subregion 1 than in subregions 2 and 3, where most soils showed an SOC increase between 0 and 0.57% per year. As mentioned above, although we cannot provide a complete explanation for this finding, it could be due to differences in climate parameters at the subregional scale. Minasny et al. (2017) also observed that high C accumulation rates (up to 10 per 1000) can be achieved for soils with a low initial SOC stock (less than 30 t C ha −1 in the topsoil) during the first 20 years after the implementation of best management practices. In contrast, and as mentioned above, areas that have reached equilibrium will not be able to further increase their accumulation.

Overcoming uncertainties and limitations of modelling
C inputs from plants play an important role in the RothC model due to the sensitivity of the model to that parameter and the fact that this is the most uncertain parameter (Poeplau 2016). In this study, a calibration of the C input from the plant was carried out during the initialization of the model. According to FAO (2019), this is a useful method for improving estimates of C input, so for future simulations, the amount of C input from plants is less uncertain. The extent of climate change impacts has shown complexity and differences over space and time (Ainsworth et al. 2019;Liu et al. 2021). For this reason, we undertook a variability analysis to ascertain the model response to a range of C inputs from plants and FYMs.
The opposite is true for climate data. Due to the simplicity of obtaining climate data from a site, this is a parameter with low uncertainty, while the sensitivity of models to climate is very high. In our study, the precision of C stock estimates relied on using specific climate data for the number of points in the region at 1 × 1 km resolution, although climate variability at the local scale is frequently small (Wiesmeier et al. 2019). This approach allowed us to identify and characterize climate parameter differences at the subregional scale that had an important effect on rate modifiers of the model, such as precipitation.
Initial SOC stock is another sensitive parameter (FAO 2019). We used initial SOC values that are based on numerous measures covering the whole region simulations. In this way, we considerably reduced the uncertainties associated with initial SOC stocks. At the regional scale, each spatial unit is treated independently, and there are no interacting processes between them. Thus, exchanges between grid cells such as water, and lateral transport of C in sediments and soil (i.e., erosion) are not included and can contribute to additional uncertainty in C stock changes at a specific location, as observed by Paustian et al. (2016).
In this study, we combined measurements and modelling approaches, and we reduced the typical uncertainties of regional modelling by calibrating the model for a study area that captures climate variability and C inputs from plants and FYMs.

Conclusions
We conclude that RothC could adequately simulate changes in C in grassland soils in Basque Country. However, more soil sampling sites for grasslands should be established to increase the capacity to calibrate and validate SOC models for Basque Country conditions. Our findings show that for future predictions, soil C change in Basque Country is highly dependent on the initial C stock, and C input from plants and climate subregions, with a higher risk of C loss in soils rich in organic matter under higher precipitation regime subregions. The overall trends showed that future climate change might decrease SOC in soils containing 110-139 t C ha −1 if no appropriate practices are taken to maintain the SOC stock.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.