Improved Stock Unearthing Method (ISUM) as a tool to determine the value of alternative topographic factors in estimating inter-row soil mobilisation in citrus orchards

Uso del Método Mejorado del Uso del Injerto (ISUM) como herramienta para determinar el valor de factores topográficos alternativos en la estimación de la movilización del suelo entre hileras para huertos de cítricos Uso do método melhorado de Stock Unearthing (ISUM) como ferramenta para determinar o valor de fatores topográficos alternativos para estimar a mobilização do solo entrelinhas em pomares de citrinos Rodrigo-Comino J. @, 1,21 rodrigo-comino@ uma.es


Introduction
Recent scientific research confirms that cultivated lands can be affected by land degradation processes (Baude et al. 2019;Gomes et al. 2019;Perovic et al. 2019). The Mediterranean areas growing fruits crops such as citrus, olive, and grapes are the ones that favour high erosion rates and losses of nutrients (García-Ruiz et al. 2013;Panagos et al. 2018). However, most of the research carried out takes place for a few months or years with soil erosion plots or rainfall simulation experiments (Ouyang et al. 2018;Wynants et al. 2019). Those data contribute to evaluating the impact of land management, but little information is found about soil mobilisation rates over the long-term, which is necessary to characterize and evaluate the problem and find solutions (Amiri 2010).
are cultivated worldwide and they are also expanding in the Mediterranean territories, where clementine varieties are growing due to the premium prices in the European markets. This expansion also results in a new allocation of the citrus plantations and as a result, a potential increase in soil erosion is being detected (Cerdà et al. 2018a).
New clementine orchards expand uphill from the bottom of the valleys as the alluvial plains and fluvial terraces are changing land use from agriculture to urban areas. Thermic inversion in winter and frosts may take place, which heavily damages the citrus production and the trees (Cerdà et al. 2018b). The new plantations are highly mechanized and chemicals, such as the glyphosate and other herbicides are applied as the main strategy to remove the herbs (Abit et al. 2012;Donkersley et al. 2018).
The farmers wish to maintain a bare soil surface during some parts of the year, and the continuous application of glyphosate takes place every few weeks to leave the soil free of weeds. However, it is well-known that vegetation cover is the key control factor of soil erosion and bare soils can enhance the highest erosion rates (Bagio et al. 2017;Chen et al. 2019). In Mediterranean fields, some farmers dislike seeing any weed and tend to prefer the soil bare as they accept this as the best practice (Taguas et al. 2015;Villanueva et al. 2015). As a result, soil losses are enhanced and achieve non-sustainable values (Marqués et al. 2015;Sastre et al. 2017).
Soil properties and erosion process are influenced by several human and environmental factors from very small to very large spatial scales. It has demonstrated its impacts over a similarly wide range of spatial scales, developing numerous methods to estimate and predict soil erosion processes such as models (Hrabalikova and Janeček 2017). Within them, the most applied model is the Universal Soil Loss Equation -USLE (Wischmeier and Smith 1978) and up to now, it is the most used model worldwide to predict the soil losses in experimental plots, hillslopes, farms, watersheds and even, regions. The roughness (random or oriented) and topography (micro and macro) constitute elements that affect the spatial organization of runoff (Govers et al. 2000;Snapir et al. 2014). The micro-topographical and macrotopographical factors are important elements that can drive the erosion process (Sun et al. 2014;Marquard et al. 2016;Ochoa et al. 2016), but there is little research about how slope angle affects soil erosion in intensively tilled agricultural fields (Meshkat et al. 2019). The USLE topographic factor is defined by the junction of the slope (S) and length (L) factors. According to Wischmeier and Smith (1978) "the slope length (L) constitutes the distance from the point of origin of the surface flow to the point where each slope gradient (S) decreases enough for the beginning of deposition or when the flow comes to concentrate in a defined channel".
In addition, there is an equation designed for determining the LS factors (usually computed together) to be obtained when the local topographic characteristics are different from those ones originally established for the traditional USLE (Qin et al. 2018;Schmidt et al. 2019). However, the equation has been adapted to accomplish regional particularities. This allowed the creation of new methods and procedures to determine the LS value for predicting the soil loss by the USLE, although it is well-known the elevated time and efforts needed to obtain LS data (Hrabalikova and Janeček 2017).
On the other hand, other alternatives to obtain topographical metrics of a specific terrain by means of drones or/and high-precise cameras have been also developed (Esposito et al. 2017;Nadal-Romero et al. 2015). However, the expensive costs, the permissions, the influence of shadows in the quality of the database, and structural and functioning limitation (short term batteries, for instance) make these procedures not available for all of the research groups and countries (Eltner et al. 2016;Beretta et al. 2018).
Another possibility to solve these issues could be the use of the Stock Unearthing Method (SUM). The SUM is a procedure developed to estimate soil mobilisation, as well as sediment deposition rates taking into account the measures of the graft union in vineyards (Brenot et al. 2008;Casalí et al. 2009). This method has currently improved in terms of detailing and precision (incorporating measurements in the inter-row areas) released by  and , updated to the Improved Stock Unearthing Method (ISUM). It needs to use graft unions for certain species (citrus, persimmons, almonds, etc.), where it provides a method for erosion measurement (Bayat et al. 2019;da Silva et al. 2019).
After obtaining the data through field work using ISUM, the database has the versatility to be used for several goals. For example, it was also suitable for mapping rills and sinks, and soil depletion and accumulation combined with a terrestrial photogrammetric scanner system (Remke et al. 2018). Such example shows that the ISUM provides a noteworthy ability to evaluate different timescales and objectives that range from single hydrologic events to decennial scales, as well as different spatial scales ranging from the plot to hillslopes (Rodrigo-Comino et al. 2019b).
To our knowledge, the use of database generated by the ISUM for computing the values for LS-USLE factor has not been tested yet, and this kind of work might facilitate the surveying of local data and in a cost-effective and accurate manner. Moreover, no research has been developed to date in order to test if the use of the graft union of the plants as a passive indicator to measure LS-USLE factor could be also used for other types of grafted plants such as the citrus orchards. Therefore, we aimed to check the feasibility of the use of the database surveyed (measurements in the graft union and inter-row areas) in the fieldwork by means of ISUM for calculation of LS factor to be used in soil mobilisation estimations by the Universal Soil Loss Equation or other soil loss-related models in a clementine plantation located in eastern Spain.

Study area
The field site (Figure 1a) is located in the Canyoles river watershed at 170 m a.s.l. in the municipality of Canals, province of Valencia (Spain). This plantation was selected because it was considered as a representative zone of the Mediterranean citrus plantation that was established in the 1960s in eastern Spain on terraces to allow the flooding of the fields. The construction of the terraces does not necessarily mean that soil erosion is completely stopped but can be disconnected, defaulted or enhanced depending on the rainfall event (Arnáez et al. 2015;Li et al. 2014). The climate is typically the Mediterranean, with a mean annual precipitation of 532 mm and a mean annual temperature of 16.2 ºC. The soil can be classified as a Xerorthent (Soil Survey Staff 2014) and the soil texture is sandy-loam: 20% clay, 29% silt, and 51% sand (USDA Classification; Soil Survey Division Staff 1993). The study area was selected in a general sloping terrain (6% slope) close to the Canyoles River, with very low rock fragment content (2.4%).
In the 2000s, the terraces were removed and transformed in a sloping terrain where drip irrigation and highly mechanized agriculture was applied. The clementine plantation (variety Tango) is a unique 30 ha farm where land levelling was done in 2008. Then, a furrow of 69 cm in height and 6 m width was established. The trees were planted at the highest point of the furrow. The crop is drip-irrigated and herbicide (Glyphosate) is applied any time the weeds germinate, which result in a bare soil surface during several parts of the year. The use of herbicides is very intense in summer.
The irrigation system pumps water from the Caroig aquifer and it uses two lines of 16 mm diameter pipe per row with drips each 1 m. Irrigation takes place every day in summer, and twice a week since March when the harvest is finished. The management until the time of soil management was similar in all plots: herbicide (Glyphosate (N-(phosphonomethyl)glycine) and inorganic fertilizers (NPK, 1.2 Mg ha -1 per year).

Experimental design
We conducted the fieldwork according to the explanations provided by Rodrigo-Comino and . Initially, we located each pairedgraft stock union, which were initially planted 7 cm above the soil surface as confirmed by the farmers during informal interviews, and we measured a control plot with recent plantations with a variation of about ± 1 cm. We then stretched a string between such elevated points on the graft unions. By means of a plastic meter stick, we measured and registered the depth from the reference string to the soil level at each 100 mm distance from the clementine tree trunk (from the two clementine tree trunks connected by the string), totaling thirty transects between each paired-tree (Rodrigo-Comino et al. 2019b). This enabled us to record all points either below or above the ground surface. The planting pattern was 6 x 2 m, the usual pattern for citrus in the new plantation in this agricultural area. Hence, we completed a total of 1,800 measurements in the study area. Figures 1b and 1c show a very common practice in fruit trees: forming ridges. The soil of the inter-row centre is placed in the row of trees. Therefore, the great furrow with gentle micro-topography of the inter-row centre could be not due to erosion. In any case, it would be tillage erosion, but since there is no soil exit from the plot, this rating is doubtful. We recorded and treated our database in digital worksheets. The data were organized in a manner that the position X corresponds to the rows, Y corresponds to the position of the trees, and Z corresponds to the depth of the soil surface in relation to reference rope. In the field, we observed some differences in the characteristics of each side of the inter-row. All transects presented an inter-row (furrow) V-shaped (see figures ahead). Hence, for each transect, we identified the deepest point of the channel and we divided the transect into the "left side" and "right side".
After that, for the calculation of the LS value for each side of each transect, we initially needed to calculate the percentage of the slope, as well as the length of the surface of each side of each transect. The percentage of the slope was determined by the quotient among the distance of the stretch of the string from the tree and the local where the deepest value was detected and the value of deep computed for this point and multiplied by 100 (Figures 2b and 2c).
The superficial length was determined by means of calculating the sum of the values of the hypotenuse of all of each segment of the side of the transect (Figure 2a). After determining the values of the percentage of slope and length of each side for each transect, we computed the LS values considering two equations: Where (for both equations) LS represents the LS-factor of the USLE, taking into account the local topography; l means the length of the segment of the transect (m) and s is the slope of the segment of the transect (in %).
These two algorithms for LS computation were chosen for testing the data surveyed by means of the ISUM used in different conditions according to the landscape features and local soil management actions (Bertoni and Lombardi Neto 2014). Hence hereinafter, the equation (1) will be presented as LS1 and equation (2) as LS2.
We also computed, for each side of each transect, the Transect Length Index (TLI), which is a variation of the Tortuosity Index published by Boiffin (1984). Firstly, by using Pythagoras's theorem equation, we computed the length of the segment. After, we summed the values of all segments. Next, we calculated the TLI through of division of the values of length of the soil surface (G) and the straight line parallel to the soil surface (H). The procedure is depicted in Considering the furrow V-shaped and taking into account the width and depth, we estimated the area of the section at the moment of the construction of the furrow. This means that the original furrow area section had already been calculated taking measures in a close plot with a recent plantation designed and tilled by the same farmer and machinery. So, we consider this reference to the shape and depth of the original furrow (zero moment) as the crucial point since all the data taken to determine the current micro-topography will be referenced to this original standard furrow (Bayat et al. 2019). Furthermore, using the database generated by ISUM, we estimated the area of the section in the moment of the data survey. We assume that the difference of the areas corresponds to the volume of mobilized soil since the construction of the furrow and planting of the trees (19 years) and by tillage.

USLE estimations
For estimating the soil loss by means of the USLE, we obtained data for each factor explained as follows: R factor-rainfall erosivity. We considered the data estimated by Panagos et al. (2015a). Such authors presented data and maps of annual values of rainfall height (P), rainfall erosivity (R), and rainfall erosivity density, this last is the ratio between R and P. These data are average values for large areas and therefore, only applicable to very small scales if the values are similar with other local and old values estimated by regional agencies. In this paper, we work at the microtopography level but using citrus plantation with similar topographical characteristics, so we consider that this is the most current reference of this study area. In addition, the use of this reference could make this paper comparable to the papers published by the above-mentioned research group in Europe in the future.
K factor-soil erodibility. Based on the granulometric composition of the soil (20% clay, 29% silt, and 51% sand) and using the USDA textural triangle, we classified the soil as loam. Afterwards, we checked that the local soils have 2.7% organic matter and bulk density of 1.3 g cm -3 . It is important to remark that the soils are relatively similar in the first 20 cm soil depth due to the tillage (Trueba et al. 1998;Bayat et al. 2019). Using all these data and using a table provided by Morgan (1986), we estimated the value of the erodibility of the local soil.
L and S factors-length and slope. They were estimated by ISUM method and using the

equations (1) and (2).
C and P factors -soil cover and conservation practices. Because the soil was almost completely uncovered (intentionally) and there were no soil conservation practices registered in the studied area, we attributed values 1.0 for both factors, following criterion presented by Morgan and Nearing (2011). When determining soil loss by USLE, we decided to take the factor C equal to 1. Fruit tree canopy provides a certain degree of coverage but during the fieldwork period, it was almost inexistent. Moreover, we mentioned that there is in some places enough herbaceous vegetation to justify the erratic orientation of runoff but only after some rainfall events and periods where herbicides are not used, which did not invalidate this C-factor value.

Statistical analysis
Using the free statistical packages R (CRAN R project, http://cran.r-project.org) and Bioestat 5.0 (Ayres et al. 2007), we calculated the main descriptive statistics. After that, we checked the pattern of distribution of our data using the nonparametric test of variance (Kruskal-Wallis test) for all metrics and the Z-transformation to normalize the data, as well as the nonparametric correlation analysis (Spearman´s rank coefficient test).

Topographical analysis by means of the use of the ISUM data, calculations of the USLE's topographic factor (LS) and Transect Length Index (TLI)
From the 1,800 points measured, we observed an overall mean depth value of 495 mm and a maximum value of 910 mm, this one being the deepest point observed among all the transects. The coefficient of variation was 38%. We observed a predominance of values ranging from 701 to 800 mm depth, with 18% of occurrences (Figure 3). However, the modal value was 630 mm.
After converting our database into line charts, we observed that the studied inter-row is concaveshaped (V-shaped -shallow and wide) along with the whole extension of the channel (Figure   4).  the Kruskal Wallis test revealed significant differences between the sides for all parameters, evidencing asymmetries and irregularities along the studied channel ( Table 2).
In the field, we initially had the impression of a perfect symmetry among the sides of the same transect and the descriptive statistics also suggest such similarity (Table 1). However, Length in m; slope in %; LS1 (l means the length of the segment of the transect and s is the slope of the segment of the transect in %; 1 corresponds to equation 1), LS2 (l means the length of the segment of the transect and s is the slope of the segment of the transect in %; 1 corresponds to equation 2), and TSI (Transect Length Index) are dimensionless; C.V. -coefficient of variation. evidenced by the higher values of variance, standard deviation, and also the coefficient of variation.
Both sides presented slope values ranging from 20 to 35% and were considered sloped. The data of the left side presented a higher variation, LS1 (l is the length of the transect segment and s is the slope of the transect segment in %; 1 corresponds to equation 1), LS2 (l is the length of the transect segment and s is the slope of the transect segment in %; 1 corresponds to equation 2); and TSI (Transect Length Index).
available: the rainfall erosivity (R factor), soil erodibility (K factor), land cover (C factor) and soil conservation practices (P factor Concerning the correlation analysis, the metric slope was significantly correlated with all the other metrics, except for TLI in both sides ( Table  3). There was no metric significantly correlated Italic: < 0.001; italic and bold: < 0.0001, values neither italic nor in bold were not significant at p = 5%. LS1 (l is the length of the transect segment and s is the slope of the transect segment in %; 1 corresponds to equation 1); LS2 (l is the length of the transect segment and s is the slope of the transect segment in %; 1 corresponds to equation 2); and TSI (Transect Length Index).
with the final values for the LS values the left side (for both algorithms). For the data of the right side, the metric length showed less correlation than the final value for the algorithm LS1 and no significant correlation for LS2.
For the two tested algorithms (Table 4), the final values of factor LS were positively and significantly correlated with the slope, while the factor length was slightly negatively correlated.
The metric length showed a higher correlation Italic: < 0.001; italic and bold: < 0.0001, values neither italic nor in bold were not significant at p = 5%. LS1 (l is the length of the transect segment and s is the slope of the transect segment in %; 1 corresponds to equation 1); LS2 (l is the length of the transect segment and s is the slope of the transect segment in %; 1 corresponds to equation 2); and TSI (Transect Length Index). area (360 m 2 ), the volume of mobilised soil, and the bulk density for the local soil (1.3 g cm -3 ), we estimated a total soil mobilisation of 8.3 mm yr -1 or 10.4 kg m -2 yr -1 ( Table 5).

By means of USLE
The extension of the studied area allowed us considering a unique spatial value of the area for 3.2. Soil mobilization estimation

Directly by means of ISUM
The difference among areas predicted in the moment of the construction of the furrows and the moment of the data survey allowed us to estimate a total volume of 56.9 m 3 of soil mobilised in 19 years. Considering the studied both for R and K factors ( Table 5). For R factor, the values that we considered for computing the soil loss were: as minimal R-value: 730 MJ mm ha -1 h -1 y -1 ; as medium R-value: 815 MJ mm ha -1 h -1 y -1 ; and as maximum R-value: 900 MJ mm ha -1 h -1 y -1 . The value of the K factor (soil erodibility) was 0.0295 Mg ha h ha -1 MJ -1 mm -1 . Soil mobilisation by ISUM (kg m -2 yr -1 ) 10.40 LS (l is the length of the transect segment and s is the slope of the transect segment in %).

Topographical analysis by means of the use of the ISUM data, calculations of the USLE's topographic factor (LS) and Transect Length Index (TLI)
The results reached by the use of soil erosion models have a high connection with the spatial resolution of DEM and the LS factor (Panagos et al. 2015b). Hence, we argue that ISUM could be considered as a useful alternative to get data for using in soil loss predictions, especially for hillslope scales as experimental plots or part of crop fields. The acquisition and analysis of topographic data carried out in the digital environment have become a custom in erosion modelling since they enable efficient analysis of a given area (Hrabalikova and Janeček 2017). For our study, the data surveyed by means of ISUM have high versatility to be handled in several ways according to the interests and goals of the users.
Concerning the values of declivity (% of the slope), we found values considered elevated. Such values could be classified as stronglywaved in terms of percentage of inclination (for example, areas with a slope greater than 20%). If we consider the statement "the runoff increases as the gradient slope increases" (Ji-Jun et al. 2012), we can infer, for our study area, that it is expected that the right side of the inter-row produces a little bit more runoff than the left side, once that the right side is more sloped than the left side.
For the LS values calculated by means of algorithm LS1, the results were higher than 1.00 for all transects in both sides (left and right). This means that, at least, mathematically, this factor could play a key role as an intensifier of the soil loss as demonstrated by other authors under field or laboratory conditions (Bryan 2000;Silva et al. 2011;Nadal-Romero et al. 2014). For the algorithm LS2, in 43% of the transects (left side), the final LS value was lower than 1.00, in 70% of the transects (right side) the final LS value was lower than 1.00. For the right side, in one transect we found an LS value of 1.00, meaning that in this area the topography plays a neutral role in the soil loss process (neither attenuator nor intensifier).
The data calculated by means of the Transect Length Index (TLI) reveal aspects of the roughness of the soil surface. In our study area, such data evidence that the soil surface presents a weak and random pattern of roughness. This pattern of roughness is influenced predominantly by the arrangement of clods and stones on the ground surface (Govers et al. 2000;Nearing et al. 2017).
We applied two models to determine the value of the LS factor and analyze their results. It would be desirable that we can conclude which of the studied models to calculate LS factor (LS1 and LS2) is the one that provides the best results for the Mediterranean conditions. However, we consider that it is difficult to decide on one of them after testing this procedure for the first time in an agricultural field and only using one transect. However, we can conclude that the LS method using a database conformed by in situ measurements is more accurate and realistic than the other methods traditionally used with a digital elevation model, which depends on the resolution, expensive costs (drones, remote sensing data, etc.) and permissions.

Soil loss predictions
All values predicted by USLE were lower than the predicted directly by ISUM. Such discrepancy of data may be explained especially by the rainfall and soil characteristics. The most interesting information can be extracted from the high rainfall intensity, which is frequent in this area. Some studies that present high erosivity density suggest that the rainfall is characterized by rainstorms of high intensity and short duration (Panagos et al. 2015a). The local value of this index is one of the highest of the continent, and Spain is one of the countries highlighted by Panagos et al. (2015a) as having highly erosive rainfalls.
For the soil erodibility, according to the interpretation table provided by Silva et al. (2007), the value is classified as medium erodibility and this situation is not worse thanks to the amount of soil organic matter that exists in the soil. Hence, it is a factor that the landowners should pay attention to in order to maintain production, i.e., correct soil management to conserve soil organic matter and reduce erodibility such as it was shown in a long-term trial by Novara et al. (2019a).
Related to soil management and soil cover management, we consider that the scarce community of herbaceous plants that grows up the sides of the inter-rows (plants up to 10 cm tall, several species) could certainly play a role in terms of providing some protection to the soil and also in determining the erratic orientation of the runoff in moments of heavy rains, although the vegetation occurs irregularly along the terrain (in some areas the soil surface is scarcely covered, in others not) (Haddaway et al. 2016;Van Hall et al. 2017). In the field, we reported that there is no occurrence of linear erosion processes such as side rills and side gullies, and this corroborates with the pattern of roughness calculated by means of TLI (Zhao et al. 2018).
Moreover, another important recommendation to be made to the local farmers in terms of soil management is the conservation of the herbaceous plant community, in order to avoid the activation of erosion process and minimize soil degradation (Chau and Chu 2018;Ruiz-Colmenero et al. 2013). For the studied region (Mediterranean Region), if the soils remain bare, they are liable to soil sealing and crusting (García-Estringana et al. 2010).
Previous studies agree that if the soil remains covered, the infiltration rate might be some order of magnitude greater than when the soil remains uncovered, since the herbaceous plants would supply litter to the soil to increase its organic matter content ). Furthermore, their complex root systems could induce improvements in soil properties, for example, the increase of soil macro-porosity (Di Prima et al. 2018). Such alterations would minimize the deleterious effects of runoff and soil losses, beyond the direct effects of their aerial cover (García-Estringana et al. 2010. Therefore, in the future, it could be interesting to see if the vegetation cover remains, whether ISUM results are the same or change from this research. Comparatively, Bagio et al. (2017) also predicted soil loss by USLE and RUSLE and when compared, they also reported underestimated values. They justify the underestimation of soil loss by the USLE model explaining the high influence of the R factor, due to the discrepancy between the rainfall regime in the conditions where the USLE was created and the prevailing climate in the research area. We confirm the statement of such authors also for our study area, although the R factor can be adapted to each study area (Benavídez et al. 2018).
Finally, we want to remark another important point to be further studied in the future. Figure 1c shows a formation of ridges (very common soil management practice in fruit trees). Therefore, it is difficult to discern how much soil has been lost by erosion due to extreme rainfall events, floods or how much soil has been displaced by tillage or compacted such as other authors found Novara et al. 2019b).
Even so, it is important to undertake more indepth research in order to distinguish how much the soil was affected by each specific natural or anthropogenic process, in order to implement the most efficient soil control protection measure.
There is a preference of the landowners to keep the soil surface uncovered, perhaps to facilitate the harvest and tillage. However, this is a mentality that should be re-orientated to conserve at least a minimum vegetation cover, in order to conserve the soil, especially in the sloped regions of the channels. This could minimize the generation of runoff, increase the water infiltration, avoid soil loss, and retain sediment or soil particles that were eventually detached from their original place.

Conclusions
The LS-factor has been successfully calculated and with great detail using a database generated by means of ISUM procedures. This was achievable because of the possibility of computing slope and length for terrain by means of ISUM. The sides of the studied inter-row are significantly different in terms of symmetry. We also verified that the sides of the inter-row present a pattern of roughness that favours the erratic formation of runoff. Therefore, we conclude that the application of the ISUM method as a tool to determine the value of alternative topographic factors can be considered for future research projects.

Acknowledgements
The research leading to these results has received funding from the European Union Seventh Framework Program (FP7/2007(FP7/ -2013 under grant agreement n° 603498 (RECARE project).
• Amiri F. 2010. Estimate of erosion and sedimentation in semi-arid basin using empirical models of erosion potential within a geographic information system. Air, Soil and Water • Rodrigo-Comino J, Keesstra S, Cerdà A. 2018. Soil erosion as an environmental concern in vineyards: The case study of Celler del Roure, Eastern Spain, by means of rainfall simulation experiments. Beverages 4(2):31.