The research area located in the Province of Quebec, Canada (37°09′–36°42′ N; 38°48′–39°12′ E). We selected 49 potato farm sites where crop sequences included potato (Solanum tuberosum L.), maize (Zea mays L.), soybean [Glycine max (L.) Merr.], wheat (Triticum aestivum L.), alfalfa (Medicago sativa L.), barley (Hordeum vulgare L.), and canola (Brassica napus L.). Soils are classified as Aquents (Entisols), Orthods (Spodosols) and Udepts (Inceptisols) in the USDA soil classification system (Soil Survey Staff, 2014). In June and July 2014, 97 large (>20 kg) samples were collected from 49 horizon A and 48 horizon B to conduct Proctor tests (ASTM D1557, 2009). Sampling depth varied from site to site with the development of genetic horizons A and B. Average sampling depth for horizon A was 10.5 cm in the range of 4–18 cm. Sampling depth for horizon B was 17–48 cm, averaging 33.4 cm. Soil samples were air dried and then stored at room temperature. Smaller samples were collected in the center of each horizon using the cylinder method (Blake and Hartge, 1986).

Soil physical properties included gravimetric water content (Topp et al., 2007), BD (Blake and Hartge, 1986), particle-size distribution, MBD (ASTM D1557, 2009), CWC (ASTM D1557, 2009) expressed on an oven-dried (105°C) basis, and DC calculated as the ratio of BD to MBD. Proportions of coarse sand (0.50–2.00 mm), medium sand (0.25–0.50 mm), and fine sand (0.05–0.25 mm) were determined by dry sieving; the proportions of silt (0.002–0.05 mm) and clay (< 0.002 mm) were determined using a modified hydrometer method (Kroetsch and Wang, 2007). Soil chemical analyses were conducted using < 2.00 mm sieved samples. C and N were quantified by dry combustion (Leco-CNS). Si, Al, Fe, Mn, Mg, and Ca were extracted using the acid ammonium oxalate method (Courchesne and Turmel, 2007) and then quantified by the inductively coupled plasma (ICP) technique. Descriptive statistics of soil data are presented as raw data in Supplementary Table 1.

Isometric log ratios are orthogonal projections of compositional data arranged in binary subsets of components displayed in a sequential binary partition (SBP) (Egozcue et al., 2003). As there are D-1 degrees of freedom in compositional vectors (Aitchison and Greenacre, 2002), log ratio transformation techniques consist of D-1 ILR (Egozcue et al., 2003).

For the 14 soil components, there are 13 ILR variables designed to represent the BD and DC field (F) files, on the one hand, and the MBD and CWC Proctor (P) files on the other (Figures 1, 2). Given that soil gravimetric water content relates to BD and CWC relates to MBD, F1 was set as the balance between soil water content and solid components, while P1 was the balance between CWC for the Proctor test and the solid components. The F2 and P2 balances contrasted soil organic matter with mineral components, showing the functional role of organic matter in soil aggregation. Balances were further elaborated by associating other groups of particles. The F4 and P4 balances contrasted mineral soil particles with mineral cementing agents.

ILR is a normalized ratio between geometric means of two subsets of functional components (tagged “+” for parts in the numerator and “−” for parts in the denominator), calculated as follows (Egozcue and Pawlowsky-Glahn, 2006):

where i = 1 to D-1, ILR_i is the ith isometric log ratio in the ith row of the SBP between g(ci+) and g(ci-), geometric means across components, and ni+ and ni- are numbers of components labeled “+1” and “−1.”

Numerical analyses were performed in the R statistical computing environment using the compositions package (van den Boogaart et al., 2014) to transform compositional data into isometric log ratios, the dplyr package (Wickham and Francois, 2015) for general data manipulation, the nlme package (Filzmoser and Gschwandtner, 2015) to develop the linear mixed model, and the pls package (Revelle, 2014) for principal component analysis (PCA). We conducted PCA across (1) sampling depths and the 13 F balances to synthesize factors affecting BD and DC, and (2) the 13 P balances and 12 P balances to identify the components most closely related to MBD and CWC. Correlation analyses were conducted between the selected principal components (PCs) and compaction indexes (BD, DC, MBD, and CWC) using site scores (49 sites for the horizon A and 48 sites for the horizon B) of selected PCs.

Both horizons provided a gradient of soil properties to model BD and DC. The 13 F balances, along with sampling depth in both horizons, were used to predict BD and DC using the linear mixed-effects (LME) model as follows:

where Y is the soil compaction index (BD or DC), X is the fixed effects including 13 ILR balances, and sampling depth, β is the fixed-effects vector, Z is the random-effects matrix, b is the random-effects vector, and ε is the observational error vector. The site was considered as random effect. To model MBD and CWC, horizons A and B were separated as different entities requiring specific diagnoses for planning farm operations using regression analysis for each horizon. Prediction performance was evaluated using the Akaike information criterion (AIC) and the coefficient of determination (R²). The AIC value is used to compare and classify multiple competing models and to estimate which is the closest to the “real” process underlying the biological phenomenon under study (Burnham and Anderson, 2003; Burnham et al., 2011; Symonds and Moussalli, 2011). The coefficient of determination is the proportion of the variation that can be explained by the set of predictor variables. Mean prediction error (MPE) and root mean squared error (RMSE) were used to measure the reliability of the models, as follows:

where σ_i and ρ_i are the observed and predicted dependent variables for the ith measurement, respectively, and n is the number of observations. The MPE indicates an average underestimation (positive bias) or overestimation (negative bias) of dependent variables. For a good prediction, R² should be as large as possible, and the AIC value, MPE, RMSE should be as small as possible (Benites et al., 2007).

The first four PCs explained 71.3% of the total variation in the variables included in PCA (Table 1). All PCs showed significantly negative correlations with BD, whereas PC1, PC3, and PC4 were significantly and negatively correlated with DC. However, DC is a more useful compaction index than BD because DC corrects BD for MBD, an intrinsic soil property reflecting the maximum effect of machines on the degradation of soil physical quality.

The value and sign of the loadings and correlation coefficients are a measure of the relationships between the original variables and soil compaction indexes. The interpretation of loadings in Table 1 is straightforward. If loadings and correlation coefficients have the same sign, the relationship between ILR and compaction index is positive; otherwise, it is negative. ILRs are noted as [+1 or numerator group−1 or denominator group], hence, larger values in the denominator lead to more negative numbers in the log scale, and vice versa. For example, if Al loads more on F12 = [AlFe], the [AlFe] balance increases. If Fe loads more, the [AlFe] balance decreases, reversing the relationship between the [AlFe] balance and compaction indexes.

The effects of PC1 and PC3 on DC were highly significant (P < 0.01). The largest loadings in PC1 were F3, F9, F11, F12, and F13, indicating that DC was positively related to the mineral cementing agents, especially Si, Fe, Al, and Ca, and positively related to F3, hence DC was larger where the C/N ratio was higher. PC3 was primarily influenced by sampling depth, F2 and F8. DC was significantly higher in horizon B (Supplementary Table 1). F2, representing the balance between organic and mineral components, was negatively related to DC. F8, the ratio between coarse and medium-size sand particles, was negatively related to DC. The effects of PC4 on DC were also significant (P < 0.05), being mainly driven by F5, the [CS,MS,FS,siltclay] balance, where larger particles tended to increase DC compared to clay.

In contrast with DC, BD was significantly (P < 0.01) related to PC2. PC2 was largely influenced by F1, F4, F6, and F7. F1 was negatively related to BD. F4, reflecting the balance between mineral soil particles and mineral cementing agents, was positively related to BD. Where mineral cementing agents loaded more, BD was lower. F6 and F7 represented coarser particles in the soil. The higher the proportion of coarser particles, the higher BD was. BD was influenced significantly, but to a lesser extent, by PC1 and PC4. Consistent with DC, BD was also positively related to F5, F9, F11, F12, F13, and sampling depth.

The first four components and their loadings for MBD are presented in Table 2. The first four PCs explained 77.5 and 72.7% of total variation in all variables in horizons A and B respectively. However, the fourth PC was not significantly related to MBD in either horizon. Mineral cementing agents loaded most on PC1. The MBD was related negatively to P9, P11, P12, and P13 balances in both horizons. In PC2, the higher proportion of soil water (P1) and soil organic matter (P2) tended to reduce MBD in horizon A, whereas the higher proportion of coarser particles (P5–P7) tended to increase MBD in both horizons. In PC3, the higher proportion of coarser particles than of medium-size particles (P8) tended to reduce MBD in horizon A, whereas the higher proportion of clay than of coarser particles (P5) tended to increase MBD in horizon B.

The first four components and their loadings for CWC are presented in Table 3. The first four PCs explained 78.0 and 76.4% of total variation in all variables in horizons A and B respectively. The four PCs were significantly related to CWC in horizon A, and only the third PC was significantly related in horizon B. The CWC was related positively to P9, P11, P12, and P13 balances involving mineral cementing agents that loaded most on PC1 of horizon A. PC2 in horizon A revealed that increasing the proportion of organic matter compared to mineral solids (P2) increased CWC, whereas coarser particles (P6 and P7) tended to reduce CWC. PC3 in horizon A displayed a relatively higher content of coarse sand than of medium sand (P8), leading to a larger CWC, whereas a higher proportion of clay than of coarser particles (P5) tended to increase CWC in horizon B. Since P10 contributed most to PC4 of horizon A, a higher proportion of sesquioxides (Al, Fe, Mn) than of divalent cations (Ca, Mg) was positively related to CWC in horizon A.

Model coefficients and their significance are presented in Table 4. The F2, F4, F6, F7, F9, and F12 balances, meaning organic matter content, particle-size distribution and mineral cementing agents, significantly influenced BD prediction. The F1, F4, F5, F9, F10, and F12 balances, meaning soil gravimetric water content, particle-size distribution and mineral cementing agents, significantly influenced DC prediction. The MPE was 0.008 and 0.002 for BD and DC, respectively, indicating some underestimation of both BD and DC (Figure 3).

As shown by significant differences in BD and DC between horizons A and B, the sub-layers were more compacted than the arable layer for the pressure generated by loading the soil increases BD values deeper in the soil (Tranter et al., 2007). Studies have concluded that a small part of the variability in BD can be ascribed to sampling depth (Calhoun et al., 2001; De Vos et al., 2005; Heuscher et al., 2005). In the present study, sampling depth was not an independent variable as it is in most BD and DC prediction models (Reichert et al., 2009), because sampling depth is confounded with soil compositions.

Previous research reported that the degree of compaction mainly depends on soil moisture, texture and organic matter content (Jones et al., 2003; Hamza and Anderson, 2005; Dexter et al., 2008). BD has been found to decrease with higher clay and clay-plus-silt contents which is in agreement with our result (Kaur et al., 2002; Benites et al., 2007; Reichert et al., 2009), whereas the effect of medium and coarse sand on BD varied among soils, regions and horizons (Calhoun et al., 2001; Kaur et al., 2002; De Vos et al., 2005). Medium sand appeared to increase BD compared to coarse sand in our coarse-textured soils, in agreement with previous studies (Suuster et al., 2011). The effect of clay on DC was found to be negatively related to clay content but must also depend on organic matter content (da Silva et al., 1997).

It has been frequently reported that the BD of high-C soils can be explained by the variation in C content, whereas soil texture has a major effect on the BD of low-C soils (Manrique and Jones, 1991; Kaur et al., 2002). Total N showed a negative relationship with BD and DC (Benites et al., 2007). In the present study, organic matter content also tended to decrease the BD and DC of coarse-textured soils, whereas the [CN] balance was related positively to BD and DC. The [CN] balance decreases from lighter particulate organic matter to heavier fractions associated with soil mineral particles in agricultural soils (Yang et al., 2012). Total C thus provided partial information on the real contribution of organic matter to soil resistance to compaction.

Manrique and Jones (1991) suggested that texture and soil properties other than organic C play a more significant role in controlling BD deeper in the soil profile. A large contribution of water content to subsoil compaction is expected primarily on coarse- and medium-textured soils and heavy-textured soils with high shrink-swell potential (Suuster et al., 2011). The [soil water contentorganic matter, mineral soil particle, mineral cementing agent] balance increased DC, facilitating particle rearrangement, whereas soil gravimetric water content reduced MBD as the soil liquid phase is of low density and is incompressible.

Conversely, although it represents just a small part of total C, dissolved organic carbon (DOC) may contribute to the formation of naturally compacted layers (Sanborn et al., 2011). The main sources of organic C in subsoils are DOC, plant roots and root exudates, and organic particulates transported from soil surface (Rumpel and Kögel-Knabner, 2011). Organic C is a complex key soil quality index characterized by two major biochemical pools (Andrén and Kätterer, 1997) and several size fractions (Six et al., 2002; Stewart et al., 2008; Tong et al., 2014). Hence, total C should be further split into several C fractions (i.e., light fraction organic carbon, dissolved organic carbon and particulate organic carbon) to determine their respective significance in BD and DC models.

Felton and Ali (1992) found that the addition of organic matter increased soil porosity and water retention and reduced MBD as determined by the Proctor test. The effectiveness of organic matter depended on soil texture and organic matter quality (Zhang et al., 1997). Organic matter may increase in significance in cases where living and dead roots provide a more filamentous network across the soil profile (Soane, 1990). In the present study, the balance that included organic matter content showed negligible effect, probably because the range of soil organic matter contents was relatively narrow.

Clay content may lead to lower MBD (Smith et al., 1997) or show little effect (Aragón et al., 2000; Ball et al., 2000). The relationship between MBD and clay-plus-silt is quadratic for specific ranges of soil texture (Nhantumbo and Cambule, 2006; Mujtaba et al., 2014). Moolman and Weber (1978) reported that increasing evenness of particle-size distribution resulted in a higher MBD value, indicating the need to include particle-size balances such as P5 and P6 in MBD models. In contrast, Van Der Watt (1969) analyzed soils where MBD could be well-predicted from coarse sand (0.5–2.0 mm) only. Well-graded sands have a higher MBD or lower CWC compared to poorly graded sands (Guerrero, 2004; Mujtaba et al., 2014).

There is usually a close relationship between CWC and soil texture or organic C (Wagner et al., 1994; Aragón et al., 2000), but no such relationship was observed in the present study, apparently due to the narrow range of organic C contents. However, clay content was found to be linearly and positively related to CWC. The CWC was poorly predictable in horizon B. Mujtaba et al. (2014) and Guerrero (2004) found that MBD and CWC were poorly related in sandy soils. When water is added to dry soil, particles absorb a film of water. A certain amount of added water thickens the water film, allowing the soil particles to slide over each other in a process known as lubrication, which is triggered by fine-grain particles (Ishibashi and Hazarika, 2010). Bruand et al. (2005) concluded that particle-size distribution and the mineralogy of silt and clay particles associated with sand can lead to variations in physical properties, e.g., water retention, resistance to penetration and permeability of sandy soil. Zhao et al. (2008) reported that CWC was closely related to the liquid and plastic limits that integrate several soil properties such as particle-size distribution, organic matter content, and clay mineralogy. Hence, CWC could be predicted more accurately using liquid and plastic limits (Soane et al., 1972).

Mineral binding agents must have positive effect on maintaining soil structure and resisting to soil compaction (Lal and Shukla, 2004; Bronick and Lal, 2005). Where the average content of Al, Fe, Mn, Mg and Ca loaded more, BD and DC values were less. The positive relationship between F11, F12 with BD and DC appeared to reflect the effect of Al, Fe, and Mn on soil mass rather than on soil strength, as soil particle density increased with content of soil oxides, agreeing with McKeague and Sprout (1975). BD and DC were positively related to F13 ([CaMg]). Increasing soil aggregation and structural stability can be observed in soil with high Ca compared to Mg attributed to the thinner hydration radius of Ca and its effect on clay flocculation (Favaretto et al., 2006).

In general, only soil texture and organic C are used to predict MBD (Aragón et al., 2000; Nhantumbo and Cambule, 2006). We showed the importance of mineral cementing agents in predicting MBD of agricultural soils, as found by Zhao et al. (2008) for forest soils in British Columbia, Canada. The Si, Al, Fe and Ca oxides increased CWC but reduced MBD. Zhao et al. (2008) found that Al and Fe oxides were positively related to CWC, but were negatively related to MBD. Concentrations of Si, Mn, Al, Fe, Ca and Mg oxides are positively related to CWC as such compounds are hydrophilic (Tschapek, 1984). The negative relationship between MBD and mineral cementing agents in our study reflected the enhanced soil strength due to the presence of soil oxides.

Mechanical, biological and chemical means can be implemented alone or in combination to enhance soil resilience to compaction. The moisture content of coarse-textured soils should be < 50 g kg⁻¹, well below the median CWC of 205 g kg⁻¹ in horizon A and 144 g kg⁻¹ in horizon B, to achieve fracturing and shattering of the subsoil while avoiding damage to the soil due to compaction and problems with excessive draft or dust (Bannan and Wrigley, 2013). It has been reported (Goldsmith et al., 2001) that plant growth is not hampered during beneficial slope stabilization for engineering works where the degree of compaction is 0.80–0.85. For the median MBD values of 1.59 g cm⁻³ and 1.54 g cm⁻³ obtained through our prediction model for the horizon A and horizon B, it means BD of 1.27–1.35 and 1.23–1.31 for the horizon A and horizon B for the median conditions of our coarse-textured soil.

Rippers can be designed to suit specific soil conditions (Godwin, 2007) and combined with cropping systems and tillage practices to avoid coalescence (Lampurlanés and Cantero-Martinez, 2003; Reintam et al., 2008). Cover crops such as ryegrass can enhance soil quality through root channeling and organic matter additions (Darby et al., 2014; McNally et al., 2015). The LME models can assess the requirement for amendments to rebalance soil compositions to increase soil resilience to compaction. Soil amendments are commonly applied in the form of animal manure, gypsum (SO₄.2H₂O), lime, and municipal waste from waste water treatment plants providing Al and Fe.

Gypsum can reduce (1) the harmful effect of Mg on the soil structure of fine-textured soils resulting from the higher ion hydration radius of Mg compared to Ca (Favaretto et al., 2006) and (2) the subsoil Al toxicity in acid coarse-textured soils resulting from Al neutralization by the sulfate ion (Noble et al., 1988; Sumner, 1993). Gypsum increased maize yield by 29–50% on sandy loams (Toma et al., 1999), not only modifying the soil [CaMg] balance but also improving plant rootability. The beneficial effects of gypsum combined with ripping may last up to 16 years (Toma et al., 1999).

There are also ways to alter the texture of surface soils, such as partially mixing upper and lower layers if other soil quality attributes, such as organic matter content, are not unduly affected. Soil texture can also be modified to increase the water storage capacity of coarse-textured soils by applying to deficient soils the residual fine soil materials adhering to potato tubers, collected in the fall before tuber storage.

The authors declare that this study received funding from Cultures Dolbec Inc.; Groupe Gosselin FG Inc.; Agriparmentier Inc.; Prochamps Inc.; Ferme Daniel Bolduc et Fils Inc. The funders were not involved in the study design or collection, analysis, or interpretation of the data.