Quantifying Soil Compaction in Persimmon Orchards Using ISUM (Improved Stock Unearthing Method) and Core Sampling Methods

: Agricultural activities induce micro-topographical changes, soil compaction and structural changes due to soil cultivation, which directly a ﬀ ect ecosystem services. However, little is known about how these soil structural changes occur during and after the planting of orchards, and which key factors and processes play a major role in soil compaction due to cultivation works. This study evaluates the improved stock unearthing method (ISUM) as a low-cost and precise alternative to the tedious and costly traditional core sampling method, to characterize the changes in soil compaction in a representative persimmon orchard in Eastern Spain. To achieve this goal, ﬁrstly, in the ﬁeld, undisturbed soil samples using metallic core rings (in January 2016 and 2019) were collected at di ﬀ erent soil depths between 45 paired-trees, and topographic variations were determined following the protocol established by ISUM (January 2019). Our results show that soil bulk density (Bd) increases with depth and in the inter-row area, due to the e ﬀ ect of tractor passes and human trampling. The bulk density values of the top surface layers (0–12 cm) showed the lowest soil accumulation, but the highest temporal and spatial variability. Soil consolidation within three years after planting as calculated using the core samples was 12 mm, whereas when calculated with ISUM, it was 14 mm. The quality of the results with ISUM was better than with the traditional core method, due to the higher amount of sampling points. The ISUM is a promising method to measure soil compaction, but it is restricted to the land where soil erosion does not take place, or where soil erosion is measured to establish a balance of soil redistribution. Another positive contribution of ISUM is that it requires 24 h of technician work to acquire the data, whereas the core method requires 272 h. Our research is the ﬁrst approach to use ISUM to quantify soil compaction and will contribute to applying innovative and low-cost monitoring methods to agricultural land and conserving ecosystem services.


Introduction
Soil degradation affects agricultural lands, which directly affects many ecosystems with clear benefits for human health such as agroecosystems. Healthy soil ecosystems can give integral benefits to the productivity of food and sources, but non-sustainable land management is drastically defaulting this. This is because cultivation results in higher soil and water losses [1,2], soil organic matter depletion [3,4], nutrients delivery [5], reduction in soil fertility [6], and soil compaction enhancement [7]. In agricultural lands, soil compaction is caused by the use of machinery, human trampling and the reduction of organic matter [8]. Then, soil structure degradation and a reduction in soil biota result in soil degradation and a loss of soil functions and services, which can be studied under different scientific points of view [9]. As a consequence, achieving sustainable development goals and land degradation neutrality objectives becomes even more challenging [10]. Stoessel et al. [11] affirmed that soil compaction is relevant to understanding the life cycle in the soil system. Furthermore, Mossadeghi-Björklund et al. [12] stated that the subsoil compaction by affected vehicular traffic determines the preferential flows in the clay soil.
Normally, soil compaction is a topic that has received minor considerations in studies related to ecosystem services when researching the sustainability of agricultural lands, and this is relevant in modern mechanized farming. Using modern machinery can cause problems due to soil compaction in groves and crops, where few types of research in agroecosystems have been conducted. Most of the research in orchards focused on the soil erosion problem and little research was conducted on soil compaction. The works of Van Dijck and Van Asch [13] on loamy soils in southern France contributed with key information from vineyards and orchards. Later, the works of Ferrero et al. [14] on the effect of the tractor and the one of Lipiec et al. [15] about the impact of the compaction on the water and heat transfer in the soil became the pioneers' works. The research developed in orchards is limited to the survey of Becerra et al. [16] in Almería, Spain, about the impact of tractor traffic in almond orchards, and Pižl [17], about the impact of compaction on earthworms in apple orchards.
During recent decades, scientists have tried to find which human-induced processes directly affect soil compaction. For instance, Liu et al. [18] studied the changes in soil bulk density under the wetting/drying cycles of irrigation regimes and found that soil compaction is dependent on the irrigation regimes. On the other hand, frequent traffic in agricultural land is one of the key issues responsible for soil compaction [19]. Long-term tillage and residue management also control soil compaction, as documented by Celik et al. [20], who determined soil bulk density using different tillage systems. There is a need to better understand the impacts of soil compaction on agricultural soils to control and mitigate its negative effects, and this has been a long-standing request of the scientific community and stakeholders [21][22][23].
The use of heavy machinery [24], cultivation methods, and human and animal (e.g., pastures) trampling [25] are the main factors of soil compaction, although human trampling used to be considered as a secondary one because of its lower impact [26][27][28]. Most of the research on soil compaction has taken place on heavily mechanized agriculture, where the machine weight per soil surface can be significant [29][30][31]. In the scientific literature, cereals and soya (Glycine max L.) crops are the ones that are the most investigated by the scientific community [11,32]. However, research on soil compaction in other agroecosystems such as citrus, olives or vineyards is still limited [33], mainly since mechanization in these crops has firstly taken place for flat areas such as terraces or large parcels, and then, for sloping areas.
The methods such as core sampling by a cylinder that are typically applied to quantify soil compaction are tedious, time-consuming, and expensive, requiring laboratory time and equipment [28,34,35]. There is a need to update the methods and to find a proper soil strategy to reduce the time and economic cost of the surveys [36], which is relevant for the development of science in developing countries. The improved stock unearthing method (ISUM) is a recently updated method that applies an easy strategy to determine the changes in the soil topography using the graft union of the plants as a passive bio-marker [37]. This method facilitates the estimation of the soil erosion rates by quantifying the differences in micro-topography between the initial planting and the time of measurement [38]. The survey of the soil topography at two different periods could also provide information about the soil compaction and mobilization during that period, to conserve ecosystem services.
In this research, the main goal is to assess the use of the ISUM to estimate soil compaction. The rationale of this investigation is to confirm if this new application is suitable in a persimmon (Diospyros kaki L.) orchard terrace. This experimental area was surveyed in 2016 after tillage and before plantation with contrasted high soil erosion rates during the first year of plantation [39]. We surveyed with the traditional core ring method to determine the soil bulk density (Bd) and perform a spatial correlation analysis among the obtained results. The comparison of both methods will shed light on the soil compaction and will elucidate whether ISUM can assist in the characterization of soil compaction in agricultural lands. Due to the direct connection between Bd and soil compaction, soil bulk density was used as a proxy measurement to assess soil compaction, as a tool to give another quality indicator for ecosystem service assessments.

Study Area
An agricultural terrace was selected as an experimental field in L'Alcúdia de Crespins municipality, Eastern Iberian Peninsula ( Figure 1). The mean annual precipitation is 532 mm and a mean annual temperature of 16.2 • C. Rainfall distribution is typically the Mediterranean, with three months of summer drought (24 mm of rainfall in July and August) and a peak of rainfall in October with 65 mm). Soil texture is clay-loam (21% clay, 39% silt, and 40% sand), soil organic matter low (1-1.2%), and soil pH reaches 8.1. We selected a traditional terraced field with a flat surface (504 m 2 -14.4 m × 35 m), surrounded by a concrete construction that prevents surface runoff and therefore the loss of water and sediments. The vegetation cover is low in autumn and winter, however, this is not the case in spring (Diplotaxis erucoides) and summer (Diplotaxis sp.) The cover was 80% in spring and summer and below 50% in autumn and winter.
information about the soil compaction and mobilization during that period, to conserve ecosystem services.
In this research, the main goal is to assess the use of the ISUM to estimate soil compaction. The rationale of this investigation is to confirm if this new application is suitable in a persimmon (Diospyros kaki L.) orchard terrace. This experimental area was surveyed in 2016 after tillage and before plantation with contrasted high soil erosion rates during the first year of plantation [39]. We surveyed with the traditional core ring method to determine the soil bulk density (Bd) and perform a spatial correlation analysis among the obtained results. The comparison of both methods will shed light on the soil compaction and will elucidate whether ISUM can assist in the characterization of soil compaction in agricultural lands. Due to the direct connection between Bd and soil compaction, soil bulk density was used as a proxy measurement to assess soil compaction, as a tool to give another quality indicator for ecosystem service assessments.

Study Area
An agricultural terrace was selected as an experimental field in L'Alcúdia de Crespins municipality, Eastern Iberian Peninsula ( Figure 1). The mean annual precipitation is 532 mm and a mean annual temperature of 16.2 °C. Rainfall distribution is typically the Mediterranean, with three months of summer drought (24 mm of rainfall in July and August) and a peak of rainfall in October with 65 mm). Soil texture is clay-loam (21% clay, 39% silt, and 40% sand), soil organic matter low (1-1.2%), and soil pH reaches 8.1. We selected a traditional terraced field with a flat surface (504 m 2 -14.4 m × 35 m), surrounded by a concrete construction that prevents surface runoff and therefore the loss of water and sediments. The vegetation cover is low in autumn and winter, however, this is not the case in spring (Diplotaxis erucoides) and summer (Diplotaxis sp.) The cover was 80% in spring and summer and below 50% in autumn and winter.
Surveys carried out during the last 3 years (2016 to 2019) after each rainfall event confirm that the surface wash has never exited from the terrace. The studied site was tilled in January 2016 to allow the planting of persimmon plants in March, and no other tillage has taken place since. Persimmon trees were planted at 4.8 m (inter-row) × 2.5 m (between trees in the same row), which is a usual pattern for persimmon orchards in this region. The plot was drip-irrigated with freshwater from the "Massís del Caroig" aquifer. Tractor passes took place for the application of chemicals (insecticides and fungicides), and the use of heavy machinery, human and animal trampling by the farmers took place during herbicide spraying and fruit harvesting.  Surveys carried out during the last 3 years (2016 to 2019) after each rainfall event confirm that the surface wash has never exited from the terrace. The studied site was tilled in January 2016 to allow the planting of persimmon plants in March, and no other tillage has taken place since. Persimmon trees were planted at 4.8 m (inter-row) × 2.5 m (between trees in the same row), which is a usual pattern for persimmon orchards in this region. The plot was drip-irrigated with freshwater from the "Massís del Caroig" aquifer. Tractor passes took place for the application of chemicals (insecticides and Agriculture 2020, 10, 266 4 of 18 fungicides), and the use of heavy machinery, human and animal trampling by the farmers took place during herbicide spraying and fruit harvesting.

Topographical Analysis Using ISUM
Three inter-row areas (four lines of trees) of 15 trees were selected. In each paired-tree, we marked the graft union following the procedures presented by Rodrigo-Comino and Cerdà [37], for the three measured rows, considering the updates based on using vine plants [40][41][42]. Following this methodology, a 1 mm nylon rope was used to connect trees 30 cm above graft unions of inter-row tree pairs. The 30 cm distance was chosen arbitrarily to facilitate easier sampling and avoid negative values due to terrain irregularities. Then, the vertical distance of the rope to the soil surface was measured at every 10 cm. We measured a total of 48 points perpendicular to each row, with 2160 points measured in the three rows used for the case study. During measurement analysis, the 30 cm offset, as well as the original distance of the graft union from the soil surface during planting, were subtracted from the measurements. We considered the original height of the graft union as a uniform for all trees and specified at 7 cm after personal communication with the farmer. Hence, we carried out the topographical measurements using ISUM in January 2019.

Soil Bulk Density (Bd)-Sampling and Calculations
We collected soil core samples twice (January 2016 and 2019) to measure Bd with 6 cm depth (100 cm 3 ) metalcore rings at 0-6, 6-12, 12-18 and 18-24 cm. We used stainless steel metal rings with 6 cm diameter and 1 mm in metal width. Those rings are also known as Kopecky rings and are widely used by the scientific community. Sampling beyond 24 cm was not considered since soil preparation tillage does not exceed this depth. Due to time and labor force constraints, the number of sampling points between each inter-row was 5 instead of the 48 sampled for ISUM; 75 in total. Nevertheless, for each sampling point, 4 soil cores were extracted at different soil depths (one core per depth, in total 4 samples), instead of a single surface measurement of ISUM. Therefore, in total, we collected 225 soil cores, to assess the spatial distribution of bulk density. Figure 2 shows the sampling positions for ISUM and Bd at the study site. After the calculation of dry soil weight (at 105 • C) and soil volume for each sample in the laboratory, Bd was determined. Then, the mean, median, maximum, minimum, standard deviation and coefficient of variation were calculated.

Topographical Analysis Using ISUM
Three inter-row areas (four lines of trees) of 15 trees were selected. In each paired-tree, we marked the graft union following the procedures presented by Rodrigo-Comino and Cerdà [37], for the three measured rows, considering the updates based on using vine plants [40][41][42]. Following this methodology, a 1 mm nylon rope was used to connect trees 30 cm above graft unions of inter-row tree pairs. The 30 cm distance was chosen arbitrarily to facilitate easier sampling and avoid negative values due to terrain irregularities. Then, the vertical distance of the rope to the soil surface was measured at every 10 cm. We measured a total of 48 points perpendicular to each row, with 2160 points measured in the three rows used for the case study. During measurement analysis, the 30 cm offset, as well as the original distance of the graft union from the soil surface during planting, were subtracted from the measurements. We considered the original height of the graft union as a uniform for all trees and specified at 7 cm after personal communication with the farmer. Hence, we carried out the topographical measurements using ISUM in January 2019.

Soil Bulk Density (Bd)-Sampling and Calculations
We collected soil core samples twice (January 2016 and 2019) to measure Bd with 6 cm depth (100 cm 3 ) metalcore rings at 0-6, 6-12, 12-18 and 18-24 cm. We used stainless steel metal rings with 6 cm diameter and 1 mm in metal width. Those rings are also known as Kopecky rings and are widely used by the scientific community. Sampling beyond 24 cm was not considered since soil preparation tillage does not exceed this depth. Due to time and labor force constraints, the number of sampling points between each inter-row was 5 instead of the 48 sampled for ISUM; 75 in total. Nevertheless, for each sampling point, 4 soil cores were extracted at different soil depths (one core per depth, in total 4 samples), instead of a single surface measurement of ISUM. Therefore, in total, we collected 225 soil cores, to assess the spatial distribution of bulk density. Figure 2 shows the sampling positions for ISUM and Bd at the study site. After the calculation of dry soil weight (at 105 °C) and soil volume for each sample in the laboratory, Bd was determined. Then, the mean, median, maximum, minimum, standard deviation and coefficient of variation were calculated.

Correlation between Bulk Density Soil Map and ISUM Map
Descriptive statistics were calculated (mean, median, maximum, minimum, standard deviation and coefficient of variation) and used for further analysis. Graphs of Bd data, which were measured using the core ring method at different soil depths, were depicted. Different methods of interpolation were used to prepare raster maps of collected data. The best method for interpolation was selected according to the root mean square error (RMSE) and the R 2 correlation coefficient. Nine methods were applied for interpolation, to prepare raster maps in ArcGIS 10.3 software (ESRI, Redlands, USA). These methods included ordinary kriging exponential (OK-XP), ordinary kriging Gaussian (OK-GA), were ordinary kriging spherical (OK-SP), inverse distance weighting (IDW), empirical Bayesian kriging (EBK), spline with tension (ST), multi-quadric (M-Q), inverse multi-quadric (IM-Q) and thin-plate spline (TPS). Empirical Bayesian kriging (EBK) was chosen as the best interpolation method based on the highest value of R 2 and lowest value of RMSE (0.8534 and 15.165, respectively).
After mapping Bd and micro-topography (topsoil level), a spatial correlation analysis was carried out among the raster maps using the band collection statistics tool in the ArcGIS 10.3 package.
The correlation coefficient Corr ij shows the relationship between the two spatial datasets i and j, and can be calculated following [43]: where δ i and δ j are the standard deviations of each dataset and Cov ij is the covariance matrix between datasets i and j defined as: where n is the number of cells in each dataset, Z k represents a particular cell value of dataset i or j, and µ is the mean value of each dataset. Correlation values range from +1 to −1, with positive values denoting a direct relationship between the two datasets, negative values denoting inverse correlation and zero denoting dataset independence.

Calculation of Soil Lowering Using Bulk Density Data Set and ISUM
When compaction takes place, the soil shrinks, and then Bd increases [44]. The following equation was used to calculate the lowering of soil by bulk density.
where Ls (mm) is the calculated amount of lowering soil by bulk density data, Hc (mm) is the height of the core cylinder, and Bd t1 (g cm −3 ) and Bd t2 (g cm −3 ) are the bulk density of the first and last measurement, respectively. Soil compaction was calculated using Equation (3) for all of the collected soil core data and each pixel of the soil density map, which was prepared from the "sum of all soil depths". Finally, the mean of each raster map obtained by ISUM was compared to the mean of the raster map of lowering soil obtained from bulk density data. Figure 3 shows the range of data distribution for four soil depths: 0-6, 6-12, 12-18 and 18-24 cm. The values of soil bulk densities are dependent on soil depth and compaction for three years ( Figure 3 and Table 1). We can observe that when soil depth increases, Bd increases. The lowest soil density was found from 0 to 6 cm depth. Differences in Bd indicate that Bd has increased so that the average for Agriculture 2020, 10, 266 6 of 18 the four measured depths of 0-6, 6-12, 12-18 and 18-24 in 2016 were 0.99, 1.05, 1.11 and 1.14 g cm −3 , and three years later, in 2019, these were 1.07, 1.11, 1.14 and 1.16 g cm −3 , respectively. The values of maximum Bd were 1.03, 1.09, 1.2 and 1.24 g cm −3 in 2016 and 1.14, 1.17, 1.2 and 1.23 g cm −3 in 2019, respectively, for 0-6, 6-12, 12-18 and 18-24 depth. The maps of soil density interpolation at different depths show that, after three years, the soil has shrunk (Figures 4-6). The data showed that the Bd increased mainly in the upper two sampling layers (0-6 and 6-12 cm), while deeper soil layers (12-18 and 18-24 cm) did not show clear changes, or even showed slight decreases in bulk density. The differences in Bd in the first and second layers during the study period are 0.075 and 0.055 g cm −3 , respectively.       Table 2 shows the descriptive statistics of the database surveyed using ISUM and soil lowering data obtained using Bd values. Regarding ISUM, the range of data changes is 453, since the minimum and maximum values of soil surface level are −339 and +114 mm, respectively ( Table 2). The average value reached −18 mm, indicating that the soil surface level has lowered by 6 mm per year on average.   Table 2 shows the descriptive statistics of the database surveyed using ISUM and soil lowering data obtained using Bd values. Regarding ISUM, the range of data changes is 453, since the minimum and maximum values of soil surface level are −339 and +114 mm, respectively ( Table 2). The average value reached −18 mm, indicating that the soil surface level has lowered by 6 mm per year on average. The ISUM map generated using the interpolation method is presented in Figure 7. It is worthy to highlight that the variation of soil lowering data obtained with ISUM, because of the larger number of measurements, is higher at shorter distances. However, we observe that it can also suppose a higher accuracy in presenting micro-topographical changes. According to this interpolation, the 2019 inter-row surface that coincides with the tractor traffic area registered a lower soil surface than the tree row soil surface. The ISUM map generated using the interpolation method is presented in Figure 7. It is worthy to highlight that the variation of soil lowering data obtained with ISUM, because of the larger number of measurements, is higher at shorter distances. However, we observe that it can also suppose a higher accuracy in presenting micro-topographical changes. According to this interpolation, the 2019 inter-row surface that coincides with the tractor traffic area registered a lower soil surface than the tree row soil surface.

Spatial Correlation Analysis between Bd and ISUM Maps
A correlation analysis between raster maps (ISUM and Bd), using band collection statistics, was performed, and the results are shown in Table 3. The relationships between micro-topography

Spatial Correlation Analysis between Bd and ISUM Maps
A correlation analysis between raster maps (ISUM and Bd), using band collection statistics, was performed, and the results are shown in Table 3. The relationships between micro-topography obtained using topsoil data were calculated, with 12 maps prepared for different depths in the surveyed years. Table 3. Spatial correlation between raster maps obtained using ISUM and Bd.

Comparison of Soil Lowering Calculated Using Bulk Density and ISUM Data Set
To measure soil compaction using the Bd data set, we applied Equation (3) to all obtained bulk density pixels. Table 4 shows the spatial descriptive statistics of soil lowering data and maps obtained using ISUM and Bd. By comparing the means, it is shown that there is a slight difference between the soil lowering obtained by ISUM and Bd. So, by comparing the raster maps, this difference is 2 mm, but for individual sampled points is 7 mm. Figure 7 shows the maps which were created using the bulk density data and ISUM. The mean value of the pixels for soil lowering using bulk density is −12 mm and for lowering soil map by ISUM is −14 mm. According to Table 4, the arithmetic mean for 225 soil density samples is −11 mm and for 2160 ISUM measurements is −18 mm. The raster maps bring a more accurate calculation of the soil lowering, due to the interpolation and calculation of the soil topography.

ISUM as an Efficient Method to Survey Soil Compaction for Ecological Services
The use of ISUM to calculate the soil compaction after planting will bring a reduction in time invested, and an improvement in the quality of data generated, as more sampling points can be taken. The use of maps adds a higher quality to the survey and allows for a comparison between core sampling and ISUM. Moreover, the core method is time-consuming and tedious [39]. In our experiment, the ring sampling was done in 2016 and 2019, and this involved the investment of four working days for three technicians during each sampling campaign; a total of 192 h. On average, each working hour allowed us to collect 9.375 ring samples. The laboratory processing (dry at 105 • C) and weighing were costed at ten labor days for a technician (900 samples); a total of 80 h. In total, 272 h were invested to characterize the Bd change in the persimmon plantation from 2016 to 2019. On the other hand, using ISUM, the measurement in 2019 required two technicians during a whole labor day (8 h), which is a total of 16 h. In total, the ISUM method required 24 h of field and computer work. We can classify ISUM as a low-cost field monitoring methodology, as it is one order of magnitude less time consuming than the traditional core sampling method.

Bd Temporal Changes
Our measurements in 2016 and 2019 allowed us to quantify the temporal changes in Bd. We found that, after planting the persimmons, the soil was compacted as a consequence of human trampling, machinery passes and the tillage needed before planting, as other researchers have also found in different crops [45,46]. Bd measurements carried out with a metallic ring showed compaction in three years, from 1.07 g cm −3 (average) to 1.12 g cm −3 (average). This means a soil lowering of the soil in 12 mm in three years. Soil compaction after preparation for planting in agricultural land or forest has been in different ecosystems [47,48], although, in recent years, the use of trees to restore degraded and compacted land has shown that trees aid recovery of soil properties and associate ecosystem service functions [49,50]. However, at the persimmon orchard studied in the Eastern Iberian Peninsula, we found that the soil is compacted after planting (4 mm of soil lowering per year-raster map), which is due to agricultural management (tractors and human trampling) during and after planting. Similar changes were found in vineyards, where high erosion rates and soil compaction resulted in environmental concern and threats for the agroecosystem services [51].
The usual method for measuring soil density is the core cylinder method, based on the use of metal rings (100 cm 3 ). According to an analysis of soil density data at the study site, it is understood that soil density has a direct relationship with increasing soil depth to 24 cm. Bd changes during the studied period have been shown with maps, such as was previously done by Raper [52]. Soil compaction is more evident in the upper layer (0-6 cm), in agreement with the results of Jien [53]; nevertheless, in our persimmon plot all the layers shown an increase in Bd.
Finally, in this study, we sampled in each place using one core per depth (in a total of four samples). It is well known, considering the literature above cited, that it is a usual and necessary a large number of samples. However, we also consider the core method to be invasive and destructive. The core method does not allow many repetitions in combination with other in situ monitoring technique, as the soil is disturbed during the sampling, which would modify the ISUM results. In this way, we state that there is still a further way to be researched, such as the deeper soil sampling; for example, 1 m would be a good option to confirm or not the hypothesis developed in this paper, since there are also several ecosystem service functions related to soil biology and hydrology that could be assessed.

The ISUM Data
According to the results of the ISUM map, soil compaction has taken place after the cultivation of persimmon, as the mean of ISUM values and the map reveals soil lowering (Figure 7). Soil topographical changes using ISUM show that soil lowering is less under the trees. This can be a consequence of the effect of the root system that reduces the soil density and lowering, as other researchers have demonstrated [54][55][56]. In the persimmon orchard in L'Alcúdia de Crespins, the key issue is the increase in Bd in the inter-row area, due to the passage of machinery, and although the roots can affect the Bd in the row area after three years, this cannot be the key factor. After the three years, the roots of the persimmon trees are still developing [57].
The sampling density (number of sampling points) of ISUM between a pair of trees is higher than with core sampling. We measured 48 points with ISUM and 5 (at 4 depths) with the traditional core method. As a result, the accuracy of topographic changes with ISUM is higher than for the core method.

Integrated Analysis of the Methods (ISUM and Bd)
The comparison of the calculated lowering soil for maps and dataset using bulk density data set and ISUM shows there is a little difference among them: 2 mm for raster maps. Those differences can be because of the mapping techniques and to the fact that the number of measurements of ISUM is larger than with the core ring methods. The core cylinder method needs time to be applied; ten times more than ISUM. In addition, the core method is invasive and destructive. The core method does not allow repetitions, as the soil is disturbed during the sampling. This reduces the chances of replication; meanwhile, ISUM allows as many repetitions in the measurements when they are necessary. The lowering soil maps developed with ISUM show more changes in the soil relief, as they better survey the micro-topography and soil roughness.
The differences in the assessments of soil lowering and compaction after the persimmon plantation using ISUM and core method could also be different since ISUM measures the total soil compaction, and the core method was applied from 0 to 24 cm of depth. The soil at the persimmon terraces is 2 m deep, and we suspect that there is also a process of compaction deeper than 24 cm depth [58]. The results of the ISUM are satisfactory, due to the low differences with the traditional core method, but the scientific community should test the accuracy of the ISUM in other terrains and under different crops. It will be also necessary to develop laboratory measurements under controlled conditions, to confirm that the soil compaction can be surveyed with ISUM.
The main question that we tried to answer in this paper was whether ISUM was the appropriate method to investigate soil compaction. Our research compares two methods of mapping soil compaction in a persimmon orchard planted on an agricultural terrace (level land surface). The traditional method of determining Bd by core sampling is time-consuming and costly. Given the use of grafted stock in the orchard, the possibility of applying a different measure of land surface change (ISUM) allowed for rapid estimation of changes in soil compaction. However, core sampling has the advantage of providing bulk density values at 6-cm intervals to a depth of 24 cm, so that the pattern of soil compaction with depth could be established, as well as modelling estimates of land surface changes in ISUM. Results of the two methods of ISUM and core sampling gave reasonably similar average results for land surface change, and it was concluded that ISUM could be used more widely as an indicator of soil compaction on level agricultural land.
It is well known that machinery and pedestrian pass increase soil compaction, which affects ecosystem services. The study uses average values and modelling to indicate approximate comparability between the results of the two methods. Either method will provide averaged results over an area, or both show considerable variability over short distances. Therefore, one question can be formulated: Which is more important-a spatial average of compaction, the detailed pattern of compaction, or identifying already known potential problem (highly trafficked) areas and comparing those with others? The hypothesis of this study was to test ISUM, which covers larger areas and even describes more spaces between rows and inter-rows. However, we cannot obviate that both methods should be applied to verify that soil movements and bulk density have the same spatial and temporal pattern.
As far as we observed, no soil loss (or minimum rates) takes place due to erosion from the terrace. Hence, the lowering of the soil surface can be attributed to soil compaction. However, some soil redistribution processes on the terrace could not have enough taken into account in specific parts of the terrace, e.g., by erosion due to the micro-topography, by splash erosion from raindrops, by soil burrowing animals, or by lateral movement of the soil material by trampling or tractor passes. These processes could cause a lateral redistribution of soil, resulting in a lowering of the soil surface at the point of measurement which could be then erroneously attributed to soil compaction. We assume that rainfall events during the first year of cultivation (bare soils in the terraces) occurred and generated the highest soil compaction and a quick increase in Bd. However, due to the flatness and structure of the terrace, the runoff discharge out of the cultivated terrace was null (fieldwork observations), and all the water was infiltrated upon our survey and the farmer's communication. In addition, the vegetation cover between the rows, despite the microtopography, could make these values even lower, so that we considered this as negligible. On the other, it is clear that lateral movement of the soil to the side by the tractor wheel must be considered, but our results (ISUM maps) demonstrated that soil compaction played a more important role on the trajectory movement than the changes caused by horizontal soil displacement by the tractor passes not showing extra differences below the trees. Our research also contributes to bringing a new method that can help in the research related to the measurement of the connectivity of flows, which are relevant to understand the changes in the soil system and the geomorphological processes [59,60], and for soil and water conservation in rural areas [61][62][63].
Ecosystem services are a key topic to promote sustainable management of the Earth, and this should focus on regions of the world that are relevant landscapes for the sustainability of the humankind, from the coastal land [64] to metropolitan areas [65], with special attention to the dynamic land-use changes we are assessing in the last century [66][67][68].
Agricultural land acts as an ecosystem that is managed by humans and produces fuel, fiber, food and landscapes for human use. The ecosystem services delivered by agricultural land produce services that the Millennium Ecosystem Assessment classifies as a provisioning service to humankind [69]. Agriculture also receives ecosystem dis-services that constrain productivity [70]. A clear example of dis-services is the one researched here due to the impact of soil compaction on the loss of infiltration rates and reduction in biota activity in soils [71,72]. The services (or dis-services) affect how agricultural ecosystems contribute to humankind's sustainability. And the diversity, functioning and composition of agriculture landscapes are definitive to achieve a sustainable society. Our research here determines that the soil compaction is an active process in agriculture land that can be surveyed with the ISUM method. Thanks to the low cost and quick survey of ISUM, the assessment of the changes in soil compaction and the consequent loss of services can be surveyed actively, and then contribute to finding sustainable solutions.

Conclusions
To achieve agriculture sustainability and conserve ecosystem services in this kind of agroecosystem, soil compaction must be assessed. This study evaluated the improved stock unearthing method (ISUM) as a low-cost and precise alternative to the tedious and costly traditional core sampling method. ISUM samples were taken at 10 cm intervals from tree to tree, and our results demonstrate that soil compaction in persimmon orchards is due to soil compression because of tractor passes, which directly affect the soil topography. The comparison of the traditional ring method and the ISUM showed a slight overestimation of the soil level with ISUM (from 12 to 14 cm of soil consolidation on average, respectively). However, ISUM was found to be very efficient, as the labor cost is one order of magnitude lower (24 h) than core sampling (272 h). The use of a simple measurement tool such as ISUM can contribute to increasing the number of measurements and assessing soil compaction in terraced landscapes, where surface wash does not take place. ISUM is a suitable alternative to assess soil compaction in agricultural land and ecological service evaluations. Using the ISUM method can help in soil consolidation and soil compaction studies. Assessments with ISUM could be combined with core samples, at least in the areas with the lowest and highest changes in soil terrain micro-topography, as they are the ones with more changes. This will allow an understanding of the bulk density changes, and ISUM will considerably decrease the number of samples taken.