Sensing and Mapping the Effects of Cow Trampling on the Soil Compaction of the Montado Mediterranean Ecosystem

The economic and environmental sustainability of extensive livestock production systems requires the optimisation of soil management, pasture production and animal grazing. Soil compaction is generally viewed as an indicator of soil degradation processes and a determinant factor in crop productivity. In the Montado silvopastoral ecosystem, characteristic of the Iberian Peninsula, animal trampling is mentioned as a variable to consider in soil compaction. This study aims: (i) to assess the spatial variation in the compaction profile of the 0–0.30 m deep soil layer over several years; (ii) to evaluate the effect of animal trampling on soil compaction; and (iii) to demonstrate the utility of combining various technological tools for sensing and mapping indicators of soil characteristics (Cone Index, CI; and apparent electrical conductivity, ECa), of pastures’ vegetative vigour (Normalised Difference Vegetation Index, NDVI) and of cows’ grazing zones (Global Positioning Systems, GPS collars). The significant correlation between CI, soil moisture content (SMC) and ECa and between ECa and soil clay content shows the potential of using these expedient tools provided by the development of Precision Agriculture. The compaction resulting from animal trampling was significant outside the tree canopy (OTC) in the four evaluated dates and in the three soil layers considered (0–0.10 m; 0.10–0.20 m; 0.20–0.30 m). However, under the tree canopy (UTC), the effect of animal trampling was significant only in the 0–0.10 m soil layer and in three of the four dates, with a tendency for a greater CI at greater depths (0.10–0.30 m), in zones with a lower animal presence. These results suggest that this could be a dynamic process, with recovery cycles in the face of grazing management, seasonal fluctuations in soil moisture or spatial variation in specific soil characteristics (namely clay contents). The NDVI shows potential for monitoring the effect of livestock trampling during the peak spring production phase, with greater vigour in areas with less animal trampling. These results provide good perspectives for future studies that allow the calibration and validation of these tools to support the decision-making process of the agricultural manager.


Introduction
The silvopastoral ecosystem characteristic of the Iberian Peninsula, known as the Montado in Portugal and the Dehesa in Spain, is a mixed system that integrates, in the same place, agronomic crops (e.g., cereals, forage, or pastures), a tree stratum (e.g., Oak trees) and livestock [1] and has been proposed to optimise economic and environmental benefits, CI values are highly influenced by diverse soil factors: intrinsic (e.g., soil moisture, bulk density, texture and structure) and extrinsic (e.g., management system) [19]. Moreover, the results of field penetrometers depend on user operating speed (penetration rate), which is often challenging to standardise; a change in operating speed alters the force the users apply to insert the equipment rod, which, in turn, changes the result [15,18]. On the other hand, the characterisation of this and other soil properties is a difficult process due to the high soil spatial variability [14] and the interaction and combination of the factors involved, particularly in silvopastoral systems [1]. The soil's superficial micro-variability is mainly controlled by soil and crop management practices, plant roots, wet/dry cycles and, in non-tillage systems, by surface-sealing processes [19].
In this context of high soil spatial variability, it is essential for management to take advantage of the technological developments associated with Precision Agriculture (PA) [21,22]). The objectives of PA are to optimise production by increasing yield or reducing costs, minimise the use of natural resources, reduce the environmental impact and improve soil quality [14]. To achieve these goals, many new technologies have been developed and used for sensing and mapping crops and soils [14]. Thus, the use of new technologies for the correct management of farms is important to prevent the degradation of new areas, since the use of pastures is a practical alternative for feeding ruminants and, concomitantly, for producing meat and/or milk [7]. One of the most widely used parameters for the management of soil spatial variability is the mapping of apparent soil electrical conductivity (EC a ), which allows the farmer to identify differentiated management zones (e.g., for differential fertiliser application, soil amendment or irrigation) [22,23]. The EC a is defined as the soil's capacity to conduct electric currents, which is influenced by many physical-chemical features of the soil, including the clay content [14,24,25]. The temporal stability of this parameter has allowed the recognition of the geospatial measurements of EC a as a valuable mapping tool that indicates the soil's potential productivity [14,26]). Simultaneously, the measurement of the EC a can be considered a relatively inexpensive, easy and fast technique with the potential to contribute to the identification and prediction of spatial variability of soil compaction [14,16]. However, in the recent literature, there is a lack of scientific papers regarding the relationship between these two parameters [14].
Other technologies that can be very useful for monitoring animal trampling are Global Positioning Systems (GPS collars). These collars have several applications in PA; for example, they can be used to monitor preferred grazing areas [27]. The geolocation of animals by remote sensing (RS) from satellites, at regular time intervals, allows the identification of grazing patterns and areas of higher grazing intensity throughout the vegetative cycle of the pasture [27]. These preferential grazing areas are potentially at greater risk of compaction by animal trampling, especially in periods of higher precipitation, such as winter or spring in regions with a Mediterranean climate. The use of RS imagery based on the Sentinel-2 satellite to obtain vegetation indices, namely the Normalised Difference Vegetation Index (NDVI), also proved a promising tool to express the response pattern of pastures' vegetative vigour [7,[28][29][30]. Therefore, with the vegetation indices that can be obtained from digital images and are sensitive to changes in the vegetation cover of pastures, before and after grazing, it is possible to monitor and to identify degraded areas with overgrazing, as well as areas that are arid or without vegetation, for example, in large or small fields and at different time scales [7]. This study aims: (i) to assess the spatial variation in the compaction profile of the 0-0.30 m soil layer over several years; (ii) to evaluate the effect of animal trampling on soil compaction; and (iii) to demonstrate the interest of combining various technological tools for sensing and mapping indicators of soil characteristics (CI and EC a ), of pastures' vegetative vigour (NDVI) and of cows' grazing zones (GPS collars).

Site Description, Field Management and Sampling Scheme
The field of the study (Figure 1) is located at the Mitra experimental farm of the University of Évora in the southern region of Portugal (coordinates 38 • 32 10"N; 7 • 59 80"W). This field of Quercus ilex ssp. rotundifolia Lam. and bio-diverse pastures has been used for extensive and rotational grazing of 60 adult cows of two native breeds ("Alentejana"-25 animals; and "Mertolenga"-35 animals) [5]. Of the total grazing area (about 100 ha), 20 ha were monitored (11 ha in "Field A" and 9 ha in "Field B"). Grazing was conducted to have a mean stocking rate of about of 0.6 head.ha −1 in both fields (A and B) between October and December. In January and February and between April and June, no animals graze in "Field B", while in March, no animals graze in "Field A". More details of this grazing management system can be consulted in Serrano et al. [5].

Site Description, Field Management and Sampling Scheme
The field of the study (Figure 1) is located at the Mitra experimental farm of the University of Évora in the southern region of Portugal (coordinates 38°32′10″ N; 7°59′80″ W). This field of Quercus ilex ssp. rotundifolia Lam. and bio-diverse pastures has been used for extensive and rotational grazing of 60 adult cows of two native breeds ("Alentejana"-25 animals; and "Mertolenga"-35 animals) [5]. Of the total grazing area (about 100 ha), 20 ha were monitored (11 ha in "Field A" and 9 ha in "Field B"). Grazing was conducted to have a mean stocking rate of about of 0.6 head.ha −1 in both fields (A and B) between October and December. In January and February and between April and June, no animals graze in "Field B", while in March, no animals graze in "Field A". More details of this grazing management system can be consulted in Serrano et al. [5]. The dominant soil type of this field is acidic and a not very fertile Cambisol [31], mainly used for mixed agrosilvopastoral systems.
The chronology of the measurements carried out in the experimental field is presented in Figure 2. In April 2018, the ECa and altimetric surveys, as well as soil sampling and CI measurements, were carried out. The CI measurements were carried out again in The dominant soil type of this field is acidic and a not very fertile Cambisol [31], mainly used for mixed agrosilvopastoral systems.
The chronology of the measurements carried out in the experimental field is presented in Figure 2. In April 2018, the EC a and altimetric surveys, as well as soil sampling and CI measurements, were carried out. March 2019, September, November and December 2021 and March 2022. Animal monitoring was carried out between January and May 2021. Vegetation Index (NDVI) time series reconstruction was performed throughout the 2021/2022 pasture vegetative cycle (between September 2021 and June 2022).
The monitoring area ("Field A" + "Field B") was sampled in two phases: the first in 2018-2019 and the second in 2021-2022. In the first phase, 24 Sentinel-2 pixels "10 m × 10 m" were georeferenced for sampling in areas without trees (outside the tree canopy, OTC; Figure 3a), with 12 in each field (A and B). In the second phase, half of these areas was sampled (12 sampling pixels, 6 in each field, A and B), as well as the area under the tree canopy (UTC) closest to each of these pixels (12 sampling trees; Figure 3b).    The monitoring area ("Field A" + "Field B") was sampled in two phases: the first in 2018-2019 and the second in 2021-2022. In the first phase, 24 Sentinel-2 pixels "10 m × 10 m" were georeferenced for sampling in areas without trees (outside the tree canopy, OTC; Figure 3a), with 12 in each field (A and B). In the second phase, half of these areas was sampled (12 sampling pixels, 6 in each field, A and B), as well as the area under the tree canopy (UTC) closest to each of these pixels (12 sampling trees; Figure 3b). The monitoring area ("Field A" + "Field B") was sampled in two phases: the first in 2018-2019 and the second in 2021-2022. In the first phase, 24 Sentinel-2 pixels "10 m × 10 m" were georeferenced for sampling in areas without trees (outside the tree canopy, OTC; Figure 3a), with 12 in each field (A and B). In the second phase, half of these areas was sampled (12 sampling pixels, 6 in each field, A and B), as well as the area under the tree canopy (UTC) closest to each of these pixels (12 sampling trees; Figure 3b).

Characterisation of the Climate
The climate of this Mediterranean region is classified as Csa (Köppen-Geiger classification) [32]. It is characterised by high inter-annual irregularity and low rainfall (<600 mm) that is more frequent in the autumn-winter period and practically nil in the summer [33].
The thermopluviometric diagram of the Évora meteorological station between July 2015 and June 2022 is presented in Figure 4. This also shows the monthly rainfall between July 2020 and June 2021 and between July 2021 and June 2022. The great irregularity of the rainfall distribution is evident: for example, 2020/2021 shows high accumulated rainfall in February (142 mm), October (141 mm), November (108 mm) and April (106 mm) and very low rainfall in March (18 mm), while 2021/2022 shows high accumulated rainfall in March (135 mm), October (113 mm) and December (97 mm) and very low rainfall in January (5 mm), February (10 mm) and November (15 mm). This irregularity and, especially, the occurrence of events with a high concentration of rainfall, associated with poorly drained soils, can lead to situations of flooding and, consequently, potentiate soil compaction by animal trampling.

Characterisation of the Climate
The climate of this Mediterranean region is classified as Csa (Köppen-Geiger classification) [32]. It is characterised by high inter-annual irregularity and low rainfall (<600 mm) that is more frequent in the autumn-winter period and practically nil in the summer [33].
The thermopluviometric diagram of the Évora meteorological station between July 2015 and June 2022 is presented in Figure 4. This also shows the monthly rainfall between July 2020 and June 2021 and between July 2021 and June 2022. The great irregularity of the rainfall distribution is evident: for example, 2020/2021 shows high accumulated rainfall in February (142 mm), October (141 mm), November (108 mm) and April (106 mm) and very low rainfall in March (18 mm), while 2021/2022 shows high accumulated rainfall in March (135 mm), October (113 mm) and December (97 mm) and very low rainfall in January (5 mm), February (10 mm) and November (15 mm). This irregularity and, especially, the occurrence of events with a high concentration of rainfall, associated with poorly drained soils, can lead to situations of flooding and, consequently, potentiate soil compaction by animal trampling.

Soil Apparent Electrical Conductivity (ECa) and Altimetric Surveys
With the aim of measuring ECa, a contact-type sensor (Veris Technologies, Salina, KS, USA) was utilised. Measurements at a depth of 0-0.30 m were performed in April 2018. An all-terrain vehicle was used to pull the sensor. The average speed of the vehicle was 2.0 m s −1 ; consecutive passages, spaced 10 m, were made across the field. The spatial resolution of the ECa measurements was a 2 by 10 m grid, since a measurement was taken every second. A global navigation satellite system (GNSS) antenna was installed near the sensor. The obtained data were used to produce the ECa map with the ArcMap module of

Soil Apparent Electrical Conductivity (EC a ) and Altimetric Surveys
With the aim of measuring EC a , a contact-type sensor (Veris Technologies, Salina, KS, USA) was utilised. Measurements at a depth of 0-0.30 m were performed in April 2018. An all-terrain vehicle was used to pull the sensor. The average speed of the vehicle was 2.0 m s −1 ; consecutive passages, spaced 10 m, were made across the field. The spatial resolution of the EC a measurements was a 2 by 10 m grid, since a measurement was taken every second. A global navigation satellite system (GNSS) antenna was installed near the sensor. The obtained data were used to produce the EC a map with the ArcMap module of ArcGIS 9.3 software (v10.5, ESRI, Inc., Redlands, CA, USA), after conducting a geostatistical analysis with the extension Geostatistical Analyst. The data of the GNSS antenna were used to create the digital surface elevation model (elevation map) using the linear interpolation TIN tool from ArcGIS 9.3 and converted to a grid surface with a 1 m grid resolution.

Soil Sampling and Laboratory Reference Analysis
In the 24 sampling areas (Sentinel-2 pixels; Figure 3a), after measuring EC a , composite soil samples (comprised of five subsamples) were collected at a depth of 0-0.30 m. These soil samples were analysed for moisture content (SMC), particle size distribution (texture: sand, silt and clay content), pH, organic matter (OM) and cationic exchange capacity (CEC). The standard processes used in the laboratory were described in detail by Serrano et al. [5,6]. All maps of the soil parameters were produced after conducting a geostatistical analysis with ArcGIS using a 1 m grid resolution. Was used the inverse distance weighting (IDW) interpolation of the georeferenced data.

Cone Index Measurements
An electronic cone penetrometer "FieldScout SC 900" (Spectrum Technologies, Aurora, IL, USA) was used to measure the soil resistance to penetration (Cone Index, CI, in kPa) [15]. The main rules for the determination of CI values are standardised by the American Society of Agricultural and Biological Engineers (ASABE; EP542 and S313.3) [15].
In each sampling point, five CI measurements were carried out between 0-0.45 m (maximum depth allowed by the device), one in the central point of the sampling area, and one in each of its four quadrants. As suggested in other works [15], to minimise possible errors resulting from the uncertainty of maintaining a constant penetration rate during the determination, measurements were always carried out by the same experienced operator. When the insertion speed changes, the equipment registers an error, and the measurement is repeated. CI measurements were carried out in the 24  After the field measurements, data processing was carried out. A preliminary analysis was conducted to remove outliers from the data set. This procedure is fundamental, since the CI is measured using portable penetrometers, with manual operation, and the roughness of the soil surface and the variation in the speed of the rod going into the soil profile can influence the results [19]. The inconsistent and unreliable readings that may occur near the soil surface due the unevenness of the soil surface, led us not to consider the readings obtained in each point at 0 m depth, an aspect also suggested by Mayerfeld et al. [2]. The mean CI value of the set of five measurements was calculated for each sampling area and each depth of determination.
Taking into account, on the one hand, the recommendations of Mayerfeld et al. [2] that soil compaction investigations in silvopastures should extend to at least a depth of 0.30 m and, on the other hand, of Pentos et al. [14] that soil compaction measured in deeper soil layers is of no practical relevance because of limitations in rooting depth, in this study, the mean CI of 0-0.10 m, 0.10-0.20 m, 0.20-0.30 m and 0-0.30 m was calculated. In shallow soils, as is the case for soils typical of the Montado ecosystem [33,34], measurements below 0.30-0.35 m may be in contact with the bedrock, as reported by Mayerfeld et al. [2].

Animal Tracking with GPS Collars and Data Analysis
To monitor the grazing patterns, five randomly selected cows were fitted with GNSS (Global Navigation Satellite System) position loggers ("Digitanimal", Madrid, Spain). The tracking system consisted of a GNSS unit, a lithium battery pack, a PVC enclosure resistant to water and dust and a communication module (GPS collars) [35]. A total of four loggers were programmed to collect geolocation data every thirty minutes between 1 January and 17 March 2021, and the fifth receiver was programmed to collect geolocation data every five minutes between 6 and 19 May 2021. Data were transmitted over the "Sigfox", a global network dedicated to the internet of things featuring low power, a long range, and small data. The devices, weighing 265 g, were adjusted to the neck of the animals using a stripe with a buckle, without affecting the animals' movements. Figure 5 shows two patterns of animal behaviour throughout the year in the Montado ecosystem: UTC at peak summer sunshine hours (a) and in preferential grazing areas in other seasons (b). loggers were programmed to collect geolocation data every thirty minutes between 1 January and 17 March 2021, and the fifth receiver was programmed to collect geolocation data every five minutes between 6 and 19 May 2021. Data were transmitted over the "Sigfox", a global network dedicated to the internet of things featuring low power, a long range, and small data. The devices, weighing 265 g, were adjusted to the neck of the animals using a stripe with a buckle, without affecting the animals' movements. Figure 5 shows two patterns of animal behaviour throughout the year in the Montado ecosystem: UTC at peak summer sunshine hours (a) and in preferential grazing areas in other seasons (b). The geostatistical analyses of the GPS collars were carried out with the ArcGIS Desktop software (v10.5, ESRI, Inc., Redlands, CA, USA). The "Optimised Outlier Analysis (Spatial Statistics)" algorithm, based on incident points, obtained by the GPS collars, creates a map of statistically significant hot spots, cold spots and spatial outliers using the "Anselin Local Moran's I statistic". This map, with s 5 m spatial resolution, includes 5 classes ("Not significant", "High-High cluster", "High-Low outlier", "Low-High outlier", and "Low-Low cluster") and serves to characterise the grazing density pattern. With this tool, statistically significant spatial clusters of high values (hot spots), low values (cold spots) and outliers were identified within the dataset. The characteristics of the input feature class of the data to establish settings that produce optimal clusters were evaluated, and, automatically, (i) incident data were aggregated, (ii) multiple test and spatial dependence were corrected and (iii) a proper scale of analysis was determined. When a high The geostatistical analyses of the GPS collars were carried out with the ArcGIS Desktop software (v10.5, ESRI, Inc., Redlands, CA, USA). The "Optimised Outlier Analysis (Spatial Statistics)" algorithm, based on incident points, obtained by the GPS collars, creates a map of statistically significant hot spots, cold spots and spatial outliers using the "Anselin Local Moran's I statistic". This map, with s 5 m spatial resolution, includes 5 classes ("Not significant", "High-High cluster", "High-Low outlier", "Low-High outlier", and "Low-Low cluster") and serves to characterise the grazing density pattern. With this tool, statistically significant spatial clusters of high values (hot spots), low values (cold spots) and outliers were identified within the dataset. The characteristics of the input feature class of the data to establish settings that produce optimal clusters were evaluated, and, automatically, (i) incident data were aggregated, (ii) multiple test and spatial dependence were corrected and (iii) a proper scale of analysis was determined. When a high positive z-score for a given feature is obtained, the features of nearby areas have similar values. The "Output Feature Class" is "High-High" or, conversely, "Low-Low", respectively, for a statistically significant cluster of high or low values. On the other hand, the "Output Feature Class" is "High-Low" or "Low-High", respectively, when the feature has a high value and nearby areas have low values or, conversely, when the feature has a low value and nearby areas have high values.

Vegetation Multispectral Measurement and NDVI Time Series Reconstruction
For this study, a multi-temporal Sentinel-2 imagery data set, free of clouds and atmospherically corrected, was downloaded from the Copernicus data hub. Band 8 (B8; NIR; 842 nm) and band 4 (B4; RED; 665 nm), both with a 10 m spatial resolution, were used to calculate the satellite vegetation index (NDVI; Equation (1)) and for the reconstruction of the mean NDVI trends (NDVI time series records). Values are the mean of the set of pixel sampling areas corresponding to "high livestock trampling" and of the set of pixel sampling areas corresponding to "low livestock trampling".

Statistical Analysis
Descriptive statistical analysis was performed for all the evaluated soil parameters. Regression analysis with a 95% significance level (p < 0.05) and the analysis of variance (ANOVA) of the data were carried out using IBM SPSS Statistics package for Windows (version 28.0, IBM Corp., Armonk, NY, USA). The multiple comparisons (Tukey's HSD test) were applied for mean separation whenever the variables presented significant differences in the ANOVA. The specific analysis of the GPS collars and the determination of the NDVI from satellite image data were described above.

Soil Parameters: Spatial Variability and Relationship with Soil Apparent Electrical Conductivity
The descriptive statistics of the soil parameters of the experimental field are shown in Table 1. The low pH (5.5 ± 0.2), the low clay content (10.5 ± 1.8%) and OM (1.5 ± 0.3%) and the spatial variability are the most important features. According to the classification proposed by Pias et al. [19], the soil parameters show low spatial variability (pH and sand, with a CV < 12%) or medium-to-high variability (all other evaluated soil parameters, with a CV between 12 and 62%). This spatial variability is the basis for site-specific management [36][37][38], which, in this case of extensive animal production, corresponds to variable livestock management [5]. This spatial variability is also shown in the maps of Figure 6.   Figure 7 shows the elevation (a) and ECa (b) maps. These highlight the slightly undulating topography, characteristic of this region, and the very low ECa values (<5 mS.m; Figure 7b), typical of coarse-textured and dryland soils [24,35].  Figure 7 shows the elevation (a) and EC a (b) maps. These highlight the slightly undulating topography, characteristic of this region, and the very low EC a values (<5 mS.m; Figure 7b), typical of coarse-textured and dryland soils [24,35]. The relationship between soil apparent electrical conductivity (ECa) and the soil parameters of the experimental field (Table 2) showed positive and significant correlations with moisture (SMC) and clay and soil silt content, and negative and significant correlation with soil sand content. The reverse behaviour of clay and sand in the relationship with ECa is visible in Figure 8. Figure 9 shows the significant relationship between SMC and ECa and between SMC and CI: on the one hand, the positive effect of the SMC in ECa and, on the other, the decrease in the CI with the increase in the SMC are visible. The positive and significant correlations of ECa with the SMC, clay and soil silt content, and negative and significant correlation with soil sand content has also been verified in other works [14,24,[39][40][41] and shows the interest of the measurement of ECa as an expedient tool for identifying homogeneous management zones [22]. The relationship between soil apparent electrical conductivity (EC a ) and the soil parameters of the experimental field (Table 2) showed positive and significant correlations with moisture (SMC) and clay and soil silt content, and negative and significant correlation with soil sand content. The reverse behaviour of clay and sand in the relationship with EC a is visible in Figure 8. Figure 9 shows the significant relationship between SMC and EC a and between SMC and CI: on the one hand, the positive effect of the SMC in EC a and, on the other, the decrease in the CI with the increase in the SMC are visible. The positive and significant correlations of EC a with the SMC, clay and soil silt content, and negative and significant correlation with soil sand content has also been verified in other works [14,24,[39][40][41] and shows the interest of the measurement of EC a as an expedient tool for identifying homogeneous management zones [22].   In terms of possible impacts on soil compaction [14,42], the results of this study highlight the important variability in OM (1 to 2%) and clay (7 to 14%) contents. In this regard, several studies state that soil compaction is affected by small changes in soil texture [15][16][17]. Mayerfeld et al. [2] and Nawaz et al. [17] conclude, for example, that silt loam soils with low colloid contents are more susceptible than medium or fine-textured loamy and clay soils at low water contents, while sandy soils are slightly susceptible to soil compaction.
According to Pentos et al. [14], many physical, chemical and biological properties of soils that affect soil compaction, namely SMC, OM or particle size distribution, also influence the soil's electrical parameters, and, therefore, there is a close relation between soil compaction and soil apparent electrical conductivity. The EC a maps frequently show a high degree of correlation with soil compaction and, thus, can potentially provide a rapid alternative for assessing soil compaction [14,16].

Cone Index (CI) of Spatial Variability
The assessment of soil compaction through the CI in spring 2018 (Figure 10a) and spring 2019 (Figure 10b) showed different vertical profiles in the two subplots under study ("Field A" and "Field B"): while in 2018 there was a trend towards higher compaction in "Field B" and only in the 0-0.20 m soil layer, in 2019, there was an inverse trend, with higher compaction in "Field A" throughout the assessed soil profile (0-0.30 m).
In terms of possible impacts on soil compaction [14,42], the results of this study highlight the important variability in OM (1 to 2%) and clay (7 to 14%) contents. In this regard, several studies state that soil compaction is affected by small changes in soil texture [15][16][17]. Mayerfeld et al. [2] and Nawaz et al. [17] conclude, for example, that silt loam soils with low colloid contents are more susceptible than medium or fine-textured loamy and clay soils at low water contents, while sandy soils are slightly susceptible to soil compaction.
According to Pentos et al. [14], many physical, chemical and biological properties of soils that affect soil compaction, namely SMC, OM or particle size distribution, also influence the soil's electrical parameters, and, therefore, there is a close relation between soil compaction and soil apparent electrical conductivity. The ECa maps frequently show a high degree of correlation with soil compaction and, thus, can potentially provide a rapid alternative for assessing soil compaction [14,16].

Cone Index (CI) of Spatial Variability
The assessment of soil compaction through the CI in spring 2018 (Figure 10a) and spring 2019 (Figure 10b) showed different vertical profiles in the two subplots under study ("Field A" and "Field B"): while in 2018 there was a trend towards higher compaction in "Field B" and only in the 0-0.20 m soil layer, in 2019, there was an inverse trend, with higher compaction in "Field A" throughout the assessed soil profile (0-0.30 m). It can also be seen that the mean CI in "Field A" showed a similar pattern in 2018 and 2019, not exceeding 2000 kPa, while in "Field B", the CI reached values above 2500 kPa in the 5-10 cm soil layer in 2018, not exceeding 1500 kPa in 2019. It is important to link this evaluation to the livestock management: the soil compaction measurement in 2018 was carried out near the end of April, after all the animals had been in "Field B" since the beginning of April; on the other hand, the soil compaction measurement in 2019 was It can also be seen that the mean CI in "Field A" showed a similar pattern in 2018 and 2019, not exceeding 2000 kPa, while in "Field B", the CI reached values above 2500 kPa in the 5-10 cm soil layer in 2018, not exceeding 1500 kPa in 2019. It is important to link this evaluation to the livestock management: the soil compaction measurement in 2018 was carried out near the end of April, after all the animals had been in "Field B" since the beginning of April; on the other hand, the soil compaction measurement in 2019 was carried out near the end of March, after all the animals had been in "Field A" between January and March. Figures 11 and 12 show the spatial variability in the CI in the two assessments carried out (2018 and 2019, respectively) in the various soil layers (0-10 cm; 10-20 cm; 20-30 cm; and 0-30 cm). These figures reflect the CV of 20-40% shown in Table 1. The spatial pattern is variable both in terms of depth and between dates.
The descriptive statistics of the SMC and CI relative to all data obtained between September 2021 and March 2022 are presented in Table 3. The results of the ANOVA applied to the CI measurements are also presented. This analysis shows significant differences for the variables "Date of measurement" and "Depth", and non-significant differences for the variables "Tree canopy" and "Fields". The compaction profiles show similar patterns for OTC and UTC for all of the four evaluation dates (Figure 13). The mean separation (multiple comparison of Tukey) for the variable "Date of measurement" showed higher CI values in September 2021, followed by March 2022, with no significant differences between November and December 2021. These results show a trend towards lower CI values as the SMC increases (Figure 14), which confirms the already presented relationship between the SMC and CI (Figure 9). The mean separation for the variable "Depth" showed higher CI values at a depth of 0.10-0.20 m, relative to the other two soil layers considered (0-0.10 m and 0.20-0.30 m). carried out near the end of March, after all the animals had been in "Field A" between January and March. Figures 11 and 12 show the spatial variability in the CI in the two assessments carried out (2018 and 2019, respectively) in the various soil layers (0-10 cm; 10-20 cm; 20-30 cm; and 0-30 cm). These figures reflect the CV of 20-40% shown in Table 1. The spatial pattern is variable both in terms of depth and between dates.  The descriptive statistics of the SMC and CI relative to all data obtained between September 2021 and March 2022 are presented in Table 3. The results of the ANOVA applied to the CI measurements are also presented. This analysis shows significant     [19]. The mean CI of 1900-2000 kPa registered in this study in 2018 and 1300-1700 kPa registered in 2019 (Table 1) are values below those limits generally considered critical for plant root growth. In the evaluations performed in 2021/2022, however, the mean CI was 2200-2300 kPa in November and December 2021, reaching very high mean values in March 2022 (>2700 kPa) and, especially, in September 2021 (>3100 kPa), showing an inverse relationship with SMC. This aspect is very important, since the vegetative cycle of dryland pastures normally begins in September (after the first autumn rains), an important phase in plant root development [42]. This trend towards higher CI values as SMC decreases, or vice versa (exponential relationship), confirms the results of other studies [15,18] and is attributed to the reduction in the cohesive forces between clay particles [43]. It should be noted, however, that this relationship between the CI and SMC restricts direct comparisons of CI values among same soils with different moisture contents [15].
In this study, in addition to the significant differences found when comparing different dates of the CI measurements, resulting basically from the change in SMC, but also from changes in the animal grazing management (an aspect discussed in the next point), it is important to note the absence of significant differences in the CI between OTC and UTC regarding the topsoil layer (0-0.30 m). Several studies report higher compaction in UTC areas [1,6]. As trees grow, their aboveground weight is transferred to the soil through surface roots, which also exert compression forces on near-by soil as they increase in diameter during radial growth [1]. However, this effect can be mitigated by the accumulation of leaves and other residues UTC [1], leading to higher levels of OM in UTC areas [34], which tends to reduce soil compaction because these open new soil channels  [19]. The mean CI of 1900-2000 kPa registered in this study in 2018 and 1300-1700 kPa registered in 2019 (Table 1) are values below those limits generally considered critical for plant root growth. In the evaluations performed in 2021/2022, however, the mean CI was 2200-2300 kPa in November and December 2021, reaching very high mean values in March 2022 (>2700 kPa) and, especially, in September 2021 (>3100 kPa), showing an inverse relationship with SMC. This aspect is very important, since the vegetative cycle of dryland pastures normally begins in September (after the first autumn rains), an important phase in plant root development [42]. This trend towards higher CI values as SMC decreases, or vice versa (exponential relationship), confirms the results of other studies [15,18] and is attributed to the reduction in the cohesive forces between clay particles [43]. It should be noted, however, that this relationship between the CI and SMC restricts direct comparisons of CI values among same soils with different moisture contents [15].
In this study, in addition to the significant differences found when comparing different dates of the CI measurements, resulting basically from the change in SMC, but also from changes in the animal grazing management (an aspect discussed in the next point), it is important to note the absence of significant differences in the CI between OTC and UTC regarding the topsoil layer (0-0.30 m). Several studies report higher compaction in UTC areas [1,6]. As trees grow, their aboveground weight is transferred to the soil through surface roots, which also exert compression forces on near-by soil as they increase in diameter during radial growth [1]. However, this effect can be mitigated by the accumulation of leaves and other residues UTC [1], leading to higher levels of OM in UTC areas [34], which tends to reduce soil compaction because these open new soil channels contribute nutrients that support the soil rhizosphere [1], increasing the resistance to soil deformation by increasing elasticity [17].

Effect of Livestock Trampling on Soil Compaction
Another aspect evaluated in this study was the effect of livestock trampling on soil compaction. Figure 15 show the grazing density map based on georeferenced information obtained by the GPS collars. In these, five classes were used: "Not significant", "High-High cluster", "High-Low outlier", "Low-High outlier" and "Low-Low cluster". , 23, x FOR PEER REVIEW 19 of 29 contribute nutrients that support the soil rhizosphere [1], increasing the resistance to soil deformation by increasing elasticity [17].

Effect of Livestock Trampling on Soil Compaction
Another aspect evaluated in this study was the effect of livestock trampling on soil compaction. Figure 15 show the grazing density map based on georeferenced information obtained by the GPS collars. In these, five classes were used: "Not significant", "High-High cluster", "High-Low outlier", "Low-High outlier" and "Low-Low cluster". For the statistical analysis, only two classes were considered: the "High-High cluster" (a statistically significant cluster of high values), which includes the sampling points A1, A5 and A9, and the "Low-Low cluster" (a statistically significant cluster of low values), which includes the sampling points A2, B1, B2, B4, B10 and B11. The other sampling points (A3, A6 and B8) were classified as "Not Significant". "High-High" areas correspond to a For the statistical analysis, only two classes were considered: the "High-High cluster" (a statistically significant cluster of high values), which includes the sampling points A1, A5 and A9, and the "Low-Low cluster" (a statistically significant cluster of low values),  Our results show a systematic trend of a greater CI in areas with "high" livestock trampling, compared to areas with "low" livestock trampling. However, the evaluation of soil compactness of silvopastoral systems (as is the Montado) is complex due to the difficulty in determining the relative importance of trees versus livestock trampling on the soil's physical properties [1]. This work shows that the compaction resulting from animal trampling was significant in OTC areas (in all evaluated dates and depths); nevertheless, in UTC areas, the effect of animal trampling was significant only in the 0-0.10 m soil layer. The tree canopy potentially provides a cushion between hoof and soil during grazing, which can mitigate the compaction by animal trampling [1]. These results seem to indicate the important role played by tree roots in the structural support of the soil in layers below a depth of 0.10 m, reducing the potential compaction effect that animal trampling tends to produce. On the other hand, a soil with higher levels of OM (as is the case of UTC soil [34]) generally has a better structure because OM helps create large, strong soil aggregates that help resist compaction [8]. In addition, this complexity also reflects the balance of restorative and compaction processes that may occur in a differentiated manner UTC and OTC throughout the vegetative cycle of the pasture. Besides the direct effect of trampling, livestock can indirectly change soil properties by consuming vegetation that would otherwise contribute to the availability of organic matter to soil microfauna by reducing the amount and extent of fine roots that open new soil channels for the soil rhizosphere [1]. The biological restorative action of soil decompaction provided by pasture roots and by the activity of soil mesofauna is enhanced when the pasture management favours the accumulation of phytomass in the aerial part and in the root system of the plants [42,45], which is what happened in "Field B" of this experiment (several months without grazing animals and without any sampling area with "high" livestock trampling). Potentially, rest periods from grazing also enhance soil recovery by physical processes (variations in soil temperature and moisture content) [1,9].
Another relevant aspect relates to the depth at which soil compaction occurs due to livestock trampling. The compaction was significant at a depth of 0-0.10 m on three of the four dates of evaluation in both OTC and UTC areas ( Table 6). At depths of 0.10-0.20 m and 0.20-030 m, the effect of livestock trampling was significant only in OTC areas. Several studies report that the highest compaction impact caused by animal trampling is confined to the topsoil layer under wet soil conditions [8], which has practical implications because deeper soil layers (below 0.15 m [17]) are generally slower to recover from compaction [1,17]. For example, Roesh et al. [8] and Debiasi [20] report that the greatest impact by livestock hooves is limited to the top 0-0.05 m soil layer, while Sharrow et al.
[1], Reichert et al. [42] and Vzzotto et al. [45] extend this influence to the top 0-0.10 m soil layer, and Mayerfeld et al. [2], Medina [12] and Donkor et al. [13] indicate that this impact can reach up to a depth of 0-0.15 m or even 0.20 m [17]. This effect depends on several factors, among them the animal load [20] or grazing intensity [2]. When the animal load is very high, an increase in compaction due to trampling may occur in deeper layers [20]. The magnitude of the differences in the CI in the top 0.20 m also appears to be variable over time; therefore, the longer-term tracking of changes in the CI may be required [2].

Effect of Livestock Trampling on Pastures
The impact of livestock trampling on pastures was evaluated by the NDVI time series obtained from satellite imagery (Sentinel-2). Figure 17 shows the reconstruction of the mean NDVI time series retrieved between July 2021 and June 2022, without the records affected by the existence of clouds. The values are the mean of the set of pixel sampling areas of "high livestock trampling" and of the set of pixel sampling areas of "low livestock trampling". It is possible to observe higher NDVI values in March, which reflects the differential grazing management ("Field A" versus "Field B"). In the period of peak pasture production (April and May), the behaviour of the NDVI is inverted, with higher values in areas of "low livestock trampling", which are potentially less compacted and, therefore, with greater vegetative vigour.
affected by the existence of clouds. The values are the mean of the set of pixel sampling areas of "high livestock trampling" and of the set of pixel sampling areas of "low livestock trampling". It is possible to observe higher NDVI values in March, which reflects the differential grazing management ("Field A" versus "Field B"). In the period of peak pasture production (April and May), the behaviour of the NDVI is inverted, with higher values in areas of "low livestock trampling", which are potentially less compacted and, therefore, with greater vegetative vigour. The composition and production of a range of vegetation are the primary vulnerable factors to grazing pressure [46]. In this study, the impact of livestock trampling on pastures was evaluated by the NDVI time series obtained from Sentinel-2 imagery. Although this approach has some limitations, namely, the constraint resulting from poor-quality images on cloudy days, which can occur frequently during the pasture's vegetative cycle [28,29], the NDVI can be used to monitor the pasture's development status since it mainly reflects the chlorophyll content, an indicator of pastures' vegetative vigour [29,30].
The results of this study show ( Figure 17) an NDVI pattern that is typical of the vegetative cycle of dryland pastures in the Mediterranean region, already presented in other works [28][29][30]. This pattern is similar in "high" and "low" livestock trampling areas and reflects the effect of the evolution of air temperature and the distribution of precipitation [30]. In Figure 17, two moments are identified where the NDVI patterns indicate significantly different behaviours in the two situations of livestock trampling ("high" and "low"): (i) "differential grazing management" in "Field A" and "Field B" (stocking rates); (ii) "pasture peak production". Higher NDVI values were found in March in "high The composition and production of a range of vegetation are the primary vulnerable factors to grazing pressure [46]. In this study, the impact of livestock trampling on pastures was evaluated by the NDVI time series obtained from Sentinel-2 imagery. Although this approach has some limitations, namely, the constraint resulting from poor-quality images on cloudy days, which can occur frequently during the pasture's vegetative cycle [28,29], the NDVI can be used to monitor the pasture's development status since it mainly reflects the chlorophyll content, an indicator of pastures' vegetative vigour [29,30].
The results of this study show ( Figure 17) an NDVI pattern that is typical of the vegetative cycle of dryland pastures in the Mediterranean region, already presented in other works [28][29][30]. This pattern is similar in "high" and "low" livestock trampling areas and reflects the effect of the evolution of air temperature and the distribution of precipitation [30]. In Figure 17, two moments are identified where the NDVI patterns indicate significantly different behaviours in the two situations of livestock trampling ("high" and "low"): (i) "differential grazing management" in "Field A" and "Field B" (stocking rates); (ii) "pasture peak production". Higher NDVI values were found in March in "high livestock trampling" areas, which reflects the differential grazing management ("Field A" versus "Field B"): during this month, grazing was concentrated in "Field B", where "low livestock trampling" areas predominate, while "Field A", where "high livestock trampling" areas predominate, was left to rest (without grazing). The fact that in March the pasture of "Field A" was not grazed, accompanied by the rise in air temperature and an important accumulated rainfall (135 mm; Figure 4), allowed the recovery of vegetative vigour, which translated into a significant increase in the NDVI. In the period of pasture peak production (April and May), the behaviour of the NDVI is inverted (Figure 17), with higher values in areas of "low livestock trampling", which are potentially less compacted and, therefore, with greater vegetative vigour. Gao et al. [47] and Jin et al. [48] also state that climatic factors such as precipitation and temperature are responsible for inter-annual fluctuations in biomass and vegetative vigour. Psyllos et al. [3] emphasise the greater yield potential of rotational grazing between fields (grazed versus not grazed). This is, however, a complex issue, as these pastures are biodiverse, comprised of grasses, legumes, forbs and other species. Different botanical species mature at different rates [49] and show different susceptibilities to animal trampling. For example, animal trampling has a positive effect on the germination of some grass and forbs species [46]. On the other hand, certain botanical species with prostrate growth (such as legumes), are more tolerant than others (grasses, for example) to livestock trampling [50], i.e., show differential susceptibility to damage by animal trampling [49], which can result in a greater abundance of trampling-tolerant plants, affecting the composition of an herbaceous plant community [50]. The assessment of the impact of animal trampling on the floristic composition of pastures is a subject that requires future long-term investigations.

Mitigation of the Negative Effects of Livestock Trampling
The various studies carried out and published on the effects of livestock trampling are unanimous on the greater susceptibility to soil compaction and pasture damage during periods of high SMC. Multiple management strategies to avoid or to mitigate the negative impact of livestock trampling should be a priority [42]. According to Reichert et al. [42], two preventative measures may be implemented: (i) the use of mobile fences to impede the entry of animals into more susceptible areas in the days following the occurrence of major precipitation events; (ii) control over grazing management (stocking rates) to permanently guarantee a minimum height of pasture to ensure soil surface protection.
One of the main causes of soil and pasture degradation resulting from livestock trampling is the natural tendency of animals to agglomerate in certain areas of the pasture. Thus, an understanding of the potential environmental effects of the concentration areas is necessary to adequately consider mitigating grazing management practices [51]. Some examples of grazing management practices that most influence grazing distribution and, consequently, livestock trampling, are the location of watering sites, supplement feeders (concentrate, hay or minerals), tree shade or fence and gate placement [51]. Some compaction at gates and waterers, for example, is inevitable, but can be minimised. The traditional advice of moving waterers and feeders, fences and gates to limit livestock concentration areas seems appropriate [51]. When it is possible to identify and circumscribe these critical compaction sites, then specific improvement interventions can be carried out at only the affected areas, minimising the cost and time necessary for the operation [42]. Furthermore, temporary livestock exclusion has been shown in previous studies to be an attractive tool to improve the physical properties of soil [9]. This site-specific management fits into the perspective of Precision Agriculture or, in this case, of Precision Grazing, where several technological tools can make an excellent contribution to the definition of homogeneous management zones, the basis for implementing differentiated management strategies [22,38]. The technologies used in this study (e.g., GPS collars, cone penetrometer, EC a sensor and satellite imagery) can be practical tools that allow us to monitor animal grazing patterns, the spatial and temporal variability of the physical properties of soil or the condition of pastures. All these are important factors that can help farmers to assess potential compaction due to animal trampling and to adopt sustainable livestock management systems [8].
The global technological approach proposed in this study to correlate the soil compaction variables related to the effect of animal trampling is shown in Figure 18. Future research on soil compaction in the Montado ecosystem, in addition to animal trampling, may include the monitoring of the effect of tree roots. Here, too, new technological approaches can be combined. One study by Xu et al. [52] is an example of this endeavour based on the development and validation of a new tool (a high-fidelity 3D root system morphological model) for simulating the mechanical analysis of root-soil composites. paction variables related to the effect of animal trampling is shown in Figure 18. Future research on soil compaction in the Montado ecosystem, in addition to animal trampling, may include the monitoring of the effect of tree roots. Here, too, new technological approaches can be combined. One study by Xu et al. [52] is an example of this endeavour based on the development and validation of a new tool (a high-fidelity 3D root system morphological model) for simulating the mechanical analysis of root-soil composites.

Conclusions
The economic and environmental sustainability of extensive livestock production systems requires the optimisation of soil management, pasture production and animal grazing. All these aspects are interdependent and linked to soil compaction. The compaction resulting from livestock trampling was significant in areas OTC in the three soil layers considered (0-0.10 m; 0.10-0.20 m; 0.20-0.30 m), but in areas UTC this effect was only significant in the 0-0.10 m soil layer. These results suggest that this is a dynamic process throughout the year, with recovery cycles associated with grazing management, seasonal fluctuations in soil moisture and temperature, or the spatial variation in specific soil characteristics (namely clay and organic matter contents). The application of technologies such as those used in this study (e.g., GPS collars, cone penetrometer, ECa sensor and satellite imagery) can be practical tools for monitoring animal grazing patterns, the spatial and temporal variability in the physical properties of the soil and the condition of pastures. This approach shows the opportunity provided by Precision Agriculture technologies, such as proximal and remote sensing, to generate knowledge, support decision making and respond to the challenge of a holistic and sustainable management system of the Montado ecosystem.

Conclusions
The economic and environmental sustainability of extensive livestock production systems requires the optimisation of soil management, pasture production and animal grazing. All these aspects are interdependent and linked to soil compaction. The compaction resulting from livestock trampling was significant in areas OTC in the three soil layers considered (0-0.10 m; 0.10-0.20 m; 0.20-0.30 m), but in areas UTC this effect was only significant in the 0-0.10 m soil layer. These results suggest that this is a dynamic process throughout the year, with recovery cycles associated with grazing management, seasonal fluctuations in soil moisture and temperature, or the spatial variation in specific soil characteristics (namely clay and organic matter contents). The application of technologies such as those used in this study (e.g., GPS collars, cone penetrometer, EC a sensor and satellite imagery) can be practical tools for monitoring animal grazing patterns, the spatial and temporal variability in the physical properties of the soil and the condition of pastures. This approach shows the opportunity provided by Precision Agriculture technologies, such as proximal and remote sensing, to generate knowledge, support decision making and respond to the challenge of a holistic and sustainable management system of the Montado ecosystem.