Introduction of Cardoon ( Cynara cardunculus L.) in a Rainfed Rotation to Improve Soil Organic Carbon Stock in Marginal Lands

: The production of a biomass as a feedstock for bioreﬁnery is gaining attention in many agricultural areas. The adoption of bioreﬁnery crops (i.e., perennial cardoon) can represent an interesting option for farmers and can contribute to increase soil organic carbon stock (SOCS). The study aimed to assess the potential e ﬀ ect on long-term SOCS change by the introduction of cardoon in a Mediterranean marginal area (Sassari, Italy). To this end, three process-oriented models, namely the Intergovernmental Panel on Climate Change (IPCC) guidelines for national greenhouse gas inventories (Tier 2), a humus-balance model (SOMBIT) and Rothamsted carbon model (RothC), were used to compare two scenarios over 20 years. The traditional cropping system’s faba bean–durum wheat biennial rotation was compared with the same scenario alternating seven years of cardoon cultivation. The model’s calibration was performed using climate, soil and crop data measured in three cardoon trials between 2011 and 2019. SOMBIT and Roth C models showed the best values of model performance metrics. By the insertion of cardoon, IPCC tool, SOMBIT and RothC models predicted an average annual SOCS increase, whereas, in the baseline scenario, the models predicted a steady state or a slight SOCS decrease. This increase can be attributed to a higher input of above- and belowground plant residues and a lower number of bare soil days (41 vs. 146 days year − 1 ).


Introduction
In the future decades, the transition from a fossil-based to a biomass-based economy is crucial for economic and environmental sustainability; therefore, Europe aims at achieving an optimal and efficient multipurpose exploitation of renewable biomasses as food, added valued chemicals and energy [1 -3]. As well, the impact in biomass production for the bioeconomy is crucial for environmental sustainability [4,5]. Land degradation caused by human activities negatively impacts the well-being of at least 3.2 billion people, costing more than 10% of annual global gross products through the loss of biodiversity and ecosystem services [6]. Soil degradation, as part of land degradation, includes the loss of soil through erosion and the depletion of soil organic matter (SOM) [3,7]; it has been estimated that nearly two billion ha (23%) of the world's used lands-agricultural lands, permanent pastures, forests and woodlands-have been degraded since the mid-twentieth century [8], causing important implications, particularly in the Mediterranean area [9]. Soil organic Carbon (SOC), conventionally assumed to be 58% of the SOM [10], is beneficial for soil health for many reasons, including improved Table 1. Thirty-year long-term monthly climate data  in Porto Torres Municipality (Sassari, Sardinia, Italy).

Farm Management and Soil Sampling of Cardoon
At sowing, the seedbed was prepared between October and November ( Figure 1) with ploughing at 20 cm, followed by harrowing at 10-15 cm. During the first year, mineral fertilization with P 2 O 5 and N was added, and weeds were chemically controlled. CAR was harvested every year when the moisture content of the seeds and stalks was about 8% and 15%, respectively. Thus, harvest was performed in late August by combine harvester, separating seeds and aboveground biomass in bales after raking. CAR was fertilized yearly starting from the second year, after the harvest (usually in October or November) and before it "covered" the field.
2018. Dominant reference soil groups of the World Reference Bases for Soil Resources classification in that area are Luvisols [46], and the climate is typical Mediterranean, with hot and dry summers. The 30-year average temperature, rainfall and evapotranspiration are 16.77 °C, 567.99 mm and 1128.98 mm, respectively (Table 1).

Farm Management and Soil Sampling of Cardoon
At sowing, the seedbed was prepared between October and November ( Figure 1) with ploughing at 20 cm, followed by harrowing at 10-15 cm. During the first year, mineral fertilization with P2O5 and N was added, and weeds were chemically controlled. CAR was harvested every year when the moisture content of the seeds and stalks was about 8% and 15%, respectively. Thus, harvest was performed in late August by combine harvester, separating seeds and aboveground biomass in bales after raking. CAR was fertilized yearly starting from the second year, after the harvest (usually in October or November) and before it "covered" the field.
In the three farms, soils were sampled at 0-30-cm depths in different growing seasons ( Figure  1). The samples were air-dried, powdered and sieved at 2 mm prior to being analyzed for the main physicochemical soil properties. Bulk density (BD; kg m −3 ) was estimated by using a pedotransfer function specific for Mediterranean regions from clay, sand and SOC content (% w/w), following Pellegrini et al. [ (1)  Farms Dec-11 Aug-12 Oct-12 Dec-12 Aug-13 Dec-13 Aug-14 Aug-15 Aug-16 May-17 Aug-17 Aug-18 Aug-18 Jun-19 Legend In the three farms, soils were sampled at 0-30-cm depths in different growing seasons (Figure 1). The samples were air-dried, powdered and sieved at 2 mm prior to being analyzed for the main physicochemical soil properties. Bulk density (BD; kg m −3 ) was estimated by using a pedotransfer function specific for Mediterranean regions from clay, sand and SOC content (% w/w), following Pellegrini et al. [47]:

Assessment of Plant Residue Inputs
To run the SOMBIT and RothC models, annual aboveground plant residues (R ag ) were measured for CAR and estimated for DW and FB. Belowground (R bg ) of all crops were estimated from aboveground biomasses. Regarding arable crops, aboveground FB was derived from Razza et al. [24], while for the DW by statistical data [44], using a harvest index (HI) of 0.40 [48]. For both crops, R bg were estimated from the R ag according to previous studies [40,49,50]. Weeds were estimated as 7% of the total aboveground biomass [40]. Regarding CAR, the main agricultural operations were recorded yearly. Three different crop-growing periods were considered to calculate the plant residues. The first period refers to the first growth year (YR1), in which grain and biomass yields are generally low; the second period goes from the second to the sixth years (YRS2-6), in which the grain and biomass yields are estimated as the mean of five years of production, and the third period is the explant year (YR7). CAR straw yields were measured, whereas R bg were estimated. The R ag of CAR was about 30% of the aboveground biomass left on the ground after harvest. Moreover, further residues came from basal leaves (8-10% of foliar biomass), stem and threshing residues of leaves produced during the annual crop cycle. Weeds were present only in YR1 before the CAR canopy development and accounted for 6% of the total CAR aboveground biomass by experimental observation.
Total CAR belowground biomass was mainly allocated in the topsoil (70-80%) and estimated considering a 1:1 shoot (as leaves, stalks and heads) to root ratio [51]. R bg included: (i) the annual decomposed fractions of roots, except for the first year (YR1), when root renewal was negligible and for the explant year (YR7) when the total root biomass was accounted for and (ii) the root exudates following Di Bene et al. [38] and Pausch and Kuzyakov [52]. Moreover, it was estimated that the CAR root biomass totally renewed in three years [53]; consequently, one-third decomposed yearly. Table 2 shows formulas used to calculate plant residues (Mg ha −1 ) according to the procedure outlined above and following previous studies. Table 2. Formulas used to calculate plant residues (Mg ha −1 ) for faba bean (FB), durum wheat (DW) and cardoon (CAR) crops.

Residues * FB and DW Crops CAR Crop
* Calculated as the sum of the aboveground residues (R ag ), belowground residues (R bg ) and residues from weeds (R w ), where Y G is the grain yield and Y S is the straw yield (Mg ha −1 of dry matter); HI is the harvest index (Mg Mg −1 ); E is the ratio between the exudates and roots: 0.11 in FB, 0.09 in DW and 0.15 in CAR; YR referred to the CAR year: YR1: year of the first harvest, YR2-6: from second to sixth harvest and YR7: the explant year.
The SOMBIT model directly used dry matter plant residues as the input, whereas, to run RothC, plant residues were previously converted to C inputs. For FB and DW, the C content of both aboveand belowground was assumed to be 45% [40]. Measured CAR C content was 45% and 49% for aboveand belowground, respectively.

IPCC for Italy at the Farm Level
In 2006, the Intergovernmental Panel on Climate Change developed a guideline for estimating the SOCS over time [19] used by countries for reporting to the United Nations Framework Convention on Climate Change. The IPCC set up the methodology for estimating the SOCS after a 20-year steady-state Agronomy 2020, 10, 946 5 of 18 management by considering three different tiers of complexity. In this study, we used the tier 2 method, applying the SOCS increase as predicted by the model in redefined standard SOCS values following the tier 1 procedure but calibrating the model from the observed SOCS data. Therefore, the yearly changes in the SOCS were calculated by the difference between the final and initial predicted SOCS values. The input data required for this model are the mean annual temperature ( • C), potential evapotranspiration (mm), frost occurrence, soil type and land-use class. In addition, for croplands, tillage and organic input classes were also considered.
The model was run to match the values of observed SOCS in each farm in the first soil sampling year (i.e., first and second year after CAR planting for farms 1 and 2 and farm 3, respectively). Since in farms 1 and 2, the first soil analysis occurred in 2013, the model was run in inverse mode to calculate the initial SOCS and starting the scenario in 2012.

The SOMBIT Model
The SOMBIT model was implemented in the BIT3G project [54] for estimating the site-specific SOM stock dynamics in cropping systems carried out to produce raw materials for biorefinery. It belongs to the humus-balancing family model referred to the SOM turnover as a function of ecological site properties and input amount and quality [55].
SOMBIT model has a year-step and uses primary or published data to estimate (i) the increase of SOCS due to the amount and type of biomass supply, the latter through the isohumic coefficient (k 1 ) relying on biomass compositions in soluble fraction, hemicellulose, lignin and cutin, crude fiber and ash contents [55][56][57][58], and (ii) the decrease of SOCS through the mineralization coefficient (k 2 ) influenced by pedoclimatic conditions and farming practices [59][60][61]. SOMBIT implements in Hénin-Dupuis model the estimation of the management coefficient using fuzzy logic to reach continuous instead of discrete values. It has already been applied in annual [62], annual-perennial [24] and perennial [63] agricultural systems.
The input data required to run the model are: mean annual temperature ( • C), bulk density (kg dm −3 ), coarse fragments (dm 3 dm −3 ), SOC (% w/w), clay and calcium carbonate contents (g kg −1 ), depth of the soil layer considered (cm), presence or absence of irrigation, organic matter supply (% w/w), tillage type, frequency (n year −1 ) and depth (cm). In addition, annual organic input supply is considered, i.e., amounts of above-and belowground residues and organic manure-if applied-(Mg ha −1 ) with respective isohumic coefficients (kg kg −1 ).
The k 1 attribution for the R ag and R bg of FB, DW and CAR was included in the parametrization. The model calibration was performed by attributing the initial SOCS equal to those determined with the first soil analysis performed for each farm, because the model cannot run in inverse mode.

RothC Model
The Rothamsted carbon model (version 26.3; RothC-26.3) is a process-based and multicompartmental model to simulate the long-term C turnover in the topsoil [25,64] for different vegetation types, including arable crops, grasslands and forests, on a monthly step base [65,66].
The model splits the SOC into five pools: one inert organic matter (IOM) pool-which does not participate in C turnover-and four active fractions, i.e., decomposable plant material (DPM), resistant plant material (RPM), microbial biomass (BIO) and humified organic matter (HUM). The decomposition rate is influenced by clay content, soil moisture, air temperature and plant cover.
The model input data are climate (monthly temperature, precipitation and evapotranspiration); soil (clay content and initial SOC S ) and land use and management information, such as soil cover (bare or covered), quantity and quality of plant residues returned to soil and organic exogenous inputs, if any.
To calibrate the model, monthly C inputs for FB-DW and CAR were estimated after Farina et al. [40] and Raccuia and Melilli [51], respectively, whereas the DPM/RPM ratio after Coleman and Jenkinson [25] and Jones et al. [67]. The model was first run iteratively in inverse mode to equilibrium (10,000 years) under constant environmental conditions and with the default pools and coefficient in order to generate the C input required to match the initial value of SOCS [68]. Then, RothC was run for calibration using the data of farm 1, with the C input estimated as in Section 2.3 from 2012 to 2018 (7 years) and by adjusting the DPM/RPM ratios up to fit the measured SOCS variation in the soil [69]. Subsequently, the 2012-2018 simulations of farms 2 and 3 were used to model the validation.

Scenarios Configuration
After the parametrization, calibration and validation processes, the three models were run for a 20-year period (2012-2031) using climate, soil and management data as the inputs to assess the potential SOCS change after the introduction of CAR in a typical Mediterranean industrial rainfed cropping system.
Two different scenarios, both conceivable, were hypothesized as follows: • Scenario S0 represented by 2-year rainfed FB-DW rotation (FB-DW) commonly cultivated in the study area under conventional tillage (annual ploughing). • Scenario S1 represented by the introduction of CAR cultivation in the conventional 2-year FB-DW rotation with this succession: CAR (7 years)-FB-DW-CAR (7 years)-FB-DW-CAR (2 years). Thus, the S1 represents a biorefinery-feed-food cropping system.

•
In both scenarios, the referred year is the year in which the crop was harvested. To explore the maximum potential for SOC sequestration, scenarios were run assuming all FB-DW and 30% CAR plant straws returned to the soil.
The potential SOC sequestration rate (SR) was used to compare the efficiency in the SOC sequestration after 20 years of S0 and S1 simulated managements. It was calculated as: where SOC final and SOC initial refer to the SOC values estimated at the end and beginning of the simulation, respectively, and I is the total C input implemented in the model during the 20-year period.
The SR values can be negative, null or positive. Value null indicates a steady state between C inputs and outputs and positive values show that C inputs are sequestered in soils, while negative values indicate that the C inputs did not offset SOC losses due to the mineralization process [70]. SR was calculated only for SOMBIT and RothC models, since the IPCC tool implements the amount of organic inputs only by qualitative classes.

Statistical Analysis
Main soil characteristics were analyzed by using one-way analysis of variance (ANOVA). If H0 data showed neither a normal distribution of error terms nor constant error variance, a nonparametric Kruskal-Wallis test was used. The analyses were performed using Statistica 7.0 (Statsoft©, Tulsa, OK, USA).
Model performances were assessed by comparing the simulated SOC dynamic against the observed measurements-excluding values utilized for calibration-to assess model performances using the following statistics: 1.
The root mean square error of prediction (RMSE) indicates the difference between observed and model-estimated values (range 0-∞). RMSE close to 0 shows a perfect agreement between observed and model-estimated values. The better is the simulation, and the lower is the value of RMSE [71].

2.
The relative error (RE) is the mean difference between measured and simulated data giving an indication of the bias in the simulation weighted as a percentage of the mean value of observed  [72].
The formulas of RMSE and RE are shown in Equations (3) and (4) below: where O is the mean values of observed data, respectively, O i and P i indicate the observed and predicted values, respectively, and n is the number of measurements.

Soil Characteristics
As a consequence of the farm size, ranging from 7 to 18 ha, the main physicochemical characteristics did not significantly differ among the farms. Soils were clay, with a common or frequent presence of coarse fragments, and medium bulk density. On average, the main differences were detected as a limestone content, i.e., scarcely calcareous in farm 3, moderately calcareous in farm 1 and very limestone in farm 2 (Table 3).

Plant and C Inputs
The total plant inputs (from R ag , R bg and R w ) and C input values of the soil are reported in Table 4 with the aim of defining the influence of different crop cultivations on SOCS changes. The total plant and C inputs varied throughout the cultivation years, depending on the amount of R ag and R bg returned to the soil.
In the three case studies, the total aboveground biomass of the CAR was quite stable, with a mean value of 11.6 Mg ha −1 y −1 . The highest and lowest productions were both observed in farm 1 during the YR2-6 and YR1 growth periods, respectively.
The CAR R bg was affected by the crop stand age and positively related to the root turnover, reaching the highest value in the explant year.
On average, the plant and C inputs observed in YR7 were around four times higher than the value obtained in YR1. The total amount of the plant and C inputs returned to the soil seven years after CAR cultivation were 94.8 and 46.3 Mg ha −1 , respectively.
The total plant and C inputs returned to the soil in the FB-DW rotation were estimated to be 11.0 and 5.0 Mg ha −1 , respectively (Table 4). a Year refers to the year of the CAR cultivation. YR1: first growth year, YR2-6: second growth period from the second to the sixth years, in which the grain and biomass yields are estimated as the mean of five years of production, and YR7: third period referring to the explant year. b Total estimated biomass observed at the harvest time as the dry matter, including the grain and straws. Belowground was estimated from the aboveground biomass; see text for details. c Residues refers to the biomass left in the field throughout the cultivation period, including basal leaves, threshing residues, roots at the explant and decomposed fraction of roots and root exudates following Table 2 formulas. d Total R is the sum of R ag + R w + R bg , and it was used as inputs to run the SOMBIT model. e C inputs are obtained by multiplying plant inputs by roots or stems' C contents, and they were used as inputs to run the RothC model. Figure 2 shows the observed SOCS values for the 2012-2018 period (grey dots), compared with the simulated SOCS trends for the same period and averaged for the three farms (mean ± standard error), reporting dots and lines in green, blue and red colors for IPCC tool (20-year step), SOMBIT (annual step) and RothC (monthly step), respectively. In the three farms, SOCS measured values changed over the seven years of CAR cultivation, reaching a mean value of 64.2 ± 2.1 Mg ha −1 at the explant year (YR7). After seven years of CAR, the final SOCS was 29% higher than the initial value.

SOC Stock Change and Model Performance
The validation of the three models was performed by comparing the simulated SOCS values with the observed data in the three farm trials. The IPCC tool (tier 2 assessment) was not able to satisfactorily reproduce the increasing trend of the observed SOCS values during the seven years of the CAR cycle (i.e., 1.3 Mg ha −1 predicted vs. 14.4 Mg ha −1 observed in 2018). Conversely, both SOMBIT and RothC models were able to adequately predict the SOCS increase in the same period (7.2 Mg ha −1 vs. 14.0 Mg ha −1 predicted with SOMBIT and RothC, respectively). In both the SOMBIT and RothC models, the highest increase was predicted at the explant of CAR (YR7) when all root biomass was completely accounted for R bg and recycled.
Models' performances are reported in Table 5. The best RSME value was observer under the RothC prediction. The lowest RE value was observed under the SOMBIT model.  The validation of the three models was performed by comparing the simulated SOCS values with the observed data in the three farm trials. The IPCC tool (tier 2 assessment) was not able to satisfactorily reproduce the increasing trend of the observed SOCS values during the seven years of the CAR cycle (i.e., 1.3 Mg ha −1 predicted vs. 14.4 Mg ha −1 observed in 2018). Conversely, both SOMBIT and RothC models were able to adequately predict the SOCS increase in the same period (7.2 Mg ha −1 vs. 14.0 Mg ha −1 predicted with SOMBIT and RothC, respectively). In both the SOMBIT and RothC models, the highest increase was predicted at the explant of CAR (YR7) when all root biomass was completely accounted for Rbg and recycled.
Models' performances are reported in Table 5. The best RSME value was observer under the RothC prediction. The lowest RE value was observed under the SOMBIT model.  Figure 3 shows the SOCS predictions by the three models over 20 years of simulation, for both S0 and S1. In all predictions, the SOCS changes under S0 showed a maintenance (IPCC tool) or slow decrease (SOMBIT and RothC), while in S1, a slight (IPCC tool) or a noteworthy increase (SOMBIT and RothC) in the SOCS values were observed.

SOC Stock Predictions by Modeling
In details, the SOCS changes predicted by the IPCC tool were unperceivable for both scenarios, with respect to the initial SOC value. After 20 years, the SOCS for the three farms did not change in S0 (49.3 ± 0.6 Mg C ha −1 ), while increased by 1.5 Mg C ha −1 in S1 (50.8 ± 0.6 Mg C ha −1 ).
The 20-year predictions by the SOMBIT model showed that the final SOCS under S0 reached 44.4 ± 1.0 Mg C ha −1 , decreasing by a mean value of 0.3 Mg C ha −1 year −1 . Contrariwise, in S1, the  Figure 3 shows the SOCS predictions by the three models over 20 years of simulation, for both S0 and S1. In all predictions, the SOCS changes under S0 showed a maintenance (IPCC tool) or slow decrease (SOMBIT and RothC), while in S1, a slight (IPCC tool) or a noteworthy increase (SOMBIT and RothC) in the SOCS values were observed.

SOC Stock Predictions by Modeling
In details, the SOCS changes predicted by the IPCC tool were unperceivable for both scenarios, with respect to the initial SOC value. After 20 years, the SOCS for the three farms did not change in S0 (49.3 ± 0.6 Mg C ha −1 ), while increased by 1.5 Mg C ha −1 in S1 (50.8 ± 0.6 Mg C ha −1 ).
The 20-year predictions by the SOMBIT model showed that the final SOCS under S0 reached 44.4 ± 1.0 Mg C ha −1 , decreasing by a mean value of 0.3 Mg C ha −1 year −1 . Contrariwise, in S1, the estimated increase in the SOCS was, on average, 0.9 ± 0.1 Mg C ha −1 per year, reaching, after 20 years, a mean SOC value of 67.9 ± 1.1 Mg C ha −1 .
Similarly to the SOMBIT model, the SOCS dynamic predicted by the RothC model showed that, in S0, the SOCS decreased annually by a mean value of 0.1 Mg C ha −1 , reaching, after 20 years, 48.1 ± 0.6 Mg C ha −1 SOCS values. Regarding S1, the predicted SOCS increased by 1.0 ± 0.1 Mg C ha −1 per year, reaching, after 20 years, the values of 69.5 ± 1.2 Mg C ha −1 . estimated increase in the SOCS was, on average, 0.9 ± 0.1 Mg C ha −1 per year, reaching, after 20 years, a mean SOC value of 67.9 ± 1.1 Mg C ha −1 .
Similarly to the SOMBIT model, the SOCS dynamic predicted by the RothC model showed that, in S0, the SOCS decreased annually by a mean value of 0.1 Mg C ha −1 , reaching, after 20 years, 48.1 ± 0.6 Mg C ha −1 SOCS values. Regarding S1, the predicted SOCS increased by 1.0 ± 0.1 Mg C ha −1 per year, reaching, after 20 years, the values of 69.5 ± 1.2 Mg C ha −1 .

Potential SOC Sequestration
As reported in Table 6, the SOMBIT simulations for the 20-year period showed that the final SOCS decreased in the baseline scenario (S0) in respect to the initial SOC value, while, in the alternative scenario (S1), the final SOCS raised by 18.1 Mg C ha −1 compared with the initial SOC value. The SR ranged from −0.1 to 0.2 under S0 and S1. Similarly, for the same period, RothC simulated a final SOCS that showed a slight decrease and a greater increase than the initial SOC value in S0 and S1, respectively. Table 6. Predicted SOC stock change, total carbon input (I) and sequestration rate (SR) in 20 years (Mg ha −1 ), including cardoon (S1) or not (S0).

Discussion
To the best of our knowledge, this paper is the first study in the Mediterranean Basin assessing the influence of CAR crop introduction in a typical industrial rainfed cropping system on SOCS change in a 20-year scenario perspective using a modeling approach.

Potential SOC Sequestration
As reported in Table 6, the SOMBIT simulations for the 20-year period showed that the final SOCS decreased in the baseline scenario (S0) in respect to the initial SOC value, while, in the alternative scenario (S1), the final SOCS raised by 18.1 Mg C ha −1 compared with the initial SOC value. The SR ranged from −0.1 to 0.2 under S0 and S1. Similarly, for the same period, RothC simulated a final SOCS that showed a slight decrease and a greater increase than the initial SOC value in S0 and S1, respectively. Table 6. Predicted SOC stock change, total carbon input (I) and sequestration rate (SR) in 20 years (Mg ha −1 ), including cardoon (S1) or not (S0).

Discussion
To the best of our knowledge, this paper is the first study in the Mediterranean Basin assessing the influence of CAR crop introduction in a typical industrial rainfed cropping system on SOCS change in a 20-year scenario perspective using a modeling approach.
Several studies have indicated that CAR-native to the Mediterranean Basin-is one of the most promising species for energy production in this area [73][74][75]. The crop characteristics that support these applications are the relatively low crop input and large biomass productivity, mainly of lignocellulosic compositions and high heating values [76].
Moreover, the introduction of these crops into the traditional agricultural system may represent an interesting alternative to increase the agroenvironmental cropping system's sustainability. This is because CAR directly enhances the soil fertility due to the reduction of the tsoil tillage frequency and low nutritional requirements [77] and reduces the soil degradation due to the protection of its dense canopy against the erosion [78,79].
CAR is a perennial herbaceous rhizomatous root species native to the Mediterranean Basin [14,80], which can grow up to 2-m-high and presents the leaves of the basal rosette petiolata very large and coriaceous [81]. It is an autumn-sown and summer-harvested crop, with an annual reproductive cycle, which is completed yearly by the end of spring [82].
In Mediterranean conditions, it is normally rainfed, and the biomass depends only on the amount of the rainwater accumulated during autumn, winter and early spring in the soil and extracted by its deep roots system [83]. Good performances occur with 450-mm water availability in the period ranging from shoot emissions to head differentiations [76]. The CAR biomass yields observed in the present study (Table 4) are comparable with the scientific literature available from similar studies carried out in other rainfed Mediterranean areas. In fact, the aboveground biomass production during the first growth year is usually low; then, the biomass production usually rises to the third and becomes more stable over the following years [84], decreasing by the seventh or subsequent years [77,85]. In this research, the average production in grain and straws was lower in YR1 and higher in YR2-6, while, in YR7, decreased slightly, hence corroborating the planning of a seven-year CAR cultivation also in the study area. Further, in a study carried out in the Basilicata Region (Southern Italy), Piscioneri et al. [86] observed a three-year average dry biomass yield of 12.5 Mg ha −1 (ranging from 10.0 to 15.0 Mg ha −1 ). Similarly, Ledda et al. [87] in Sardinia (Italy) observed a dry biomass yield of 10.2 Mg ha −1 (ranging from 6.8 to 11.6 Mg ha −1 ), while Francaviglia et al. [88] and Neri et al. [79] in a marginal hilly area of the Latium Region (Central Italy) reported a yield of 12.9 and 11.4 Mg ha −1 as the average, respectively. In Portugal, Gominho et al. [89] reported a total dry biomass yield of 9.7 Mg ha −1 (range 4.4-18.4 Mg ha −1 ), while in Northern Greece, Vasilakoglou and Dhima [90] observed dry biomass yields of 11.0 Mg ha −1 .
Since the cropping system's purpose aimed to maintain or increase the SOCS, in both scenarios, straws of FB and DW after the harvests were supposed to be left on the ground and incorporated by the next tillage. Indeed, in the study area, in 70% of the cropping systems, straw is left in the field [91], notably in marginal areas where yields are lower. CAR straws were harvested as biorefinery feedstock, although 30% of the aboveground biomass was left in the soil to increase sustainability, instead of 10%, as reported in previous trials [24].
The models required specific parametrizations. Following the IPCC classification (2006 to 2019), in the study area, (i) the climate fell in a "warm temperate, dry" region (comparing Table 1 data), (ii) soils were classified as "high-activity clay soils" (Table 3), (iii) land use was "long-term cultivated" and "cropland remaining as cropland". However, the SOCS predicted by the IPCC model changed in S0 and S1 due to variations in the soil management. The steady-state (S0) management fell in the "full tillage" class due to annual arable crops with generally full inversions, and organic inputs resulted "medium", according to the IPCC definition, because the annual crop residues were returned to the field, and N-fixing crop (FB) was present in the cropping system. Conversely, as a result of semi-perennial CAR cultivation (S1), tillage occurred only six times in the 20-years scenario (0.3 frequency), so the cropping system could be attributed to the "no till" IPCC class, and the estimated residues supplied in 20 years increased from 104 Mg ha −1 in S0 to 229 Mg ha −1 in S1, taking into account both agricultural cropping systems and the data reported in Table 3. For that reason, the organic inputs in S1 were attributable to the "high without manure" IPCC class.
The SOMBIT parametrization provided k 1 from the literature, except CAR R bg , that has been estimated. DW and FB R ag was 0.10 and 0.08, respectively, whereas R bg was 0.15 for both crops [92]. CAR R ag was 0.18, as experimentally calculated by Razza et al. [24], whereas R bg was estimated considering that roots of perennial plants have higher decay resistance than annual ones, due to the higher lignin contents [93]. In detail, k 1 was estimated as 0.27, increasing the k 1 R bg value applied for annual crops by the perennial and annual RPM ratio applied in RothC (1.80).
In RothC, the arable crop default value (1.44) of the DPM/RPM ratio was used for FB and DW residues, whereas, for CAR residues, 1 and 0.35 ratios were used for R ag and R bg , respectively. The DPM/RPM R ag value means 50% DPM and 50% RPM and is intermediate between the typical value for arable crops and scrubs (1.44 and 0.67, respectively). The DPM/RPM R bg value was lower following Kätterer et al. [93] and was assumed as the ratio of 0.26% DPM and 0.74% RPM.
The introduction of CAR (S1) led to a 75-kg C ha −1 annual increase applying the IPCC tool (2019) refining default values (Figure 3). This value seems rather underestimated compared with the observed data reported in Figure 2. Indeed, for both RMSE and RE statistics (Table 5), the IPCC tool showed the lowest performance, because it did not accurately predict either the SOC trend or its change. In the SOMBIT and RothC models, the simulated values remained within, or very close to, the standard error of the measured values. The high model performances in simulating the measured data in the studied area gave high reliability to the simulation in the perspective of the next 20 years. However, the IPCC tool results were preserved mainly to compare the predictions obtained with other models. Indeed, the increase in SOCS predicted by the SOMBIT and RothC models resulted higher than the IPCC model (Figure 3), because they considered even the annual residue quantities supplied and the specific soil features. This improvement in input quantification allowed to assess the different field conditions, reducing the uncertainty due to IPCC wide classes tailored for a national scale.
In S0, SOCS predicted by SOMBIT and RothC decreased, on average, by 267 ± 14 and 61 ± 0.3 kg C ha −1 per year. These findings are consistent with a decrease observed by Mazzoncini et al. [94] in a 28-years trial in similar conditions. Contrariwise, in S1, SOCS predicted by SOMBIT and RothC increased in average by 905 ± 110 and 1009 ± 53 kg C ha −1 per year, respectively. The higher increase occurred in the last year of CAR cultivation for both models when, at explant, all roots were left in the soil.
The long-term prediction highlighted the positive effect of CAR introduction in the rainfed two-year FB-DW rotation in terms of SOC sequestration. In this study, the SR was about 70% greater than that estimated in different cropping systems [67], resulting from the high C input from plant residues and a simultaneous reduction of C loss by no soil disturbance. Furthermore, the predicted increase is remarkable if related to the carbon footprint of CAR cultivated without organic input in the same area, as estimated by Razza et al. [24] and Cocco et al. [95], in 1.7 and 2.1 MgCO 2 eq ha −1 , respectively. Indeed, it is possible to convert the change of the SOCS in potential CO 2 emissions by the ratio of CO 2 and carbon molecular weights: 1 Mg ha −1 of SOCS corresponds to 3.66 Mg of CO 2 [96,97]; hence, the predicted SOCS increase in S1 can probably offset the carbon footprint of CAR. The 4p1000 initiative [98], despite the criticism on a global-scale application [99], could be an interesting target at the farm scale, applying the formula suggested by Priori et al. [97]. It allows to compare the 4% target SOCS increase with results of the model predictions. In the investigated fields, to reach the target, the SOCS after 20 years should be 53.9 ± 1.4 Mg C ha −1 , a result widely achieved in S1 following the RothC and SOMBIT models, while not in S0 (Figure 3).
An explanation to better understand the increase of SOCS, both observed and predicted in 20 years, might be attributed to the reduced bare soil introducing CAR (S1-72% bare soil days than S0), consequently reducing the SOCS mineralization.
Thus, the S1 scenario represents a possible option for farmers to adopt a more sustainable cropping system able to decrease C losses and to increase C inputs in comparison with the traditional systems.

Conclusions
The transition from a fossil-based to a bio-based economy in future decades is of crucial importance. The demand for more sustainable materials from agriculture to replace oil-derived supplies is ever-increasing due to ecological concerns [100]. The CAR, a perennial plant adapted to the Mediterranean low rainfall and hot, dry summer conditions, is a promising nonfood agricultural crop for marginal or set-aside lands suitable to feed a vegetable oil-based biorefinery. Within this research, the effects of the introduction of CAR on SOCS dynamics were assessed by different models, namely the IPCC tool, SOMBIT and RothC, whose outcomes were compared against the experimental farm data on SOCS in order to test the robustness of the provisional estimations. The observed results showed that the RothC and SOMBIT models provide the most accurate previsions. RothC slightly