Seismic signatures reveal persistence of soil compaction

A well‐developed soil structure is a key attribute of a productive and functioning soil. Evidence shows that subtle changes in the spatial arrangement and binding of soil constituents impart large changes in soil mechanical and hydraulic properties and associated ecological services. However, these features remain difficult to quantify at spatial scales relevant for agricultural management. In this work, we propose a pedophysical model to interpret macroscopic seismic properties in terms of soil structure. The model captures subtle soil mechanical traits accounting for soil plastic deformation due to compaction. In order to evaluate the model, we use data from field monitoring at an experimental site that revealed elevated seismic velocities in plots that were compacted 5 yr prior to our measurements. Our results show that P‐wave velocities carry a strong imprint of soil compaction and are well predicted by the proposed model. The model infers contact areas between aggregates that are nearly threefold larger for compacted than for non‐compacted soils, indicating that soils have not recovered from compaction. The study illustrates the potential of seismic methods to identify chronic compaction at field scale.


INTRODUCTION
Soil structure refers to the spatial organization of a soil's solid constituents (minerals and organic matter) and voids resulting from bioturbation, abiotic processes (e.g., wetting-drying and freezing-thawing), and soil management practices (Gregory et al., 2009;Oades, 1993;Or et al., 2021). Soil structure is critical for determining a soil's macroscopic transport and mechanical properties. Well-developed and stable soil structure plays a central role in supporting soil ecosystem services, including climate regulation (Fatichi et al., 2020;Hirmas et al., 2018), groundwater recharge (Keesstra et al., 2012), (Nawaz et al., 2013). The quantification of compaction remains a challenge, particularly at the field scale. Moreover, favorable soil structure takes decades to develop primarily by biological processes, whereas compaction may occur in a few seconds, with its adverse impacts on soil functions lasting for many years thereafter (Håkansson & Reeder, 1994;Webb, 2002). Thus, gaps in our understanding of long-term dynamics of soil structure can be associated with (a) the disparities in the timescale of degradation and generation processes, and (b) the lack of characterization methods capturing the relevant temporal and spatial scales . For these reasons, the global extension of soil compaction in all terrestrial surfaces remains unknown and challenging to assess. Three decades ago, Oldeman (1991) estimated that 68 Mha of arable lands were compacted globally. Recent estimations indicate that about 25-40% of all arable land is compacted in the United Kingdom (Graves et al., 2015), Denmark (Schjønning et al., 2015), and the Netherlands (Brus & van den Akker, 2018). Characterization of the structure and mechanical status of a soil often relies on pointwise on-site measurements (e.g., cone penetrometer testing; Lunne et al., 1997), or on highly invasive soil evaluation methods (Guimarães et al., 2017), often involving time-intensive characterization based on soil sampling (Peng & Horn, 2008). Such invasive approaches provide fragmentary information in time and space and offer a limited capacity to evaluate soil structure evolution under natural conditions. Recently, geophysical methods have been proposed to add complementary information to traditional techniques and to bridge the gap from point to field scales at relevant temporal scales (see review by Romero-Ruiz et al., 2018). Because of the direct link between seismic velocities and soil elastic moduli, seismic methods hold a particular promise for minimally invasive characterization of variations in soil mechanical status at field scales and at highly resolved timescales. On-site applications of shallow seismic methods and monitoring have increased in the last decades including monitoring of thawing in arctic environments (Stemland et al., 2020;Wu et al., 2017), landslide monitoring (Chen et al., 2018;Whiteley et al., 2020), and characterization of soil compaction (Donohue et al., 2013;Keller et al., 2013;Uyanik, 2011;Uyanik & Ulugergerli, 2008). These studies highlight that seismic methods can provide valuable information pertaining to soil structure. However, existing studies that analyze soil structure based on seismic velocities are largely empirical, soil specific, and lack descriptions of arrangement of soil components and the mechanisms behind the soil structure state (for linking deformation, macroporosity, and mechanical state with seismic signals). The development of seismic-based practices for characterization of soil structure at relevant scales will largely rely on continued improvements in our understanding of the relationships between the macroscopic seismic properties of structured

Core Ideas
• We monitored seismic velocities in an experimental field 5 yr after a compaction event. • We measured higher P-wave velocities in compacted than in non-compacted treatments. • A novel pedophysical model was used to link seismic velocities to inter-aggregate contact areas. • We estimated a nearly threefold increase in contact area in compacted relative to non-compacted treatments. • Seismic methods offer great potential for soil structure characterization at management scales.
soils and the spatial arrangement and binding of soil constituents.
In this study, we seek to improve mechanistic understanding regarding the potential of seismic methods in quantifying the mechanical impacts of soil compaction. The objectives of this study are (a) to develop a soil-structure-informed pedophysical model that accounts for soil structure deformation; (b) to monitor seismic signatures of different states of soil compaction at a long-term soil structure observatory (SSO) established in Zürich, Switzerland ; and (c) to apply the model to interpret seismic signatures of persistent soil compaction 5 yr after a compaction event at the SSO.

ELASTIC PROPERTIES OF STRUCTURED PARTIALLY SATURATED SOIL
Seismic induced compressional wave (P-wave, V p ) and shear wave (S-wave, V s ) velocities can be expressed as functions of soil's elastic moduli and bulk density ρ b . For isotropic elastic media, the seismic velocities are expressed as where K soil and G soil are the macroscopic soil bulk and shear moduli, respectively. Relating these properties to their microscopic counterparts remains challenging for structured soils as they depend on both (a) the elastic properties and Step 2 (modified from Ghezzehei & Or, 2003) volumetric fractions of the soil phases (solids, water, and air), and (b) the spatial distribution of the soil constituents and their corresponding mechanical bonds. In this work, the incorporation of salient features of soil structure for predicting seismic response of unsaturated soils is based on self-consistent mixing of soil constituents which, in turn, is based on a conceptual model of their arrangement. We conceptualize the soil as a composite dual-domain medium (Figure 1a). Such a dual-domain representation of rocks and soils has been successfully used to model seismic (Dvorkin et al., 1999), hydraulic (Durner, 1994;Tuller & Or, 2002), dielectric (Blonquist et al., 2006), and electrical properties (Day-Lewis et al., 2017). The dual-domain pedophysical model proposed herein accounts for soil structural features by considering four mixing steps (Figure 1b). These steps are described in the sections below.

Elastic properties of soil aggregates and soil frame
The elastic properties of the soil frame are obtained by applying the Hertz-Mindlin (HM) model in two consecutive steps. First, we obtain the elastic properties of soil aggregates (Step 1), and then we consider a geometrical arrangement of aggregates (Step 2). The dual-porosity soil with total porosity of ϕ T comprises an intra-aggregate domain and an inter-aggregate domain occupying volumetric fractions w m and w M = 1− w m , respectively. The aggregates are composed by pores and solid particles and have an intra-aggregate porosity of ϕ m . The inter-aggregate domain (space between aggregates) is composed by void space only (thus, making ϕ M = 1). The total porosity can thus be expressed as The HM mechanical contact model is used to derive the elastic moduli of an aggregate (Step 1 in Figure 1b) that depend on the normal and tangential contact stiffness between adjacent particles (pp. 245-248, Mavko et al., 2009). The elastic moduli are expressed in terms of the elastic properties of the particles and a confining pressure P e as agg = where K agg and G agg are the bulk and shear moduli of the aggregate, respectively. Furthermore, N m is the average number of contacts per particle, which can be considered a function of the intra-aggregate porosity (pp. 232-234, Mavko et al., 2009). In addition, G s and ν s are the shear modulus and Poisson's ratio of the particles, and f m is the fraction of nonslipping particles (see Bachrach et al., 2000). The elastic properties of the soil frame ( Step 2) are obtained by considering an assembly of spherical aggregates, whose homogeneous isotropic and elastic properties are given by Equations 4 and 5. For this purpose, we use the HM model once again. The elastic properties of the soil frame are then expressed as where K frame and G frame are the bulk and shear moduli of the soil frame, ν agg is the Poisson's ratio of the aggregates, N M is the average number of contacts per aggregate, f M is the fraction of non-slipping aggregates, and P a is the HM confining pressure for soil aggregates. The confining pressure P e is used as in traditional geophysical applications of the HM framework, attributed to overburden pressure and suction stress for partially saturated media (see Section 2.1.1).
In this study, we expand the interpretation by introducing P a to relate HM elastically deformed soil aggregate contacts and permanently (viscous) deformed soil aggregate contacts due to compaction induced by a transient vertical stress (see Section 2.1.2). Hence, these signatures of soil compaction are interpreted via their elastic equivalence in the HM framework that, in turn, affects the seismic response.

Effects of hydration status on the elastic properties of the soil frame
The elastic moduli of the soil frame (Equations 6 and 7) are strongly dependent on the elastic properties of the aggregates, which depend on the confining pressure P e (Equations 4 and 5). Research shows that P e and, therefore, the elastic moduli of the frame are affected by the hydration status in soils (Brutsaert & Luthin, 1964;Lu, 2011;Shin et al., 2016). A frame softening effect is related to the decrease in the aggregate matric suction through the concept of effective stress (Bishop & Blight, 1963;Lu et al., 2010;Nuth & Laloui, 2008).
As proposed by Shen et al. (2016), we account for soil frame softening by adding a suction stress term [χ( a − w )] to the HM pressure appearing in Equations 4 and 5: where σ = ρ b 0 is the overburden pressure, a = ρ a 0 is the air pore pressure, u w is the water pore pressure, ψ agg = a − w is the matric suction of the aggregates, g is the gravitational acceleration, and Z 0 is the depth of investigation. The effective stress parameter χ is a function of the effective saturation of the aggregates e agg and varies between 0 and 1. In its simplest form, it is expressed as χ = e agg (Nuth & Laloui, 2008), yet the dependence of χ with saturation may vary with the soil texture (Khalili & Khabbaz, 1998). The bulk density is defined as where ρ s , ρ w , and ρ a are the densities of soil particles, water, and air, respectively, and (= θ∕ϕ T ) is the total soil water saturation defined as the ratio between total water content θ and total porosity. Consequently, the effective moduli of the soil aggregates become dependent on water content through the water retention curve of the aggregates. The matric suction of the aggregates ψ agg can be expressed as a function of effective water saturation of the soil aggregates e agg as (van Genuchten, 1980) where α m is related to the inverse of the air-entry pressure of the aggregates, and m is the Van Genuchten exponent, which is related to the pore size distributions of the intraaggregate space (Mualem, 1976;van Genuchten, 1980). The effective saturation in the aggregates is a function of the residual water content θ r , the intra-aggregate porosity ϕ m , and the intra-aggregate water content { e agg = (θ agg − θ r )∕[ϕ m (1 − M ) − θ r ]}. In our model, the intra-aggregate water content is obtained from the soil total water content (θ; e.g., from timedomain reflectometry [TDR] data) and the volumetric proportion of intra-aggregate pore space. If the total soil water content is lower than or equal to the volumetric proportion of intra-aggregate pore space [θ < ϕ m (1 − w M )], the water content in the aggregates is the total water content in the soil (θ agg = θ). For a total soil water content higher than the proportion of intra-aggregate pore space [θ > ϕ m (1 − w M )], the water content in the aggregates is equal to the volumetric proportion of intra-aggregate pore [θ agg = ϕ m (1 − M )].

2.1.2
Soil compaction effects on the elastic properties of the soil frame In the framework of the HM theory expressed by Equations 6 and 7, the soil frame is represented as an assembly of soil aggregates with contacts elastically deforming in the presence of a confining pressure P a . In reality, the aggregate contacts in soils deform plastically due to soil stresses induced by abiotic processes (e.g., wetting-drying and freezing-thawing cycles), anthropogenic activities (e.g., passages of agricultural vehicles), and/or biological activity (e.g., penetrating earthworms and growing roots) producing aggregate coalescence. Below, we describe how to account for permanent deformation of the soil frame using the HM framework through Equations 6 and 7 via the relationship between P a and the corresponding contact areas between aggregates.
The application of a transient load (e.g., passage of a tractor) results in an elastic (temporary) and a viscous (permanent) deformation of the soil frame producing an axial strain ε = ε e + ε v (see Ghezzehei & Or, 2003), where ε e and ε v are the elastic and viscous strains, respectively. The lasting effect of a soil compaction event is associated with irreversible deformation comprised in ε v , which can be modeled using information on the initial (prior to compaction) strain ε 0 , the axial load and duration of stress application (i.e., tractor weight and speed), and the soil rheological properties represented by the Bingham model (Ghezzehei & Or, 2001). For small deformations (ε v ≪ 1), the post-compaction viscous component of the strain (reflecting irreversible soil deformation) ε v is related to the contact area between aggregates c through the approximation where a and R are the radii of the contact area and the aggregate, respectively. We do not seek to reconstruct the deformation process (Ghezzehei & Or, 2001;Or & Ghezzehei, 2002) but to infer from the primarily elastic nature of seismic methods what is the representative contact area between aggregates based on the HM contact theory. To build such a pedophysical model, we rely on the contact area (Equation 11) that is permanent in nature, yet for convenience, it is treated here as an elastic deformation in the context of the HM theory. This allows expressing P a in terms of the contact radius a by (p. 246, Mavko et al., 2009) To assess compaction independently of the aggregate size, we define the relative contact radius rel = ∕ as the ratio between the radius of the contact area and the radius of the aggregates, which can be related to the viscous strain using Equation 11: ε ≈ 1 2

Elastic properties of the water-air fluid mixture
The effective bulk modulus of the fluid mixture (Step 3 in Figure 1b) is required for evaluating the elastic moduli of a partially water-saturated structured soil. In the low frequency limit, wave-induced pore fluid pressure gradients induced by the propagating seismic wave-field have sufficient time to equilibrate in the pore space (Müller et al., 2010). This allows representing the elastic properties of the fluid mixture as properties of an effective fluid (composed by water and air) using the Reuss isostress average, which is based on an harmonic mean (Reuss, 1929) and relates the bulk modulus of the fluid mixture K fluid to the saturation ( = θ∕ϕ T ) by where K w and K a are the bulk moduli of water and air, respectively. Due to the low frequencies used in this study (i.e., wavelengths are much larger than the inter-and intraaggregate pore radii), the pore fluid pressures can be considered as constant, and thus the effective fluid is considered to be occupying the full pore network (comprising both interand intra-aggregate pores).

Elastic properties of structured partially saturated soil
Gassmann's low-frequency fluid substitution relationships (Gassmann, 1951) are widely used to systematically incorporate the effects of a saturating fluid in the elastic moduli of porous materials (Bachrach & Nur, 1998;Pasquet et al., 2016). Gassmann's relationships can be used to predict the elastic moduli of soils by considering the elastic moduli of (a) the soil frame (Equations 6 and 7), (b) the soil solid constituents forming the soil frame K s , and (c) the saturating fluid (Equation 14). Gassmann's relationships using the total porosity are valid in the low frequency limit even for complex pore networks (e.g., cracked porous media in Gurevich et al., 2009). Effects of partial water saturation are modeled by considering an effective saturating fluid, whose properties are given by Equation 14 (e.g., Johnson, 2001). In this framework, the elastic moduli of the partially saturated structured soil are expressed as (p. 273, Mavko et al., 2009) and Equations 15 and 16 are valid in the so-called relaxed state. This implies that pore fluid pressures are equilibrated throughout the pore space (Berryman, 1999). Dissipative effects associated with the relative motion between the saturating fluids and the frame, which are commonly considered in the context of Biot's poroelasticity theory (Biot, 1962), are thus assumed to be negligible in the frequency range used in our study. Note that the macroscopic elastic properties (and seismic velocities) of the soil are sensitive to changes in water content (a) through the water retention relationship of the aggregates, which determines their effective stress without dependence on the elastic moduli of water (Equation 8), and (b) through the impact of water saturation in the elastic properties of the fluid mixture (Equation 14).

FIELD MONITORING OF SEISMIC DATA
We conducted seismic monitoring to evaluate long-term effects of soil compaction and analyze potential recovery for different soil cover. Monitoring was carried out in the spring and summer of 2019 at an experimental field site located in the vicinity of Zürich, Switzerland (47˚25′39 N, 8˚31′04 E; Keller et al., 2017). This SSO is a long-term experiment designed to study the evolution of soil structure, after a compaction event in the spring of 2014, for different types of post-compaction management. The SSO has a strip-plot design with three blocks (replicates, see Figure 2a). We monitored the seismic response of experimental plots (located in Block A) with different covers (bare soil and vegetated soil) and compaction treatments (compaction on the full sur-face and no compaction). The corresponding soil treatments are referred to as full compacted ley (grass-legume mixture), non-compacted ley, and full compacted bare soil (Figure 2b). Since the soil properties (and texture) prior to the compaction event were similar at all monitoring sites , we attribute differences in the corresponding seismic signatures to the different soil covers and treatments.
The seismic array comprised a line of 18 geophones (30 Hz SM-11 IO) deployed on both sides of an impact source (12 on one side and 6 on the other side, Figure 2c). Two geophone spacings were used: Δ = 10 cm for sensors with offsets shorter than 80 cm, and Δ = 20 cm for sensors with offsets longer than 80 cm. The sensors were connected to a Geode Exploration Seismograph located in an operation box at the edge of the soil plots. The Geode was controlled by a laptop operating the Seismodule Controller Software continuously during the full monitoring campaigns. The geophone array was not symmetric due to the limited length of the seismic cable used in our experiment. As shown in Figure 2b, the operation box controlling the seismic monitoring experiment was placed at an untreated soil area located a few meters away from the seismic lines. In order to connect the geophones to the Geode (located in the operation box), part of the seismic cable was used as an extension, which prohibited us from using six outtakes of the cable.
An electromagnetic hammer (piston and base plate system driven by a 300-W audio speaker) was harnessed to provide an impact-like source that was activated at regular time intervals. The source frequency content was centered between 75 and 150 Hz. The hammer was controlled by a waveform generator (Agilent 33210A, located in the operation box) that was programmed to trigger an impact every 15 min. The function generator emitted a pulse waveform, with a period of 15 min, high level of 4 V, low level of −1 V, width of 300 ms, and edge time of 100 ns. The seismic source is considered a point source, as shown in the schematic representation of the seismic array presented in Figure 2c. Having only one source, the data were collected during three different time periods: from 17 May to the 26 June (compacted ley), from 26 June to the 8 July (non-compacted ley), and from 25 July to 9 August (compacted bare soil). All these periods contain significant rainfall events followed and preceded by dry periods with associated changes in soil hydration status. Figure 2c presents a photograph from May 2019 of the seismic monitoring layout deployed in the compacted ley.
Time-domain reflectometry (TDR 100 by Campbell Scientific with MDX multiplexers) probes for soil water content measurements were installed in all plots and were continuously collecting data at four different depths (10, 20, 40, and 70 cm). Due to a technical issue, TDR data were not available from Block A in 2019. For this reason, we considered TDR data collected at the TDR banks in Block C  Figure 2a). The TDR data from previous years in Block C are in good agreement with data from Block A (see Supplemental Section S1). The soil properties are assumed to vary only as a function of depth within each of the soil treatments.

EFFECTS OF COMPACTION AND WATER CONTENT ON MEASURED SEISMIC SIGNATURES
The seismic velocities derived from the seismic monitoring data are contrasted with predictions made by feeding the pedophysical model with water content θ data derived from TDR sensors at 10-cm depth. The P-wave velocities were obtained by first-break picking (Supplemental Section S2). These seismic velocities are representative of the first soil layer with an estimated thickness of 20 cm obtained from first arrivals by solving for the depth to the first refractor. Below, we present the mean values and standard deviations of the estimated seismic velocities at different offsets.

P-wave slowness variations with soil water content
We compared time series of estimated P-wave slowness together with TDR-derived water content using Topp's equation (Topp et al., 1980) for compacted ley (Figure 3a), non-compacted ley (Figure 3b), and compacted bare soil (Figure 3c). The P-wave slowness and the water content showed a reasonable correlation regardless of the studied soil treatment with corresponding correlation coefficients of r = .7, r = .95, and r = .85 for compacted ley, non-compacted ley, and compacted bare soil, respectively (see Figure 4).

Pedophysical parametrization at the SSO
Our pedophysical model was developed to deduce information related to soil structure from the measured seismic signatures. In applying the model, we first identified the properties that are linked to compaction states and those related to soil F I G U R E 3 Time series of estimated P-wave slowness and water content for (a) full compacted ley, (b) non-compacted ley, and (c) full compacted bare soil. Note that slowness and water content correspond to different time periods for the different soil treatments F I G U R E 4 Scatter plot of P-wave slowness as a function of water content for all soil treatments texture. The different soil plots at the SSO present negligible differences in soil texture and the same pre-compaction history . Thus, we assumed that differences in soil structure for the studied soils are manifested primarily via differences in the compaction-induced contact area A c (Equation 12) and the inter-aggregate space w M (Equations 6, 7, and 12). We considered that the aggregate properties (Equations 4 and 5) are unaffected by compaction and the same for all soil treatments (see Berli et al., 2008;Or & Ghezzehei, 2002). Below, we present the values of the model parameters used (see also Table 1).
The properties of water (K w = 2.2 GPa and ρ w = 1 g cm −3 ) and air (K a = 0.101 MPa and ρ a = 1.29 mg cm −3 ) were taken from the literature (pp. 176 and 468, Mavko et al., 2009). The elastic properties of the grain mixture K s = 36.5 and G s = 27.6 GPa were chosen to be representative of the soil texture with roughly 25% sand, 25% clay, and 50% silt . As suggested by Gurevich and Carcione (2000), we modeled the grains as a mixture of quartz (sand) and clays using the Hashin-Shtrikman lower bound (Hashin & Shtrikman, 1963), implying that the softer material (clays) acts as the primary load-bearing material. The silt fraction was split in equal parts between clay and sand leading to a mixture containing 50% clays and 50% quartz. The dominant clay minerals at the SSO are illite and smectite; thus we used the elastic properties of an illite-smectite mixture (K clay = 36 and G clay = 18 GPa) reported by Wang et al. (2001). The elastic moduli of quartz (K quartz = 37 and G quartz = 44 GPa) were taken from the literature (Table A.4.1 in Mavko et al., 2009). The Poisson ratio of grains and aggregates (ν s and ν agg ) were computed from their corresponding elastic moduli assuming homogeneous isotropic elastic properties. The highest values of water content measured in 2019 were 0.44, 0.46, and 0.44 for compacted ley, non-compacted ley, and compacted bare soil, respectively. To account for air entrapment (i.e., the total porosity is higher than the maximum in situ water content derived from TDR measurements; Faybishenko, 1995; T A B L E 1 Values of the parameters used for pedophysical model predictions for compacted ley (CL), non-compacted ley (NL), and compacted bare soil (CB). Only one equation is reported when the model parameter depends on water content measurements. Column "Eq." indicates the equation where the parameter is used. The comments specify if the parameters were taken from literature values, calculated from a given model or assumed 16 a The table includes the bulk modulus K clay and shear modulus of clay G clay ; the bulk modulus K quartz and shear modulus of quartz G quartz ; the bulk modulus K s and shear modulus of the soil particles G s ; the intra-aggregate porosity ϕ m ; the number of contact per soil particle N m ; the fraction of non-slipping particles f m ; the bulk densities of soil particles ρ s , water ρ w , and air ρ a ; the total porosity ϕ T ; the water saturation S; the depth of investigation Z 0 ; the effective stress parameter χ; the van Genuchten properties α m and n m ; the effective saturation of soil aggregates e agg ; the saturated θ s and residual water content θ r ; the inter-aggregate pore space w M ; the number of contacts per aggregate N M ; the number of non-slipping aggregates f M ; the relative radius R rel ; the aggregate radius R; the viscous strain ε v ; the HM pressure for soil aggregates P a ; the contact area A; the bulk modulus of water K w , air K a , and fluid mixture K fluid ; and the bulk K frame and shear modulus of the soil frame G frame . Sakaguchi et al., 2005), we assumed the total porosities to be 0.47, 0.49, and 0.47 for compacted ley, non-compacted ley, and compacted bare soil, respectively. The velocities reported herein are to be considered as apparent velocities (i.e., averaged velocities over the first soil layer). For this reason, we selected a depth of investigation that is representative of the first soil layer. Consequently, the depth of investigation Z 0 was set to 10 cm, which is half the depth of the first soil layer and is in agreement with the depth of TDR measurements. For simplicity, we took the effective stress parameter χ to be equal to e agg . As suggested by Bachrach et al. (2000), the fraction of non-slipping particles and aggregates were taken as m = M = 0.5 for all soil treatments to account for the fact that some contacts have zero tangential stiffness. The water retention properties of the aggregates are based on the measurements reported by Carsel and Parrish (1988) on clayloam soils sharing a similar texture as the soil studied here ( m = 1.25 and α m = 0.02 cm −1 ). The intra-aggregate porosity ϕ m = 0.46 was assumed based on TDR measurements. The fraction of inter-aggregate porosity w M was calculated from the intra-aggregate porosity and the total porosity (Equation 3). The average number of contacts per particle N m = 6.5 was chosen according to data by Manegold and von Engelhardt (1933).
The average number of contacts per aggregates were calculated using the equation for dense packings proposed by García and Medina (2006) as N M = 4.46 + 9.7(0.384 − w M ) 0.48 , resulting in N M = 10.44 for the compacted treatments and N M = 10.14 for the non-compacted treatment. We used the relative contact radius (R rel , see Equation 12) as a fitting parameter. For each soil treatment, a grid search was performed to minimize the L 2 norm of the misfit between modeled and observed P-wave velocities. Minimum values were obtained when the relative contact radii were 0.162, 0.095, and 0.166 for compacted ley, non-compacted ley, and compacted bare soil, respectively, which correspond to viscous strains of ε v = 1.3%, ε v = ε 0 = 0.4%, and ε v = 1.4%. Assuming a cubic packing of aggregates (Or, 1996), we inferred a volumetric strain (non-compacted vs. compacted) of roughly 3%. When considering an aggregate radius of R = 5 mm (e.g., Márquez et al., 2004), the corresponding contact areas were 2.06, 0.70, and 2.16 mm 2 (see Equation 11). The average confining pressures P a (Equation 12) were about 1.23 and 0.21 MPa for compacted and non-compacted ley, respectively, with corresponding standard deviations of 0.5 and 0.1 MPa.

Pedophysical predictions of seismic velocities
All soil covers and treatments display gradual increases in seismic velocities during drying periods associated with increasing matric suction and decreasing bulk density (Equation 9) (Figures 5a, 5b, and 5c). Conversely, seismic velocities decrease sharply after rain events which, in turn, are associated with a rapid increase in water content and bulk density, and a decrease in matric suction. P-wave velocities are higher for the compacted soils (see Figure 4a). In particular, the highest P-wave velocities were measured in the compacted bare soil treatment, and the lowest were measured in the noncompacted ley treatment. Considering periods with overlapping water content values, the P-wave velocities of compacted ley are, on average, 31% higher than non-compacted ley and 1.7% lower than compacted bare soil. At water content close to field capacity (θ ∼ 0.35 cm 3 cm −3 ), the P-wave velocities of compacted treatments were higher (41.1%) than for the noncompacted treatment. These results provide evidence concerning the effects of persistent soil compaction on the seismic data. For saturation values close to 1, the P-wave velocity is expected to increase sharply in response to a rapid increase of the bulk modulus of the effective fluid (Equation 14). However, our measurements never reached such high saturation levels. For this reason, we only observe a decreasing trend in seismic velocities with increasing water content due to a decreasing suction stress and an increasing bulk density.
Our pedophysical predictions of P-wave velocities reproduce the main trends in the dynamics of the measured data such as those related to large rain events or drying periods (Figures 5 and 6). We evaluate the capacity of the model for reproducing the trend and fitting the observations with correlation coefficients and the weighted RMSE (WRMSE) values between observed and modeled P-wave velocities for each soil treatment. The correlation coefficients between predicted and measured P-wave velocities were r = .71 for compacted ley, r = .95 for noncompacted ley and r = .87 for compacted bare soil. The corresponding WRMSE were 1.6, 1.08, and 1.3, assuming relative data errors of 5%. The slope of the seismic velocity-water content curve is controlled by the van Genuchten exponent n m , which is the same for all treatments. As our model does not consider hysteresis in the water retention curve, a unique value of seismic velocity is predicted for each water content (see Figures 6a, 6b, and 6c).

DISCUSSION
This study presents a method for quantifying the compaction state of a soil from its seismic response. Particularly, we present a pedophysical model for interpretation of seismic measurements in terms of soil aggregate contacts. The model was tested in a field experiment where we monitored seismic properties and water content, associated with rain events and subsequent drying. We compared seismic velocities of different soil plots at the same water content (or saturation) over a wide range of values. Consequently, effects of soil compaction Our results suggest that compacted soil structure has not yet recovered from the effects of compaction as indicated by higher P-wave velocities measured in plots that were compacted 5 yr prior to our measurements. For a given water content, the seismic velocities are similar in compacted ley and compacted bare soil (only 1.7% difference, on average), suggesting that the presence of ley has not had a significant effect on the recovery of mechanical properties towards noncompacted conditions. This agrees with on-site point measurements of dry bulk density (soil cores sampled at 10-cm depth) and penetration resistance (average from 5-cm to 15cm depth) at the SSO presented in Figure 7. Figure 7 includes observations in the non-compacted bare soil, where seismic monitoring was not carried out. Figures 7c and 7d present the ratios of compacted vs. non-compacted dry bulk density and penetration resistance data by soil cover measured at the same time after compaction. These relative values provide insights of how mechanical properties recover from compaction with respect to the soil cover and suggest that (a) the soil cover does not have a strong influence on the recovery of soil mechanical properties towards non-compacted states, and that (b) soil mechanical properties still show a significant imprint of soil compaction 4 yr after the compaction event.
Five years after compaction, using our pedophysical model and the retrieved seismic data, we inferred contact areas that are 2.9 larger for compacted than for non-compacted soils. This corresponds to a seismic-inferred volumetric strain of 3% for compacted soils. In this context, Figure 8 presents F I G U R E 7 (a) Dry bulk density data as a function of the time since the compaction event for compacted ley, non-compacted ley, compacted bare soil, and non-compacted bare soil. (b) Penetration resistance as a function of the time since the compaction event for compacted ley, non-compacted ley, compacted bare soil, and non-compacted bare soil. (c) Relative dry bulk density (DBD) as a function of years since compaction for ley and bare soil. (d) Relative penetration resistance (PR) as a function of years after compaction for ley and bare soil. Error bars correspond to the standard error of the measured dry bulk densities and standard deviation of measured penetration resistance. The relative values were calculated as the ratio between the compacted treatment and its corresponding non-compacted treatment measured at the same time. For comparison, note that the seismic monitoring results were acquired 5 yr after the compaction event F I G U R E 8 Comparison of soil volumetric strain derived from dry bulk density data (ε ρ ) as a function of the time since the compaction event and the volumetric strain estimated from inferences of axial strains obtained from the pedophysical model and seismic data (ε M ). The volumetric strains from bulk densities and pedophysical model were obtained by comparing compacted ley (CL) and compacted bare soil (CB) with non-compacted ley at the same time after the experimental compaction event long-term volumetric strains derived from dry bulk density data measured in the compacted ley and compacted bare soil and data measured in their corresponding non-compacted treatments at the same time after compaction. For further information regarding the volumetric strain calculations, we refer the reader to Supplemental Section S3. The volumetric strains resulting from pedophysical modeling 5 yr after compaction (3%) fit nicely the post-compaction decreasing trend observed in volumetric strains computed from dry bulk density data (7% for 4 yr after compaction and 1.5% for 6 yr after compaction). Our model considers the aggregates as elastic spheres, for which the confining pressure P a (1.2 MPa for compacted soils) is required to maintain the inferred contact areas. This pressure should not be confused with the stress applied during visco-elastic deformation of the aggregates during compaction, which was much lower (150 KPa at 10cm depth, see Keller et al., 2017).
Our results show that combined seismic monitoring and pedophysical modeling can provide information related to soil compaction. In well-controlled field experiments such as the SSO, our framework of combined measurements and modeling could be used to evaluate the evolution of the mechanical status of soils at the plot scale by studying the long-term evolution of inferred aggregate contact area. Due to the low frequencies used, the small scale of the experiment and the target soil volume (first soil layer), our approach does not consider vertical variations of the overburden pressure and rather offers a simplified description of the effective compaction properties (i.e., contact area) of the first soil layer. Still, this already provides useful insights on the compaction state of soils. A detailed one-dimensional analysis of soil compaction would require the water content at all depths. Soil-structure-based concepts might be further developed to advance other applications of seismic monitoring.

CONCLUSIONS
We monitored seismic signatures of soils with different compaction treatments and covers at a controlled experimental field site, in which a compaction event took place 5 yr prior to the presented seismic monitoring. We found that P-waves carried a strong imprint of soil compaction. Furthermore, no significant differences were observed between soil cover (ley vs. bare soil), suggesting that the ley did not play a crucial role on the recovery of soil mechanical properties. Significantly higher P-wave velocities observed in the compacted plots suggested that the soils are still appreciably affected by the compaction event, thereby, confirming that soil structure recovery is a slow process. The seismic velocities were well reproduced with a newly proposed dual-domain pedophysical model accounting for soil plastic deformation due to soil compaction events. The model inferred contact areas between aggregates that are 2.9 times larger for compacted than for non-compacted soils. Based on our results, we suggest that seismic methods are suitable for characterizing soil structure, partly by offering a link to soil hydration status and partly by providing a direct link to soil mechanical properties to which other geophysical methods do not respond.

A C K N O W L E D G M E N T S
The geophysical data used in this study are available at Romero-Ruiz et al. (2020) (see https://doi.org/10.5281/ zenodo.3906467). The authors are grateful to Simón Lissa (University of Lausanne) for his help in the field monitoring campaigns and discussions on data processing. Nicolás Barbosa (University of Lausanne) is thanked for useful discussions on rock-physics models in the early stages of this work. Viktor Stadelmann, Valerio Volpe, and Selina Lutz (Agroscope) are thanked for their help with logistics at the SSO.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.