Large‐scale detection and quantification of harmful soil compaction in a post‐mining landscape using multi‐configuration electromagnetic induction

Fast and accurate large‐scale localization and quantification of harmfully compacted soils in recultivated post‐mining landscapes are of particular importance for mining companies and the following farmers. The use of heavy machinery during recultivation imposes soil stress and can cause irreversible subsoil compaction limiting crop growth in the long term. To overcome or guide classical point‐scale methods to determine compaction, fast methods covering large areas are required. In our study, a recultivated field of the Garzweiler mine in North Rhine‐Westphalia, Germany, with known variability in crop performance was intensively studied using non‐invasive electromagnetic induction (EMI) and electrode‐based electrical resistivity tomography (ERT). Additionally, soil bulk density, volumetric soil water content and soil textures were analysed along two transects covering different compaction levels. The results showed that the measured EMI apparent electrical conductivity (ECa) along the transects was highly correlated (R2 > .7 for different dates and depths below 0.3 m) to subsoil bulk density. Finally, the correlations established along the transects were used to predict harmful subsoil compaction within the field, whereby a spatial probabilistic map of zones of harmful compaction was developed. In general, the results revealed the feasibility of using the EMI derived ECa to predict harmful compaction. They can be the basis for quick monitoring of the recultivation process and implementation of necessary melioration to return a well‐structured soil with good water and nutrient accessibility, and rooting depths for increased crop yields to the farmers.


| INTRODUCTION
Soil compaction has been considered as one of the main soil threats by the Food and Agriculture Organization of the United Nations (FAO & ITPS, 2015), and mapping soil compaction is a key issue in studying land degradation (Alaoui & Diserens, 2018). According to Batey and McKenzie (2006), the changes in soil properties induced by soil compaction can be grouped into primary and secondary effects. Primary effects include the alteration of soil physical properties such as bulk density, porosity and strength. Secondary effects are induced by the primary effects and include a reduced permeability of air and water (e.g. Delgado et al., 2007), an increase in resistance for root growth (e.g. Haigh & Sansom, 1999;Stone, 1988;Taylor et al., 1966), a reduction in nutrient uptake by plants (e.g. Batey & McKenzie, 2006) and increased surface run-off and erosion (e.g. Haigh & Sansom, 1999;Schack-Kirchner et al., 2007).
Soil compaction of agricultural land is caused by management operations (Hamza & Anderson, 2005), mainly by the use of heavy machinery or intensive grazing (Mulholland & Fullen, 1991), which is even more detrimental on wet soils (Batey, 2009). In opencast mining areas, where soil is reclaimed for agricultural practices, the use of heavy machinery endangers soil restoration by inducing local subsoil compaction. While topsoil bulk densities and compaction may vary in short-term because of seasonal land management (e.g. tillage, seedbed preparation, harvest) and atmospheric impact (e.g. rainfall and frost), subsoil compaction cannot be reversed by classical soil cultivation and may therefore persist over decades (Batey, 2009). As the recultivated land will be handed over to farmers after a certain 'recultivation period' of several years, and as the farmers do expect farmland without negative conditions, it is of uppermost importance to recognize soil compactions early at the restoration sites. This would allow initiation of soil management measures that ensure optimum future agricultural use.
Classical methods to assess soil compaction, apart from a visual morphological description of the soil structure and strength, are point-scale measurements of bulk density, penetration resistance, shear resistance, permeability of air and water and image analysis of thin sections (Batey & McKenzie, 2006). However, these are destructive, timeand cost-intensive as well as restricted to a few sampling points within the field or landscape, which might not be able to cover the wide lateral and vertical variation in the degree of compaction. In contrast, comparable low-cost non-destructive geophysical methods such as electromagnetic induction (EMI), electrical resistivity tomography (ERT) and ground-penetrating radar (GPR) can help to assess and understand the spatial variability of soil properties (e.g. Binley et al., 2015;Parsekian et al., 2015). ERT offers great potential to study small scale soil functions (e.g. Garré et al., 2013;Mary et al., 2020;Weigand & Kemna, 2017), while being restricted to static layouts of electrode arrays. Field-scale GPR measurements can be performed relatively quick with challenging data interpretation because of reflection, refraction, wave-traps and general 3D wave propagation patterns. Non-invasive EMI measurements offer great potential to measure areas of several hectares (e.g. Brogi et al., 2019;Corwin & Lesch, 2003;Doolittle & Brevik, 2014). By measuring apparent electrical conductivity (ECa), EMI demonstrated the potential to estimate soil texture variations (e.g. Anderson-Cook et al., 2002;Heil & Schmidhalter, 2012;Mertens et al., 2008), soil water content (e.g. Altdorff et al., 2017;Sheets & Hendrickx, 1995), and pore water conductivity (e.g. Dakak et al., 2017;Kaufmann et al., 2020;Rhoades & Corwin, 1990), but also for crop growth studies (e.g. Brogi et al., 2020;Stadler et al., 2015;von Hebel et al., 2018). For a more detailed overview on the influence of soil and soil conditions on the ECa measurements, we recommend the literature review in Al-Gaadi (2012).
Irrespectively of the larger number of studies analysing the soil properties listed above, only little attempt has been made to estimate soil compaction. However, soil electrical properties are influenced by bulk density changes. Compacted soils exhibit more soil particles per unit volume and a denser packing with increased contact between the particles, both leading to a better bulk electrical conductance (Brevik & Fenton, 2015;Sauer et al., 1955;Wyllie & Southwick, 1954). Compaction additionally changes the soil hydraulic properties because of changes in porosity and pore geometry (Dedousis & Bartzanas, 2010;Zhang et al., 2006), whereby more or less water can be held against external forces (e.g. drainage, evaporation, and/or root water uptake) and thereby the electrical properties change compared with uncompacted soil. Overall, the soil structural changes because of compaction imply electrical changes that should be detectable with EMI measurements.
In a study of Al-Gaadi (2012), EMI was employed at different heights above ground and different soil water contents of a sandy field plot to assess three different levels of soil compaction, induced by a vibrator plate. Reference measurements were made by a compaction meter. The results showed that EMI readings correlate with soil compaction. On-ground measurements obtained a R 2 of .90, while increasing the measurement height up to 0.4 m decreased the correlation to R 2 = .47. Additionally, lower soil water contents showed best correlation between EMI and compaction meter readings. Further, a study of Brevik and Fenton (2015) revealed a linear relationship between increasing bulk density and the corresponding increase in soil ECa within the ranges of bulk density changes that they observed on three different test sites. Here, it has to be noted that only the compaction in the upper soil layer (up to 0.25 m) was investigated. Sudduth et al. (2010) used EMI data to identify soil compaction at two different sites in central Missouri (USA) and observed a moderate correlation between the EMI-ECa and depth to a compacted claypan layer. Further, a medium correlation between maximum penetrometer cone index and ECa (R 2 = .48) was detected, whereby soil water content was used as an additional covariable in the regression. Irrespectively of the low correlation found, the authors stated that ECa can be useful to locate areas within the field, which are most likely to exhibit high levels of compaction.
Based on the sparse information found in literature that EMI can be used for the detection of compaction and the urgent need to detect subsoil compaction for recultivation purposes at an early stage of the recultivation process, a field study was performed in 2019 on a recultivated agricultural site with observed variability in crop performance. Over the course of five months, the field was intensively investigated using EMI and ERT along with ground truth soil measurements of soil water content, bulk density, and soil texture. These data were used to validate that the main indicators for soil compaction are linked to electrical conductivity changes and EMI offers a cheap methodology to recultivate land for optimal agricultural cropping conditions.

| Site description
The study was conducted on an agricultural field near the RWE Power AG opencast lignite mine Garzweiler, which is located approximately 35 km northwest of Cologne in North Rhine-Westphalia, Germany (51°04′11.5″N 6°29′34.9″E). The study area is located centrally within a 23.5 ha field with low elevation and a flat surface (90 to 92.5 m a.s.l.). Its climate is characterized by an average precipitation of 805 mm and a mean annual temperature of 10℃. The geological conditions of the Garzweiler opencast mine reveal sandy-gravelly Tertiary and Quaternary sediments as well as Saalian and Weichselian loess (Dumbeck, 2014). For the recultivated soil, the stackers spread a mixture of calcareous loess and the former solum of loess-derived Luvisols at the top of the refilled mining areas. After a certain settling time, the loess-solum material is levelled in a soil-sparing fashion using caterpillars. The final result is a recultivated soil with at least 2 m thickness that refers to a Calcaric Regosol according to WRB (IUSS Working Group WRB, 2015). The soils typically have a silt content of 70%-80%, a clay content of 10%-22%, and a sand content of 2%-9% and can be assigned to the soil textural class silt loam according to WRB (IUSS Working Group WRB, 2015), whereby the coarse fraction (>2 mm) of the soils and the organic carbon contents are very small.
The investigated field was recultivated in 2013 and has since been managed by the mining company RWE Power AG farm. For the first three years, alfalfa as a deep rooting plant was grown. After alfalfa, winter wheat was cropped, and in October 2018, winter barley was sown, which was harvested in July 2019. Soil management was performed using mouldboard ploughing and harrowing for seedbed preparation.

| Soil sampling and analysis
Within the 23.5 ha field, an area of 150 × 150 m was selected for intensive investigation. Two transects were defined, which cross structures that attracted attention by high ECa values during a coarse hand-held pre-EMI survey in January 2019 and by an observed low plant productivity during the last growing season (see aerial photograph Figure 1a). Measurements of electrical conductivity (EMI, ERT), and soil sampling for gravimetric and volumetric water content, dry bulk density and soil texture determination were conducted at different times during the growing season ( Figure 1c). As a full coverage EMI survey of the area was not feasible during the growing season, these measurements were performed after barley harvest. The same holds for the large undisturbed soil core sampling, which was also carried out after harvest.
For soil texture characterization, soil samples were collected by augering at eight locations along each transect (at transect position 0, 2.5, 5, 7.5, 10, 15, 22.5 and 30 m from start). The soil in the 1 m auger was divided into predefined depth intervals of 0-0.3, 0.3-0.6 and 0.6-0.9 m, whereby the 0-0.3 m depth corresponds to the plough horizon. Soil texture was analysed in the laboratory according to DIN ISO 11277 (2002) by wet sieving and sedimentation using the SEDIMAT 4-12 (Umwelt-Geräte-Technik GmbH). Here, it has to be noted that the sand fraction in DIN ISO 11277 (2002) is defined between 2 and 0.063 mm, according to IUSS Working Group WRB (2015).
Gravimetric water content was determined at two dates (8th of May and 27th of June) from auger samples at same depth intervals as soil texture (0-0.3, 0.3-0.6, and 0.6-0.9 m). The samples were transferred into plastic bags after sampling and stored at 4℃ until measurements. The gravimetric water contents were determined by the weight loss after oven drying at 105℃ for 3 days. Based on the knowledge of the dry bulk density, the gravimetric water contents were transformed to volumetric water contents.
As the measurement of the dry soil bulk density required undisturbed samples, large soil cores (0.1 m in diameter and 1 m length) were taken 17th of July after harvest using an automated caterpillar based corer. The cores were sampled in Plexiglas tubes, sealed, and stored until measurements at 4℃. In the laboratory, the tubes were cut into equal increments (0-0.3 m for topsoil (Aphorizon) and afterwards in 0.1 m increments) dried at 105℃ for 3 days and weighed after that. By the knowledge of the dry weight and sampling volume, the dry bulk density was calculated.

| Electromagnetic Induction (EMI) measurements and calibration
Frequency domain EMI systems use an alternating current with a fixed frequency to generate and transmit a primary magnetic field. This field induces eddy currents in the electrically conductive subsurface, which in turn generate a secondary magnetic field (Keller & Frischknecht, 1966). Both, the primary and secondary field, create currents in a receiver coil, where the ratio is related to the apparent electrical conductivity (ECa) of the subsurface (Ward & Hohmann, 1988). The measured ECa represents certain integral depth ranges that depend on the coil orientation and separation. The depth range at which 70% of the relative EMI response accumulates is defined as the depth range of investigation (DOI) (see Figure 2 for an example of our measurement setup). While vertical coplanar (VCP) orientations are most sensitive to the shallow subsurface with DOIs of 0.75 times the coil separation, horizontal coplanar (HCP) configurations are most sensitive at a depth of around 0.4 times the coil separation and a DOI of around 1.5 times the coil separation (McNeill, 1980). In this study, we simultaneously used a three-and a six-coil EMI device (GF instruments) with operating frequencies of 30 and 25 kHz, respectively, to measure VCP and HCP configurations at the same time. Thereby, the CMD-MiniExplorer, measuring in VCP orientation, has coil separations of 0.32, 0.71 and 1.18 m (VCP 0.32 m to VCP 1.18 m) and the custom-made CMD-MiniExplorer special edition, measuring in HCP orientation, has separations of 0.35, 0.50, 0.71, 0.97, 1.35 and 1.80 m (HCP 0.35 m to HCP 1.35 m). With this setup, depth ranges between 0-0.24 to 0-2.70 m were achieved (see Figure 2).
Electromagnetic induction transect measurements were performed 8th of March (during tillering of the plants), 27th of June (during stem elongation) and 12th of September 2019 (after harvest and power harrowing) together with the areal EMI survey. The reason to perform repeated EMI surveys at different dates over the season, associated with different soil water conditions along the transects, was to avoid a 'lucky shot' and to find out if there are environmental conditions (mainly different soil water status), where compaction can be detected better or worse.
During the measurements, the EMI instruments were mounted on two plastic sledges, which were separated by 1.5 m. For transect measurements, the sledges were pulled by hand, while for the area measurement the sledges were pulled by a lawn mowing tractor with a distance of 4 m to the first sledge and a driving speed of 5 to 8 km h −1 . By using the sledges, effects of the operator handling and the influence of terrain roughness were reduced. At the used driving speed, a sampling frequency of 5 Hz resulted in an inline measurement resolution of around 0.3 to 0.4 m on tracks with an approximate separation of 2 m between the tracks for the areal measurements. Center-point RTX DGPS systems (Trimble Inc.) with high accuracy were used to provide spatial positions of the EMI systems during all measurements.
The raw ECa measurements were filtered using the filtering strategy suggested by von Hebel et al. (2014) and calibrated following the approach of Lavoué et al. (2010) and von . Thereby, the vertical electrical conductivity distribution of inverted electrical resistivity tomography (ERT) data along a transect was inserted in an electromagnetic forward model for a horizontally layered half-space (van der Kruk et al., 2000) together with the EMI system specifications (coil separation, coil orientation, frequency) to predict ECa values. These modelled ECa values were plotted against measured ECa values and a linear regression provided shifting and scaling factors for each coil configuration (Lavoué et al., 2010;von Hebel et al., 2014). Such calibration was done to overcome physically illogic (negative) ECa values and to account for any influences of the sledges on the EMI signal. Additionally, as raw EMI measurements are often prone to systematic errors, they can only be interpreted qualitatively as already stated by Binley et al. (2015), and the calibration therefore transforms the raw ECa signal to quantitatively meaningful data (von Hebel et al., 2014Hebel et al., , 2018. For detailed information on the calibration of EMI measurements using ERT, we refer to the paper of Lavoué et al. (2010). The ERT measurements were performed together with the EMI measurements three times during the investigation ( Figure 1c). Therefore, a Syscal Pro Switch resistivity meter (IRIS instruments) was used on the two 30-m transects with 120 electrodes (0.25-m electrode spacing) measuring in dipole-dipole mode. The raw ERT readings were inverted using the RES2DINV software (Loke & Barker, 1996) with the default damping parameters. Here, it has to be noted that neither the EC (measured by ERT) nor ECa (measured by EMI) has been temperature corrected as the vertical temperature profile is assumed to be the same along the transects for both measurements and as the repeated measurements from different days were not qualitatively compared with each other.

| Statistical analysis
For the statistical analyses, MATLAB (MathWorks Inc.) was used. To test on significant lateral differences between soil sampling locations, boxplots were calculated and a Wilcoxon rank-sum test using the MATLAB ranksum function with a probability of p = .05 was performed to assess statistical dependencies. Linear regression was used to describe the relationship between soil characteristics and ECa values from sampling and the squared Pearson correlation coefficient R 2 was used as statistical measure.

| Soil properties
Soil texture analyses (n = 48) of the three predefined depth intervals along the two transects ( Figure 3) showed varying silt contents between 58% and 83%, clay contents between 10% and 26% and sand contents between 2% and 17%, indicating that the material used for recultivation was not homogeneous. In general, sand and clay contents decreased along both transects, while the silt content increased towards the end of the transects (at 30 m). Only the uppermost (0-0.3 m) layer of transect 1 showed increasing sand and clay and decreasing silt content in the first 12.5 m. Textural changes along transect 1 were most visible from a distance of 15 m. Along transect 2 differences in textures were detectable between measurements up to 7.5 m and those beyond, whereby the clay contents were highest at 2.5 m associated with lower silt contents.
Dry bulk density (hereafter named bulk density) measurements along both transects (see Figure 4)   the topsoil (0-0.3 m) had the lowest bulk densities between 0.81 and 1.35 g cm −3 . This was expected as this layer is the plough horizon, which is regularly tilled. Additionally, no clear spatial trend was found in this uppermost layer, as any natural differences in bulk density were overruled by the tillage practice. The intermediate (0.3-0.6 m) and bottom (0.6-0.9 m) depth intervals showed similarly high bulk densities ranging between 1.43 and 1.72 g cm −3 , which slightly decreased along the transects indicating a higher compacted zone at the start of the transects and lower compaction towards the end. This trend was more evident for transect 1 than for transect 2, where isolated higher bulk densities in a depth of 0.5-0.7 m even led to an increase in average bulk densities at a distance of 10-15 m along the transect (see Figure 4d). To test on significance between the horizontal locations along the transects, all bulk densities below the plough horizon were used to calculate box plots, shown in Figure 5. Different labels at the boxes indicate significant differences (p < .05) between locations, whereby same letters indicate no significant difference from each other. While for transect 1 all sampling locations up to a distance of 15 m could be grouped into one class and only the locations 22.5 and 30 m were statistically different, transect 2 showed a more scattered picture, in which only those points located at the start of the transect differed significantly from those taken at the very end.
In comparison with the static soil properties, soil texture and bulk density, the soil water content (SWC) changes over time because of precipitation, evapotranspiration and deep drainage. Gravimetric SWC was measured twice (8th of March and 27th of June) along transect 1 and only once along transect 2 (27th of June) and converted to volumetric SWC ( Figure 6). As the aim of the study was to detect areas with harmful compaction within the field and as the previous transect measurements gave confident that bulk densities can be detected, water contents were not measured on the day the areal EMI measurements were performed. Along transect 1 (Figure 6a), drier soil conditions were detected for the June measurements, which was expected, as only little rainfall was recorded from May on and as the plants extracted water during their growth. Additionally, deeper zones were generally wetter. For both sampling dates, the SWC of the plough horizon (0-0.3 m) was quite variable along the transect and there was no detectable trend. In March, the SWC below the plough horizon (0.3-0.6 and 0.6-0.9 m) slightly decreased along the transect with higher SWC at the first 15 m and lower ones after. In June, a wetter zone between 7.5 and 22.5 m distance could be found, which was in contrast to the general decrease of SWC along the transect found before. The SWC trend measured in June along transect 2 ( Figure 6b) was comparable to the trend found in March along transect 1, with a general decrease of volumetric water contents from 0 to 30 m.

| Electrical conductivity measurements
The electrical conductivity distribution along the two transects was measured three times using ERT and EMI. The inverted ERT measurements revealed a higher average electrical conductivity for transect 1 compared with transect 2. As an example, the measured bulk electrical conductivity (EC) from 12th of September is shown in Figure 7. As can be seen, along transect 1 a zone of higher EC up to a distance of 20 m was detectable, up to a depth of about 2 m. In comparison, high EC values were only found in a fairly small volume at the beginning of transect 2 (up to 5 m). Additionally, transect 2 was more variable with areas of low and higher EC over depth and length, especially after 7.5 m transect length, which was in a good agreement with the bulk density measurements, shown in Figure 5. The inverted ERT EC values were used to calibrate the apparent electrical conductivity measurements of the two EMI devices, whereby this calibration had to be performed for each individual measurement day separately. The corresponding results of the calibrated EMI data are shown in Figure 8. Transect 1 was measured three times (in March, June, and September), while transect 2 was only measured in June and September. Over all measurements, EMI derived ECa varied between a maximum value of 41 to a minimum of 7 mS m −1 , with a maximum range per measurement date of 23 mS m −1 . Generally, ECa increased with depth of investigation, and in almost all measurements (dates and configurations), a decrease of ECa along the transects was T A B L E 1 Squared Pearson correlation coefficient (R 2 ) between measured soil properties (soil texture and bulk density) and volumetric water content for the depth intervals 0-0.3, 0.3-0.6 and 0.6-0.9 m. High correlation coefficients are coloured in green, while low coefficients are red. Note that the volumetric water content was measured in March and June detectable. However, the measurement of transect 1 in June showed an increase of ECa until a distance of 15 m and a decrease after, which was in good agreement with the water content measurements discussed before. EMI measurements showed a higher variability over depth in June and September, which for June, corresponded to a higher variability of water content over depth. Note that in September no water content measurements were performed, which could be compared with measured EMI derived ECa. Therefore, a more detailed analysis of the temporal changes in ECa over time is restricted. Additionally, no ECa data were temperature corrected as only ECa data from the same day were regressed to soil properties. Nevertheless, the two dates for which water content information and ECa readings were available indicated that the ECa pattern always follows the changes in bulk density. This was also confirmed by the correlation analysis discussed below, giving confidence that the bulk density changes can be detected even under different soil water contents. From a pure visual inspection, the EC and ECa distribution along the transects already corresponded to the observed soil textural changes with higher silt contents at the end of both transects associated with lower clay contents (see Figure 3). Also, the bulk density patterns followed the same trend with higher bulk densities at the start of the transects (Figure 4).

| Correlation between soil characteristics and ECa
The squared Pearson correlation coefficient (R 2 ) was used to analyse the linear relationship between different soil parameters and volumetric soil water contents (Table 1) as well as the relationship between soil parameters/soil water contents and measured ECa for the nine different depths of investigations (Table 2). In general, the correlation between soil texture and bulk density as well as volumetric SWC of the top layer (0-0.3 m) showed small correlations (R 2 < .4), which could be explained by regular tillage of the plough horizon. Nevertheless, a relatively high correlation was found between bulk density and volumetric SWC for both measurement dates of soil water content. With R 2 of .85 for March and R 2 of .76 for June for 0-0.3 m, the correlations were even slightly higher as the correlation between bulk density and water content at greater depths, where no direct impact of tillage occurred.
Additionally, correlation between bulk density and SWC was higher for the SWC measurements in March (0.3-0.6 m R 2 = .79; 0.6-0.9 m R 2 = .75) compared with the SWC measurements in June (0.3-0.6 m R 2 = .36; 0.6-0.9 m R 2 = .54). The change in the correlation strength seemed to be attributed to the degree of saturation, whereby the SWCs were in general higher in March compared with June (see Figure 6), because of low rainfall and large evapotranspiration between the two dates.
Correlation between bulk density and soil texture was low to medium for the lower sampling depth with R 2 of .42 (0.3-0.6 m) to .64 (0.6-0.9 m) for silt and R 2 of .48 (0.3-0.6 m) to .62 (0.6-0.9 m) for sand. However, the correlation with clay was weak (R 2 < .34) and no correlation between soil texture and bulk density was found for the topsoil (0-0.3 m) likely because of the tillage practice. Unfortunately, different authors reported either good or no correlation between soil texture and bulk density. On the one hand, Jones (1983) showed no correlation between both variables, irrespectively of a large variability of clay and bulk density. On the other hand, Chaudhari et al. (2013) found high correlation for bulk density with clay (r = −.63), sand (r = .90) and silt (r = −.73). Bernoux et al. (1998) for example, reported only low correlation between bulk density and soil texture (especially clay) and stated that additional soil properties such as organic matter or pH do also have an impact on bulk densities.
In a next step, the soil properties and volumetric SWC were correlated against ECa values obtained from the nine different EMI configurations. As EMI measurements were conducted three times (8th of March, 27th of June and 12th of September) on the transects, also the temporal changes in the correlation between bulk density, as the target of interest, and measured ECa can be analysed. As could be clearly seen, the correlation coefficients for the plough layer (0-0.3 m) were low for all measurement dates (R 2 < .33) and in most cases even below .2. As a result of tillage, the plough horizon was rather heterogeneous, and therefore, the point bulk density measurements could not be captured by the EMI measurements, which were integrated along the distance of the coils of the instruments. Additionally, only the VCP 0.32 m configuration was sensitive to the plough horizon only, whereas all other configurations had deeper DOIs. For the deeper layers, relatively high R 2 values between .58 to .92 for the 0.3-0.6 m layer and .75 to .97 for the 0.6-0.9 m layer were found for the March measurements. For the June measurements, the correlation was slightly worse and no correlation was found for the VCP 0.32 m configuration, which was expected, as this configuration measured quite shallow (DOI 0-0.2 m). In comparison with the June measurements, the correlations slightly increased again for September, but were in general lower as those found in March, especially for the shallow sensing configurations (VCP 0.32, VCP 0.71, HCP 0.35, HCP 0.49 and VCP 1.18 m). Overall, the results revealed that the correlation between bulk density and ECa varied between the different measurement dates but revealed acceptable results for all measurement dates. Differences in the predictive power (R 2 ) can be T A B L E 2 Squared Pearson correlation coefficients (R 2 ) between EMI measured ECa (VCP 0.32 m to HCP 1.8 m) and soil properties (soil texture and bulk density) as well as volumetric water contents for the depth intervals 0-0.3, 0.3-0.6 and 0.6-0.9 m measured at three different measurement dates. High correlation coefficients are coloured in green, while low coefficients are red. Correlations that were used to create the probability map of soil compaction zones are written bold explained by changes in SWC between the different dates. These changes influenced the ECa readings as well as the pore water conductivity, which will also change if the pore saturation will change (Corwin, 2008;Doolittle & Brevik, 2014). As the June measurements were performed under the driest conditions, it can be recommended to perform EMI measurements for bulk density detection not in the dry soil state. This is in contrast to the study of Al-Gaadi (2012), who found a decreasing correlation with increasing water content in sandy soils. The contradictory findings between the study of Al-Gaadi (2012) and our results could be explained by the differences in soil texture (sandy soil compared with loamy soil) and the water holding capacity, especially at drier conditions. Whereby it is generally known that sandy soils exhibit only small water contents at dry conditions, whereas loamy soils still contain relatively large amounts of water.
As the areal EMI measurements were performed on 12th of September after harvest, along with the transect measurements, these data will later be used on for the areal prediction of the bulk densities. Therefore, these data listed in Table 2 will be discussed in more detail. As already discussed, correlation was weak for the plough horizon. For the deeper layers, the predictive power increased to moderate or even high. Thereby, the VCP 0.32 and VCP 0.71 m showed lowest correlation with R 2 < .53, which was expected, as these configurations were mainly sensitive to the shallow soil (DOI of 0-0.2 and 0-0.5 m). Nevertheless, these R 2 values were in the range found also by Sudduth et al. (2010) for the correlation between EMI ECa and penetrometer cone index. For all other configurations, the R 2 exceeded .65 and reached up to .84 (HCP 1.80 m, DOI 0-2.7 m). The relatively low R 2 for the HCP 0.49 m (DOI 0-0.7 m) was somehow surprising, as the next less shallow configuration (HCP 0.35 m, DOI 0-0.5 m) as well as the next deeper one (VCP 1.18 m, DOI 0-0.9 m) indicated higher predictive power. The reasons for this behaviour are unknown.

| Areal EMI measurements
The areal calibrated EMI measurement conducted on the 12th of September showed ECa values varying between a minimum of 11.3 and a maximum of 51.7 mS m −1 and increasing ECa with depth ( Figure 9). As expected and already discussed for the EMI transect measurements (Figure 8), not only lateral, but also horizontal variability was detected. Looking at the entire area, several zones exhibiting higher ECa values over all depths were detectable, especially in the centred and south-eastern part. By visual inspection, these zones almost matched the zones with observed low plant productivity during the previous growing season.

| Bulk density prediction
As a result of the high correlation between bulk density and EMI measured ECa for the layers below the plough horizon along the transects, bulk density and ECa were regressed. The plough layer was not considered in the regression because of low correlation and the fact that the study focused on deep soil compaction. For the regression, the ECa data measured along the transects on 12th of September were used as on this date also the areal EMI measurements were available. For an exemplary prediction, the regression between bulk density and HCP 0.97 m was used as this EMI configuration yielded high correlation for both layers to be predicted (0.3-0.6 m R 2 = .83, 0.6-0.9 m R 2 = .83). Here, it has to be noted that also other EMI configurations such as HCP 1.35 or HCP 1.80 m could have been used as they also showed high R 2 values for both depths. The regressions between bulk density and HCP 0.97 m for both layers are shown in Figure 10 a and c. Based on these regressions and the calibrated areal ECa data shown in Figure 9, the bulk densities for the two layers (0.3-0.6 and 0.6-0.9 m) were calculated and plotted in Figure 10 b and d. The predicted areal bulk densities varied between 1.5 and 1.8 g cm −3 and, as expected, they showed the same pattern as the ECa readings plotted in Figure 9 across the field and a slightly increase of bulk densities with increasing depth.
As a result of measurement (ECa and bulk density) and prediction errors in the regression, quantitative predictions F I G U R E 1 0 Exemplary linear regression of EMI measured ECa [mS m −1 ] using HCP-oriented coils with a separation of 0.97 m vs. bulk density [g cm −3 ] for the depth intervals (a) 0.3-0.6 and (c) 0.6-0.9 m. Bulk density predictions for the depth intervals (b) 0.3-0.6 and (d) 0.6-0.9 m using the linear regression equations of (a) and (c), respectively. High bulk density predictions are coloured in red, while low ones are blue and the locations of transects are given in black of the bulk density using this method also contains uncertainties. Therefore, instead of providing quantitative predictions, we calculated probability maps to delineate zones of harmful soil compaction below 0.3 m. To generate these probability maps, all areal bulk density predictions of the individual EMI configurations were used, which exceeded a R 2 value of .75 (see Table 2). Based on the regressions for each configuration, individual maps were generated for an equally spaced grid. In a next step, a compaction threshold was defined and applied to the individual maps. Finally, all information from the individual thresholded maps for each grid point were overlaid and the probability was calculated, as shown in Figure 11. In this study, the threshold value was chosen to be 1.7 g cm −3 , as it is a good compromise of ranges of bulk densities for harmful soil compaction in silt loams reported in literature. Literature values of minimum bulk density from which harmful, root growth effecting or even restricting soil compaction starts in a silt loam vary between 1.60 and 1.75 g cm −3 (Arshad et al., 1997;Eckelmann, 2005;Jones, 1983;Usaborisut & Ampanmanee, 2015). In general, harmful soil compaction depends on several soil properties like soil texture and structure but also on organic carbon content (e.g. Reichert et al., 2009;Usaborisut & Ampanmanee, 2015). For example, critical bulk densities for root growth have been found to decrease with higher clay and clay-plus-silt contents (e.g. Jones, 1983). Therefore, a global single bulk density value as threshold for harmful compaction cannot be determined. Instead, a soil-specific value needs to be defined, which accounts for the given soil textural and sorting information. In general, the defined threshold should be higher as bulk densities classically found for the soil investigated.
The probability map revealed that the bulk density values of a large proportion of the field were below 1.7 g cm −3 (green zones), which should not affect root growth and crop performance. Nevertheless, a small proportion (yellow zones) indicated that the bulk density likely exceeded the threshold of 1.7 g cm −3 , as the probability of all regressions from Table 2 was larger than 50%. Finally, the red coloured zones indicated those areas where a harmful subsoil compaction for the two different depths could be estimated with high certainty, as the probability of all regressions was 100%. Looking at the areal statistics, 91% and 83% of the area was not affected by harmful subsoil compaction (green zones) for the layers 0.3-0.6 and 0.6-0.9 m, respectively. Harmful compaction (red zones) was predicted for 7% and 6% for the two depths and only a small proportion was in the transition zone with 50%-100% probability (yellow zones) with 2% and 11% for the two depths.

| SUMMARY AND CONCLUSION
The present study investigated soil compaction by means of electromagnetic induction (EMI) and electrical resistivity tomography (ERT) validated with ground truth soil samples along transects to assess large-scale usability of EMI to identify compacted zones in a recultivated field. The geophysical methods measured the spatially distributed apparent electrical conductivity (ECa), while the soil samples were analysed regarding bulk density, volumetric soil water content and soil texture. Linear regression was used to describe the statistical relationship between soil characteristics and ECa and the squared Pearson correlation coefficient was determined as statistical measure to express the quality of the relationships. The EMI transect measurements showed good correlations to the bulk density measurements taken along the same transects. The correlations improved with soil depth and increasing volumetric water content. The retrieved linear F I G U R E 1 1 Probabilistic map of high compaction risk zones for the depth intervals (a) 0.3-0.6 and (b) 0.6-0.9 m. Zones that reveal high probability to be harmfully compacted are coloured in red, medium probability yellow and zones that have a high probability to be uncompacted are coloured in green 180 BD < 1.7 g cm -3 BD ≥ 1.7 g cm -3 BD ≥ 1.7 g cm -3 BD < 1.7 g cm -3 BD ≥ 1.7 g cm -3 BD ≥ 1.7 g cm -3 UTM-N[m] + 5660442 (b) Depth 0.6-0.9 m (a) Depth 0. relationship between ECa and bulk density was further used to predict a spatial distribution of bulk density based on the areal ECa values measured with EMI. By this, we were able to upscale the local measurements of bulk density and to create a probability map of zones that are likely affected by harmful soil compaction. For applications within the recultivation process, we recommend the following protocol based on our findings. (1) Areal multi-coil EMI measurements under bare field conditions should be performed, best one before first sowing and one after harvest to cover two different soil moisture conditions as compacted zones seem to be best detected under dry or wet soil status, depending on the soil type.
(2) Based on the areal measurements, consistent locations/patches with strong ECa contrast (areas with low and areas with high ECa) should be identified and used for ground truthing. At these locations, sparse bulk density measurements (or as an alternative soil resistance measurements) should be performed, which will allow to validate bulk density anomalies. If required, also pits can be opened for ground truthing. As also changes in soil texture will impact the ECa readings, soil texture should be determined either in the laboratory or by finger probing to exclude textural changes from compacted zones.
(3) Based on the soil texture information and sorting, a soil-specific bulk density threshold should be defined from which root growth effecting or even restricting soil compaction starts. (4) Based on the soil bulk density measurements and known ECa, regressions should be performed and compaction probability maps should be established, using the information of all EMI coils having good correlations to bulk density, to identify harmfully compacted zones. This cost-efficient and fast method would allow a quick intervention for planning soil melioration measures at early recultivation stages, which are of immediate importance for, for example, mining companies that have to make recultivated land available to farmers for long-term and fertile cultivation.