## 2 Materials and Methods
### Data Processing and Statistical Analysis  
As a first step, the data were summarized in Microsoft Excel sheets by exporting them from the penetrometer to a personal computer using the dedicated software. Then, the data were arranged based on the project planned with the software, and then the measurements taken in the field were labeled with codes to indicate the extraction means and soil disturbance category to which they belonged. As such, the measurement datasets were grouped based on the extraction means used (H, F, S), type of data (control or pass number, hereafter C or P) and the number of the sample point (1 to 10). Table 2 summarizes the features of the initial datasets. Figure 2: Main features of the testing areas: CP1 to CP10 stand for control sampling points that had the same location in all of the \(n\) testing areas, where \(n=3\) (H, F or S); SP1 to SP10 stand for sampling points on the skid trails and had the same location in all of the \(n\) testing areas, where \(n=3\) (H, F or S); g stands for a physical ground guiding mark; DT stands for the distance between the centers of the wheels’ footprints which depended on the extraction means used and increased in the order H, S, F; m stands for the means used for extraction, m = H, F or S; p stands for the number of passes counted until a measurement on penetration resistance was taken, \(p=1,2,3,4,5,10,15\). To be able to detect natural changepoints in the control datasets (HC, FC and SC), while accounting for the variability in penetration resistance data at the same depth, subsets of 10 observations each (CP1 to CP10) taken at each centimeter in depth (from i = 0 to 65) were used to estimate the mean (hereafter x), standard deviation (hereafter s) and coefficient of variation (hereafter cv) of each centimeter in soil depth. Variation in values of x, s and cv was then plotted against the soil depth for the three control datasets, along with the absolute differences between For each of the studied extraction means, based on the data subsets standing for the mean values in control (C) and compacted soil, absolute differences were calculated to compare visually the changes induced in the soil penetration resistance following passes P1, P2, P3, P4, P5, P10 and P15. In a similar way, changes induced by passes were visually compared by plotting the absolute differences in penetration resistance between P1-C, P3-P2 and so on, up to P15-P10. Based on the graphical representation of the control data (x, s and cv) against soil depth (0 to 65 cm), segments of relatively uniform behavior in penetration resistance were identified visually and used to divide the full soil depth into five characteristic layers: from 0 to 10, from 11 to 20, from 21 to 30, from 31 to 50 and from 51 to 65 cm in depth. Then, all the point observations of these layers were grouped based on the type of data and extraction means, forming the subsets for layer-level data comparison. The statistical workflow used to compare the data included a check for normality in data which was performed by a Shapiro-Wilk test [48] to be able to select the most appropriate type of statistical comparison test. However, most of the layer-level datasets failed the normality assumption (Tables A2-A4), hence a Mann-Whitney [49] comparison test was used to check if there were significant differences between the control and pass-level datasets. The Mann-Whitney test assumes that (i) there is a dependent variable measured on a continuous or ordinal scale, (ii) there is one independent variable that consists of two categorical groups and (iii) there is independence in observations [50]. In addition, the test requires a check for the distribution of scores to determine if the distributions of compared pairs have the same or a different shape [50]. As such, the test