Determination of minimum number of samples allowing to detect long term soil organic carbon changes in Mediterranean arable lands using paired-sites

Background Legacy data are unique occasions for estimating soil organic carbon (SOC) concentration changes and spatial variability, but their use can pose limitations due to the sampling schemes adopted and improvements may be needed in the analysis methodologies. When SOC changes is estimated with legacy data, the use of soil samples collected in different plots (i.e., non-aligned data) may lead to biased results. In the present work, N=302 georeferenced soil samples were selected from a regional (Sicily, south of Italy) soil database. An operational sampling approach was developed to spot SOC concentration changes from 1994 to 2017 in the same plots at the 0-30 cm soil depth and tested. Results The measurements were conducted after computing the minimum number of samples needed to have a reliable estimate of SOC variation after 23 years. By applying an effect size based methodology, 30 out of 302 sites were resampled in 2017 to achieve a power of 80%, and an a=0.05. A Wilcoxon test applied to the variation of SOC from 1994 to 2017 suggested that there was not a statistical difference in SOC concentration after 23 years (Z = -0.556; 2-tailed asymptotic signicance = 0.578). In particular, only 40% of resampled sites showed a higher SOC concentration than in 2017.


Introduction
Soil organic carbon (SOC) is a main contributor to fertility in agricultural soils, which also guarantee water accumulation and biodiversity (Lal, 2008). Baseline SOC estimates and maps are generally built on legacy data (Odeh et al., 2012), whereas any new collection of soil samples in the same legacy locations are often scarce. Up to date SOC data and assessments are on the global agenda (Vermeulen et al., 2019) and are necessary to evaluate many ecosystem characteristics such as resilience, productivity, ability of soil to provide a wide range of ecosystem services (Williams et al., 2020), and to gain precious insight into policy measures for soil preservation .
Prior to a SOC assessment, a sampling campaign would be needed, and the number of samples would affect obtained results. Many sampling size determination strategies have been proposed in the last decades to spot SOC changes (Biswas and Zhang, 2018;Köhl et al., 2011;Lark, 2009;Parras-Alcántara et al., 2015;Wadoux and Brus, 2020). Most of these have been suggested for new data acquisition and to overcome the problems of di cult terrain. There is now a strong and growing need to utilize legacy soil data sources for monitoring SOC changes (Boubehziz et al., 2020;Chen et al., 2019;Francaviglia et al., 2012). Regions that have soil monitoring networks need periodic recollection of soil samples to evaluate changes over time. Such resampling could be minimized to contain costs, but should be large enough to produce reliable estimates.
The need for long-term temporal paired sites is essential when aiming to depict SOC changes (Bellamy et al., 2005;Zhang and McGrath, 2004). Few countries have such monitoring schemes. The European Union (EU) started around 10 years ago with a Pan-European monitoring network to improve sustainable farming solutions and monitor soil pollution (Ballabio et al., 2014;Lugato et al., 2014).
Soil and SOC conservation is critically important in semi-arid areas, where deserti cation risk is increasing . For these areas, a recent increase in the output of literature regarding SOC accounting and spatial modelling (Smith et al., 2005), legacy databases (Schillaci et al. 2018), and digital soil mapping (Minasny and McBratney, 2016) has been noted.
Legacy soil data and soil maps could be integrated into a uni ed database. This would provide special insight into hardto-sample areas, past and present trends, and insight into the application of proper modelling procedures.
Recently, the development of digital soil mapping, or pedometrics, and the presence of an ample archive of historical soil data has allowed for the assessment of country scale SOC patterns with relatively high accuracy (Krol, 2008;. However, models can amplify uncertainty when the assessment is based on multiple predictors (Grunwald et al., 2018;Shi et al., 2018).
In some areas of the world, a lack of recent SOC measurements is prompting a rediscovery of legacy data which is in the process of being fully integrated into mapping methods at an operational level (Arrouays et al., 2017). SOC distribution is determined by multiple factors (Lacoste et al., 2014), the importance of which vary mainly with bioclimatic conditions. It is therefore hard to delineate general functions that explain the world SOC distribution using only geographical position, although a general inverse correlation was found with average annual air temperature on a regional scale between 52° N and 40° S and a direct correlation beyond this region .
Land use and land use change is also a main driver of SOC stocks, although mechanisms of SOC dynamics seem to be independent of the ecosystem type or land use (Giannetta et al., 2018). Guo and Gifford (2002) showed that around 50% of SOC is gained in the transition from cropland to secondary vegetation communities with a meta-analytic approach, and recent papers con rmed such a trend (Wei et al., 2014;Zhou et al., 2015). Sommer and Bossio (2014) hypothesized that SOC sequestration in arable land can show a 0.012-0.027% annual increase in the rst two decades after the establishment of SOC preservation practices, after which a saturation occurs and the increase ceases. Following the same hypotheses, Zomer et al. (2017) presented a global assessment of cropland SOC under the aforementioned scenarios and found that the potential SOC sequestration in cropland is below 53% of the 4p1000 target (Minasny et al., 2017).
High spatial variability and temporal trends induced spatial modellers to design reliable sampling strategies (Biswas and Zhang, 2018;Schmidt et al., 2008;Shiwen et al., 2017) and develop e cient methods to compare intra-eld and inter-eld variations (Pezzuolo et al., 2017) with similar agro-ecological conditions over the course of two decades (Kühnel et al., 2019). Application of this technique was carried out to determine the effects of sampling density on interpolation accuracy  and uncertainty assessment (Szatmári et al., 2018;Veronesi and Schillaci, 2019).
Cropland covers 12.6% of the world's surface (FAOSTAT data, accessed 2019). Cropland SOC has been mapped on a global scale using the WoSiS database (Batjes et al., 2017;Hengl et al., 2014;Zomer et al., 2017) and SOC maps were obtained by applying Generalized Additive Models (GAM) and machine learning methods.
At the European level, cropland plays an important role. Due to the large area covered, cropland acts as a potential carbon (C) sink. If considering a biomass return of up to 45 Mg C per year in raw materials, the biological potential of cropland for C storage is on the order of 90-120 Mg C per year (Freibauer et al., 2004;Smith, 2004). In particular, Smith (2004) demonstrated that models of SOC changes should be used with statistical power analyses for planning sample design to determine density and time of sampling during experiments.
Little is known about long-term SOC changes in Mediterranean semi-arid arable lands, which are frequently dominated by winter-growing species (mostly cereals and legumes) in rotation with fallow periods characterized by various crop residue management practices. Field crop production in these areas can cause SOC depletion (Rodríguez Martín et al., 2019) and soil loss by erosion, especially when conventional tillage (e.g., ploughing) is continuously applied. Conversely, no-till has been shown to be strongly bene cial compared to conventional tillage in semi-arid climates with an aridity index lower than 0.52 ± 0.03 (Gristina et al., 2018). Land management practices such as reduction of tillage intensity (Sperow, 2020), addition of manure and sowing cover crops could help to increase SOC contents and keep signi cant amounts of nutrients (Álvaro-Fuentes et al., 2012b;Barbera et al., 2012).
In addition, in Mediterranean areas, frequent re events can burn tons of biomass; this lowers the yearly C input derived from crop residues utilization while increasing SOC permanence, and affect the cycle of several nutrients including nitrogen (Knicker, 2007). In these areas, short-term SOC changes due to management practices (and especially land use) can temporarily override background changes (Lou et al., 2011) since length of cultivation is a main driver of SOC variation (Wang et al., 2019). In particular, time-space expected variation in soil traits can be differentially modelled if a priori information is available. Such information can allow for the determination of a minimum sample size to test a hypothesis effectively. De ning the sample size and location is required to enhance the power analysis while reducing laboratory costs and maximizing the accuracy of the representation (Confalonieri et al., 2009).
This experiment aimed at verifying whether or not a legacy estimation of SOC changes (1994-2008 model results) from non-paired data (Schillaci et al., 2017) matches the SOC variation measured in paired soil samples after 23 years . Thirty temporal paired sites from Sicily (South of Italy) and under continuous crop cover were resampled  and included in the present study. The analysis focused on arable land as it represents the main land use in the study area. The land cover of these sites was veri ed using historical remote sensing imagery to con rm that each site was continuously cultivated during the intervening period. Minimum sample size was determined and locations were randomly selected. Topsoil SOC contents were determined using the same laboratory method as in 1994 and 2008 (Walkley and Black, 1934).

Study area
Sicily (25,286 km 2 ) is a semi-arid to arid region of Italy. Sicily is characterised by prolonged droughts from mid or late spring to early or mid-fall, with the addition of high energy storms in fall and winter. Rainfed arable land selected using the CORINE code 211 is the most common land cover class in the area under study with roughly 300,000 ha yearly under cultivation with durum wheat (Istat, 2013).
Rainfed arable lands represented the target land cover in the study area as they represent 60% of the surface in this region.
Thus, it is a primary candidate for C sequestration and mitigation of the anthropogenic impact on the landscape. The land is predominantly under private ownership and the average farm size is around 6 ha, in general family-run businesses (all farm types), with approximately 10% of foreign labour (Istat, 2013).

Sampling campaign
In the rainfed arable lands of Sicily, the seedbed is generally prepared by soil ploughing during late summer and one or two harrowing in early fall. The amount of nitrogen (N) applied in non-legume eld crops is usually between 80-100 kg N ha − 1 year − 1 and durum wheat yield is between 2-4 Mg ha − 1 (with a harvest index ranging from 45 to 55%, and therefore with a similar straw yield). Land cover data was derived from the Geographic information system of Sicily (SITR; http://www.sitr.regione.sicilia.it/). Availability of a legacy databases and determination of the minimum number of samples and their locations for SOC determination.
In this work, data obtained from a sampling campaign carried out in 1993-94 (SMSLD) were used. In particular, data about the rainfed arable land (CORINE code 211) were identi ed using the optimization procedure shown in Schillaci et al. (2019). Brie y, the whole SMSLD consists of 6674 checked samples. Samples from the CORINE Land use 2 (agricultural lands) were 5471 from 2886 locations, with 2, 2, and 3 layers sampled at the 1st, 2nd, and 3rd quartile, respectively. Within these 2886 locations, samples from the sole CORINE Land use 2.1 (arable land) were 2162 from 880 locations, with 2, 2, and 3 layers sampled at the 1st, 2nd, and 3rd quartile, respectively.
A power analysis (Lakens, 2013) was used to nd the minimum number of samples needed to determine the SOC change with time (from 1994 to 2017). The expected change with time was derived from a modelled SOC variation (1994 to 2008) at the regional level, as observed in Schillaci et al. (2019). In particular, samples from the same soil layer, land use and sub-area were used (Fig. 1).
To de ne the effect size, means and standard deviation of the oldest SOC survey (1993) and the hypothesized change after 15 years were used. The effect size was computed at different degree of con dence (0.1, 0.05, 0.01) using the G-Power software (Faul et al., 2007). The estimated change comes from the difference between the average topsoil SOC concentration (0-30 cm) measured in 1994 and the estimated values at the same location in 2008 (predicted) at a 1-km spatial scale (Schillaci et al., 2017). This allows for six sets of data (Table 1) used as raw or log 10 transformed data, none of which was normally distributed after a Shapiro-Wilk test of normality (Table 2). Table 1 Datasets use in the present experiment. Log 10 of each of these databases were also computed.

SOC17-SOC94
30 Differences between data measured in 2017 in the same locations of data selected in 1994 The expected difference in fourteen years from 1994 (measured data in LEG-SOC94) to 2008 (estimated data in EST-SOC08) was subjected to the effect size computation (i.e., Cohen's d) by means of the Wilcoxon signed-rank test (Wilcoxon, 1945) in IBM SPSS software. Standard deviation of legacy data was computed. Cohen's d provided the minimum number of sites to be sampled in the 2017 survey on the original locations of the 1994 (LEG-SOC94) to ascertain if a change in SOC content occurred. When this number was achieved, the sampling consisted in the collection of topsoil (0-30 cm) in 30 randomly selected sites for which coordinates were reported with up to six digit precision (thus with an error < 30 m; Schillaci et al., 2019).
A two-tailed paired Wilcoxon-test was performed because, despite previous ndings (Schillaci et al., 2017), there was uncertainty regarding the pattern of differences between SOC2017 and SOC94. Such uncertainty derived from the internal variability of the model used in Schillaci et al. (2017) suffered by the lack of paired locations. Indeed, dynamics of the difference from the data in 1994 (measured in LEG-SOC94) to 2008 (estimated in EST-SOC08) was on average positive for each land use, but with a high variation. Determination of sampling sites was carried out via random selection. This random selection was performed within the set of locations determined in 1994 (SOC94).

Calculations of the legacy topsoil SOC
Original legacy SOC data were collected at various depths. To compare SOC data sampled in the 2017 at a uniform depth to the ones sampled in 1993, the Hobley and Wilson (2016) method was used to uniform the SOC concentration of the former sampling campaign. Such a method is based on an exponential generalized function, as follows: where d is the depth (expressed in meters), SOC 0 is the concentration of SOC at the soil surface (%), and k is the depletion constant (m − 1 ). The Hobley and Wilson (2016) method rst ts this model and nds an optimal k and SOC 0 for each location, then computes SOC at any depth (d). Note that in the original work (Hobley and Wilson, 2016), Eq. (1) also contains an additional term (i.e., SOC∞) modelling the concentration of residual SOC to a soil depth tending to in nity. In this study, SOC ∞ was assumed to be null.
In this work, the Hobley and Wilson (2016)

Historical land use and soil type in resampled sites
Web Mapping Services (WMS) from the regional geodata service (http://www.sitr.regione.sicilia.it/), consisting of aerial photograph surveys, were used to check the historical land cover. Sites were checked to ascertain that these sites maintained the same land use during the intervening period. To ensure this, aerial photographs of at least 3 surveys carried out from 1994 to 2017 (Fig. 2.) were consulted. To avoid sites being temporarily converted to short-term grassland or bare soil or abandoned, local farmers were also interviewed.
Soil samples from the selected sites, rainfed arable land (N = 302) derived from the 1994 sampling campaign, and the resampled soils (N = 30) had a consistent distribution among the main soil systems according with the Soil map of Italy (Costantini et al., 2014) (Supplemental materials Fig. 1). Around 80% of sites were Vertic soils or Cambisols or Regosols. All of these soils frequently showed calcic or calcaric subgroups.
2.5 Soil organic carbon determination SOC17 measurements were carried out using the Walkley-Black method (Walkley and Black, 1934). The SOC values so obtained were within the range of LEG-SOC94. In LEG-SOC94, soil data were georeferenced, and this has allowed for the resampling at exact ( eld site) location.

Computation of SOC differences between 1994 and 2017 by means of ad hoc sampling
Since the residues from the means of the logarithm 10 of SOC94 and SOC17 were normally distributed (Table 2), the Log 10 -SOC variation from 1994 to 2017 was tested for difference from zero by means of a paired t-test. This does not assume that observations within each group are normal, but only that their residuals are normal (Hsu and Lachenbruch, 2014;Zar, 1999) (Table 2). Such assumptions were only partly met, since the Log distribution were normal, whereas the differences between log were not normal but its skewness was comprised between − 1.96 and + 1.96. To overcome such problems, a bootstrap ANOVA was performed (Krishnamoorthy et al., 2007). Changes in cropland SOC content have previously been accounted for with the same statistical approach (Wang et al., 2018) where no subsampling was done to assess the SOC change.

Power analysis
Farmer interviews con rmed that the land use of the sampling sites have been continuously cultivated during the intervening period. The crop rotation was generally made with durum wheat (Triticum durum Desf.) followed by broad beans (Vicia faba L.), or other pulses alternated with fallow. The sampling campaign LEG-SOC94 involved 302 sites, with a SOC concentration in the CORINE land use 2.1 (arable) of 1.01 ± 0.59% (mean ± s.d.; Table 3) after normalization to 0.3 m depth. This value was signi cantly lower than the mean SOC in the 30 sites selected in 1994 (i.e., the SOC94, 1.31%; Table 3), but very close to its median (1.05%). Mean (predicted) SOC of the EST-SOC08 database was 1.38 ± 0.39% (mean ± s.d., with a median of 1.25%). SOC was expected to vary in the cropland between the 1994 and the 2008 values when using original data or Log 10 data (Table 4). Such variation was due to an increase in SOC in 75.8% of the sites (2008). The calculated effect size of the SOC was 0.54 for the original data and 0.69 for the Log 10 data (Table 5). According to Cohen (1988), such effect sizes correspond to 'medium to high' effect, which needed a minimum sample size ranging from 15 to 45 samples to be able to detect a SOC variation. Given this effect size and the power chosen for the Wilcoxon test, which is by default set to 80%, and a signi cance level of 5%, the calculated sample size required would be 30 samples. These were identi ed in LEG-SOC94 and collected in their respective paired locations SOC17. In 2017, only 12 sites showing a SOC concentration higher than 1994 were found, so SOC variation from 1994 to 2017 depended more on the SOC difference within each pair of samples than on the % of samples in 2017 having a SOC higher than in 1994.  *for logs, mean of log represent the arithmetic mean of the logarithms 10 of the SOC values expressed in %; ** such an estimation is from the 302 sites in which a measure was available 1994 and applying the 1-km resolution estimation process used in the 2008 described in (Schillaci et al., 2017). Test statistics based on positive ranks (Wilcoxon) Z based on negative ranks −9.119 −10.243 Asymptotic signi cance (2-tailed) < 0.001 < 0.001 Table 5 Output parameters of the a priori power analysis computation process at varying the α (0.10; 0.05; or 0.01). The process was carried out through the G-Power software with the Wilcoxon test for non-normal distributed datasets (Faul et al., 2007). Input data were from 302 measured samples of SOC in 1994 and modelled SOC in 2008, each of which expressed as either raw or log10 data.  (Table 3). SOC94 had a mean of 1.31 ± 0.96% (mean ± s.d.), whereas in the LEG-SOC94 sites, mean was 1.01 ± 0.59% (Fig. 3). The SOC17 mean was on average signi cantly lower than EST-SOC08, but with similar distribution properties to both the SOC94 and LEG-SOC94. Also, the SOC17 mean was slightly lower (−6.1% relative change, − 0.08% absolute change) than the mean of SOC94, but not than the EST-SOC08. The accuracy of the tted depth functions was expressed with the RMSE and MAE resulting in 0.36 and 0.23 respectively. Transformation to logarithms improved the distribution properties in term of skewness and kurtosis, especially for SOC17 (Fig. 4). According to the Wilcoxon test, the change in SOC between SOC94 and SOC17 did not appear different from zero (Table 6 and Fig. 5). The t-test applied to the Log 10 of SOC17 and SOC94, in which residues were normally distributed, provided similar results compared to the Wilcoxon test carried out on the same data (SOC difference = + 0.03300 ± 0.27659, mean ± s.d.; C.I 95% : −0.07028 to + 0.13628; t = + 0.653; 2-tails signi cance = 0.52). The bootstrap ANOVA also provided consistent results compared to the Wilcoxon test (Table 7), but notably, the difference in the log-SOC between the sampling campaigns had 95% con dence intervals marginally overlapping zero. This work was aimed at assessing the reliability of estimations of topsoil SOC changes with non-paired data using timepaired sampling. Compared to a previous work on temporal variation in the same area (Schillaci et al., 2017), various improvements in methodology were applied. These improvements consisted of the use of temporal paired sites instead of estimates based on different locations. Locations were further checked for land use continuity. The previous analyses based on modelled data over a 15-year span (from 1994 to 2008) predicted a mean relative increase of around the 21% of SOC content in arable lands, but such an increase was affected by a strong variability when survey data was reviewed. By using a monitoring network spanning 30 years, Gubler et al. (2019) found that SOC dynamic is more determined by a change in land use than other predictors in a colder climate (Switzerland).
According to a hypothetical linear growth, the increase expected from 1994 to 2017, which was the target period of the present study, was predicted to be a relative SOC increase of around 30%. In this study, no evidence of this magnitude of increase was found, and no such difference was seen in neither the original data nor in the log transformed data. When using a similar amount of data from meadows, Gubler et al. (2019) found that the minimum detectable change in 10 to 100 years (at a power of the 80%, such as in the present study) spanned from approximately 2 to 6%, respectively. Bellamy et al. (2005) showed that variation in time strongly depended on the initial SOC content, e.g., sites with low SOC had more opportunity to increase their SOC than sites with high SOC. However, both these latter studies were conducted in wetter climates than that of the present study. Using data from 20 regions in the world, Minasny et al. (2017) showed a tendency of a higher C sequestration potential on croplands with low initial SOC stock (≤ 30 t C / ha at 0-30 cm) compared to grasslands, which already have a high initial SOC stock, although a general decrease of the stock rate with time was also reported.
In addition, most of the samples in the present study were derived from the thermo-Mediterranean bioclimatic area (Rivas-Martínez et al., 2011), so these results can also re ect the latency of SOC variation in these areas, which even under soil abandonment showed limited increases in SOC when compared to other bioclimates in the same region (Novara et al., 2017).
The lack of an increase found from 1994 to 2017 when compared to the estimated increase from 1994 to 2008 was partly expected given that a SOC mean of the SOC94 extracted from the LEG-SOC94 was higher than that of the complete LEG-SOC94 dataset. In few exceptions (3 samples) there was an increase in SOC concentration that can also be due to soil deposition following erosion from the sites at higher altitudes (Lizaga et al., 2019;Navas et al., 2017;Panagos et al., 2015), changes in tillage depth or soil compaction (Álvaro-Fuentes et al., 2012a). In addition, in this area soil respiration due to increasing temperatures (Viola et al., 2014) may have offset any potential increase in SOC, and the 23-year time span may have not been su cient to detect a SOC change, as pointed by Saby et al. (2008). Soil management with ploughing and the increasing mean temperature (Viola et al., 2014) are not conducive to SOC accumulation (Goidts and van Wesemael, 2007;Kämpf et al., 2016). In particular, Goidts and van Wesemael (2007) showed that ploughing may override the increase in SOC over time. Also, an increase in SOC in ploughed soils is hard to achieve unless high quantities of organic residues and N are provided (Mazzoncini et al., 2011). These two latter conditions are very limited in Sicily due to the low crop yield, the low amount of residues returned to the soil and scarce fertilization. Other authors showed that slight increases in temperature and water availability ratio may contribute in reducing SOC (Davidson and Janssens, 2006). In the present work, this may have occurred and may be what can be seen when comparing the expected change from 1994 to 2008 to the measured change from 1994 to 2017. In addition, the ratio between water availability and temperature may be increasing in the area under study, despite the fact that no direct report is available (Viola et al., 2014).
The lack of increase in SOC concentration found here (from 1994 to 2017), compared to those estimated in the period from 1994 to 2008, could also be due to a SOC reduction from 2008 to 2017 which cannot be excluded using the present data. Such a reduction from 2008 to 2017 might possibly be connected with the decoupling of EU Common Agricultural Policy (CAP) payments regarding agriculture in (Regulation EEC 1782/2003. Before this year, wheat was the continuous primary crop for arable lands and after this point crop rotation with legumes or fallow land was encouraged by new regulations. Notably, continuous wheat has been shown to favour SOC accumulation when compared to a high percentage of wheat-legume rotations or wheat-fallow rotations, in the same or similar environments (Barbera et al., 2012;López-Bellido et al., 2010;West and Post, 2002). Lastly, differences between the previous estimation (1994 to 2008) and the present measurements (1994 to 2017) can also depend on the differences between direct and indirect measurement (Smith et al., 2019), or transient changes in cultivation history that may have in uenced the previous estimates (Wang et al., 2019), the latter of which was discarded here by an ad-hoc sampling in continuously ploughed soils with eld crops.
Variability and con dence intervals regarding SOC suggest that the estimated change in the previous work (Schillaci et al., 2017) could have been affected by outliers or errors in the measurements related to analytical methods used (Gessesse and Khamzina, 2018), especially when used for highly alkaline soils with scarce SOC (Bisutti et al., 2004;Walkley, 1935), such as in the present study. These issues may have produced an apparent pattern of SOC accumulation from 1994 to 2008 that was not detected between 1994 and 2017. Problems in the estimate due to the sampling strategy (Brus, 2019;Wadoux and Brus, 2020) were excluded, since the selection method used here provided a dataset with similar statistical properties to the original SOC distribution, and thus allowed for a maximum reduction in sampling sites.
Mitigation strategies and international projects have contributed to a debate regarding SOC sequestration. These traits are increasingly being taken into account in EU subsidies to the agricultural sector (Ciliberti and Frascarelli, 2015;European Comission, 2017).
There were different ongoing discussions following the paper by Sommer and Bossio (2014) and more recently after the Soil 4x1000 initiative (Minasny et al., 2017), all of which are above all pivotal to a reliable estimation of SOC dynamics. Zomer et al. (2017) modelled the estimated increase in SOC concentration and stocks at a global scale in croplands and found that an average increase of approximately + 26%, which is similar to the estimated value found from 1994 to 2008.
However, the Zomer et al. (2017) estimation may not be suitable for small scale assessment and mapping and the use of legacy information can be crucial to con rm these trends.

Conclusion
The estimated SOC change in arable lands occurred from 1994 to 2008 in the studied area was not con rmed by a direct measurement using paired sites in the 1994-2017 timespan. This result has a direct implication for the SOC monitoring network in the mid-term (e.g., 15-25 years). Results also suggest that measurements of soil properties based on legacy data need to be supported by reliable information on the land use, land use changes and soil management practices. Both these aspects were indirectly taken into account here by choosing sites with no change in the land use or soil management. These variables could be used to correct the effect of other predictors (e.g., rainfall, slope). Only when many limitations are found, there would be the need for large scale sampling campaigns.
The debate regarding the structuring of subsidies should take into account the current information to assure on the one hand, their competitiveness, and on the other hand, their environmental sustainability and ability to provide ecosystem services. Further works should aim at: i) increasing the number of sites to be resampled which would be derived from both the collection in 1994 and 2008, ii) gathering information regarding land use dynamics in paired sites, and iii) quantifying the effect of soil erosion in the ow of SOC within and among catchments. Lastly, direct information regarding the variation in SOC concentration and stock in topsoil and subsoil in the target areas are urgently needed to drive policy making. The discrepancies between the present data compared to the previously published results may depend on various factors, including:

Declarations
Availability of data and materials The legacy data that support the ndings of this study were made available from [the Regional Bureau for Agriculture, Rural Development and Mediterranean Fishery, the Department of Agriculture, ARTA. Palermo]. The new data samples taken in 2017 can be obtained upon request. The legacy data are subjected to licensing; in this study they were used under an agreement.

Competing interests
The Authors declare that they have no competing interests

Funding
The work is partially founded with the H2020 program Landsupport, ID: 774234 Authors   Study area, original soil data (dots) and original compared to new soil samples (histograms, background mask the Arable land cover from CORINE 2000, code 211 (yellow). Blue bars for SOC94, red bars for SOC17.

Figure 4
Distribution of the 30 samples of SOC 94 (blue bars) and SOC17 (red bars) and relative probability distribution function. Upper panels are for raw data, lower panels for Logs Figure 5