A Unified Model for Soil Shrinkage, Subsidence,and Cracking

Many clayey soils shrink when drying and swell when wetted. Such soils—often classified as Vertisols or vertic intergrades—are found across the globe, including within numerous agricultural and urban regions. Vertic soils include significant variation in physical properties, such as bulk density (Peng and Horn, 2007), pore size distribution (Kutílek, 1996), and field-saturated hydraulic conductivity (Messing and Jarvis, 1990; Jabro, 1996; Lin et al., 1998; Chertkov and Ravina, 2000; Chertkov, 2002; Das Gupta et al., 2006) with changes in the soil water content. During periods of drying, individual clay particles shed hydration layers, causing compaction of the soil aggregates. This process results in subsidence and cracking of the soil (Kutílek, 1996; Chertkov et al., 2004).

Numerous conceptual and mathematical models have been developed to describe the shrinkage of soil aggregates. Such models identify up to four distinct shrinkage phases, including: (i) structural shrinkage; (ii) proportional (i.e., basic or normal) shrinkage; (iii) residual shrinkage; and (iv) zero shrinkage (Peng and Horn, 2013). To date, however, there has been limited success in extending soil shrinkage curves to predict or describe field conditions. Indeed, while field cracks have been often characterized using in situ measurements of either crack volumes (Návar et al., 2002; Abou Najm et al., 2010; Stewart et al., 2012b, 2014) or vertical subsidence (Arnold et al., 2005; te Brake et al., 2013), these observations have not resulted in models that can predict dynamic variations in field crack volumes and/or dimensions. A recent exception provided a set of dual-permeability models that uses aggregate-scale shrinkage behavior to determine the dynamic porosity and/or flow-weighted distribution of the fracture vs. aggregate domains (Coppola et al., 2012, 2015). However, these models were derived for the specific case where horizontal cracking is predominant, which ignores subsidence effects. At the same time, the models require independent parameterization of the aggregate, fracture, and bulk domains, substantially increasing the challenge of obtaining the required values.

Shrinkage and cracking can be expressed as functions of water content (u) or suction head (h) (Salager et al., 2010; Coppola et al., 2015); however, these hysteretically linked properties have previously relied on independent sets of parameters for each relationship (i.e., the parameters used to describe the porosity–water content relationship are independent of those used to describe soil porosity and suction). This increases the effort required to parameterize these models and may be a reason why relatively few studies have explored the relationship between soil suction and soil shrinkage.

To provide a consistent linkage between soil aggregate shrinkage (as measured in the laboratory) and soil cracking and subsidence (as observed in the field), we propose a set of equations that describe the three major porosity domains that can exist in shrink–swell soils (i.e., aggregate, cracking, and subsidence). These porosity domains are described using parsimonious equations parameterized from laboratory and field measurements and can be used to predict how porosity shifts among these domains as a function of soil water content and/or soil suction. We then utilize these equations to develop expressions for the idealized geometries of pores within each domain.

Naturally occurring shrink–swell (vertic) soils are characterized by crack networks that extend from the soil surface to a depth at which either the soil texture changes, the soil water content level remains sufficiently high (thus preventing soil shrinkage), or the overburden pressure becomes large enough to inhibit soil swelling and shrinking (Jones and Jefferson, 2012). For this derivation, we define a control volume, ∀, that spans this “active” zone from the bottom of the lowest cracking layer to the maximum height of the soil surface, with length and width selected to be larger than the size of an individual pedon. We take erosion and other mass movements to be insignificant, and we assume that the total mass of solid particles within this volume (m_s) is constant. We focus on near-surface soils, avoiding complications that can arise from changing overburden pressures.

The total porosity within the control volume, ϕ_total, can be divided based on hydrological and morphological behaviors. In a typical shrink–swell soil, individual particles cluster together in aggregates, which in turn group together as a matrix. Assuming that the particles themselves are not porous, the smallest pores are those that occur between individual particles. The packing together of aggregates also forms pores that are typically larger and more connected than the “interparticle” pores; these “interaggregate” pores have also been defined as interpedal pores (Kutílek, 1996), capillary cracks (Chertkov and Ravina, 2001), interaggregate macropores (Abou Najm et al., 2010), or structural pores (Cabidoche and Ruy, 2001). Interparticle and interaggregate pores decrease in size as the soil dries (due to the loss of hydration layers in the case of the former and to increased suction in the case of the latter; Li et al., 2014), and increase in size as the soil wets. Here, we consider both of these porosity domains together and refer to the combination as the aggregate porosity, ϕ_aggr.

As the soil dries, typically, shrinkage cracks begin at the soil surface (Konrad and Ayad, 1997; Greco, 2002). As the soil profile dries, these cracks will extend vertically from the surface to some finite depth and often will become interconnected in the horizontal direction. As such, shrinkage cracks can be considered to demarcate individual pedons and have previously been referred to as border (Favre et al., 1997), inter-prism (Cabidoche and Guillaume, 1998), or planar (Bouma et al., 1977) cracks. Also during shrinkage the soil surface can vertically subside, which can be considered to be an independent porosity domain that forms external to the soil matrix. With these distinctions in mind, we partition ϕ_total into the following domains: (i) aggregate (ϕ_aggr); (ii) shrinkage cracks (ϕ_crack); and (iii) vertical subsidence (ϕ_sub) (Fig. 1). The distribution of porosity between these three domains varies as a function of soil water content (Kutílek, 1996; Fityus and Buzzi, 2009).

There is an ongoing debate about the importance of hysteresis in shrinkage and swelling processes. Some studies, primarily focused on laboratory samples, have contended that soil shrinkage and swelling are fundamentally irreversible processes (Braudeau and Mohtar, 2006; Peng and Horn, 2007; Chertkov, 2012), whereas other studies conducted in structured mineral soils have found that repeated shrinking and swelling cycles produce an equilibrium shrinkage curve with well-defined endpoints (Ring, 1966; Tripathy et al., 2002). Thus, while recognizing that hysteresis in shrinking–swelling processes may be an important consideration under certain conditions, we do not attempt to capture hysteretic behaviors at this stage. Likewise, the model considers the soils to be at thermodynamic equilibrium for any given water content, thus ignoring any time-dependent swelling processes that may exist (e.g., Braudeau and Mohtar, 2006).

To demonstrate the distinction among these three domains (aggregate, cracking, and subsidence), Fig. 2 provides an example of the porosity distribution at different water contents. The aggregate porosity is represented by the curve beginning at the minimum porosity ϕ_min and extending to the maximum porosity ϕ_max; as the aggregate porosity decreases (due to shrinkage), the shrinkage crack and/or subsidence porosities necessarily increase. This extra-aggregate (i.e., cracking and/or subsidence) porosity can be seen as the space between the soil shrinkage curve and ϕ_max. We use ϕ_ext to signify the extra-aggregate porosity of the soil.

Note that in this framework ϕ_max is equivalent to ϕ_total, and ϕ_min represents porosity that always remains within the interparticle + interaggregate domain and thus is never exchanged.

Thus, if we assume a constant shrinkage geometry factor, such as the common assumption of isotropic shrinkage, χ = 3 (Arnold et al., 2005; te Brake et al., 2013), we can reduce the number of fitting parameters by one. The derivation of Eq. [13] is presented in the appendix.

In this proposed set of equations, gravimetric water content (u) is clearly the critical variable that determines the distribution of soil porosity. However, six additional parameters are also required to fully describe these domains. Five of these parameters—u_max, ε, q, ϕ_min, and ϕ_max—can be directly estimated using the shrinkage curve for a given soil and thus are relatively easy to evaluate. The remaining parameter, ϕ_pedon, represents the division of extra-aggregate porosity into subdomains and can be estimated from the soil shrinkage curve if horizontal vs. vertical shrinkage is measured (Hallaire, 1984; Chertkov et al., 2004) or else can be calculated from any given shrinkage geometry factor, χ, using Eq. [13]. That said, future research is warranted to better understand and predict these parameters at the field scale.

Overall, the proposed set of equations is both parsimonious, in that only a handful of well-constrained parameters are required, and tractable, as those parameters can be directly determined from laboratory and/or field measurement.

Using the porosity models developed above, we can predict relative changes in the sizes of soil pores and cracks. First, we will idealize the total control volume ∀ as a collection of N soil blocks (Fig. 3a). Each soil block possesses dimensions of L_j by L_j by H (Fig. 3b) that on drying will shrink such that the geometry becomes (L_j− w_j) by (L_j− w_j) by (H − ΔH) (Fig. 3c), where w_j is the crack width of the jth soil block and ΔH is the subsidence. In the overall hierarchy of our framework, these soil blocks and cracks together constitute a pedon; this pedon then combines with any subsidence to comprise the total control volume. Note that this implies that all soil blocks will possess the same height (and experience the same amount of subsidence) for any given water content.

Within each idealized soil block, all aggregate porosity is represented by a total number of cylindrical soil pores M. Each cylindrical pore has some finite radius r_i. As the soil wets and dries, these pores all vary between some maximum radius, r_i,max (Fig. 3d), and some minimum radius, r_i,min. Here we make a couple of notes:

1. Even though many interparticle pores are likely to be irregular in shape and tortuous in structure, given the platy shape of typical clay particles, studies have demonstrated that structural interaggregate pores often have tubular shapes that can penetrate deep into the soil profile (Cabidoche and Ruy, 2001); moreover, these structural pores can be hydraulically dominant relative to the much smaller interparticle pores (Bouma et al., 1977). Therefore, for now we will use the cylindrical pore model for the aggregate domain.

2. Structural pores have traditionally been considered to have “stable” geometries, with less shrinkage and swelling than interparticle pores (Cabidoche and Guillaume, 1998). Thus, future model development may focus on understanding and integrating how the geometry of different sized pores vary with water content.

Returning to our idealized soil model, during the shrinkage process, any subsidence will cause the soil blocks to contract in the vertical dimension by ΔH such that the height of the soil blocks at any water content, U, will be H − ΔH. Also during shrinkage, a crack forms on the outside (i.e., border) of each soil block, with approximate dimensions of w_j by (2L_j− w_j) by (H − ΔH) (Fig. 3e). We assume that all shrinkage cracks have rectangular geometry.

Equation [21] thus signifies that, by using the porosity distribution parameters ε and q, we can calculate the relative width of any given soil crack for any given water content. Moreover, by measuring soil crack widths at a single (known) water content, we can predict the maximum widths of the cracks. Note that ε and q are considered to be macroscale constants in both space and time, whereas w_max can vary in space but not time (and water content u can vary in both time and space). Also note that this formulation implies that all crack widths change in proportion to one another, which causes the crack size distribution to maintain a constant shape and variance for all water contents u < u_max. Finally, in the case where information about the overall crack size distribution is not required or available, w_max can be considered to be a single effective macroscale parameter.

where h_b is known as the bubbling pressure and λ is a pore size distribution parameter.

Soil samples were collected at two study sites. The first was located near Corvallis, OR, in a 2- by 3-m field plot. The soil was classified as a Waldo silty clay loam, a fine, smectitic, mesic Fluvaquentic Vertic Endoaquoll with moderate to high shrink–swell potential (Knezevich, 1975). The soil was found to have 5% sand, 37% silt, and 58% clay (from hydrometer tests), giving it a clay texture. The second site was located near the community of Ninhue in south-central Chile, in a set of 3.5- by 11-m field plots. The near-surface (0–20-cm) soil was found to have 44% sand, 31% silt, and 25% clay, giving it a loam texture. Stewart et al. (2015) provided a detailed description of the site, procedures, and analyses.

Three intact soil clods collected from the Oregon site (Samples W6, W7, and W8) and three soil clods from the Chile site (Samples L5–000, L4–085 and L5–060) were analyzed for shrinkage and swelling using the imaging method of Stewart et al. (2012a). The Oregon samples were collected from a depth of ∼20 cm. The sampling depth (in centimeters) of the Chile samples is indicated by the three-digit sequence after the dash (e.g., L4–085 was collected within Plot L4 at a depth of 85 cm).

For the Oregon samples, shrinkage measurements were taken as the clods were allowed to air dry from a water content near field capacity (as presented by Stewart et al., 2012a). For the final shrinkage measurement, the cores were oven dried at 105°C for 24 h. Individual measurements of volume and water content were binned and averaged to fit the Tariq and Durnford (1993) soil shrinkage model. At the conclusion of the shrinkage measurements, three of the samples were rehydrated to measure their swelling. Sample W6 was hydrated by being slowly placed in 0.5 cm of standing water and allowed to wet by capillary rise. The volume was measured after 1 d and again after 2 d, with no change being observed after the first day. Sample W7 was hydrated by being placed in a sealed humidity chamber with relative humidity >99%. Volume measurements were taken every day until the water content level stabilized (determined as being when the sample weights and/or volumes did not increase between successive days). Sample W8 was hydrated by being slowly immersed in water and was then left submerged for 1 wk, after which time the volume was measured.

The Chile samples were initially wetted from oven-dry (105°C) conditions by being placed in a humidity chamber with relative humidity >99%. Once the water content stabilized (again, determined as being when the sample weights and/or volumes did not increase between successive days) the samples were allowed to re-dry. Sample volumes and weights were measured once or twice per day during the swelling and shrinking processes. Note that Samples L4–085 and L5–060 lost their structural integrity and fell apart during the shrinkage process, so only data from the swelling direction is included for those two samples.

During January of 2012, we performed an irrigation experiment at the Chile site, in which water was applied to the field plots during a 3-wk period and the soil hydrologic response was quantified (Stewart et al., 2015). As part of that experiment, we collected data on crack closure using surface-based digital imaging. The digital images were collected using a Pentax K-x dSLR camera mounted on a tripod and positioned approximately 0.6 m above the soil surface. Digital imaging was done within four of the plots; in each plot, three large representative cracks were delineated using a 50- by 50-cm quadrat. We collected images at each quadrat before, during, and at the conclusion of the irrigation experiment. At the same time, we also collected duplicate soil cores from each plot to determine the soil water content at the time of each image. The surface area of cracks within the quadrats were quantified by first converting the digital images to a gray scale and then renormalizing them using a custom MATLAB script such that cracks were black and soil was white. The number of black pixels were then counted and converted to an area. Then, to estimate representative (i.e., effective) crack widths, we assumed that each crack had a length L = 50 cm (the length of the imaging frame). To fit Eq. [21] to the data, we then scaled the widest observed crack width to U = 0, thus predicting w_max for each individual crack. A set of example crack images, before and after processing, is shown in Fig. 4.

To demonstrate the versatility of the proposed porosity model, we fit it to several sets of measurements collected in previous studies. First, Reeve and Hall (1978) measured shrinkage and water retention (i.e., the water content–soil suction relationship) for six clayey subsoil samples; we reanalyzed three of those samples (Wrye Bw horizon, Wrye Bg horizon, and Fladbury Bg1 horizon) using the proposed model. Next, Hallaire (1984) measured the shrinkage of soil cores taken from a clayey soil, while also quantifying subsidence using a displacement transducer and the cracking of the samples by filling the voids with small glass beads. Hallaire (1984) also estimated crack volumes in a field setting by collecting and analyzing images. Finally, Salager et al. (2010) varied the initial compaction of a clayey silty sand and measured the relationships among sample volume, water content, and suction as the samples were equilibrated in a pressure plate system. We chose three of the samples from Salager et al. (2010) to reevaluate for purposes of this study (e₀ = 1.01, e₀ = 0.86, and e₀ = 0.68, where e₀ is equivalent to e_max).

Nonetheless, it should be noted that Eq. [10] does not have the ability to describe the various shrinkage domains, and thus other models are recommended if the ability to separate individual shrinkage phases is necessary.

For the Oregon soil clods (W6, W7, and W8), Eq. [10] was compared with the four-phase shrinkage model of Tariq and Durnford (1993) using the parameters ε = 3.2 and q = 2 (Fig. 6a). Across the entire water content range, the two models differed by a maximum of 3.9%, with a mean difference of 2.0%. The second derivative κ(U) (from Eq. [27]) had a maximum value at U ≈ 0.6, thus identifying the approximate water content at which shrinkage transitions from structural to proportional. Shrinkage measurements for the Chile soil clods (Samples L5–000, L4–085, and L4–060) were also modeled using Eq. [10] (Fig. 6b). In this case, Eq. [10] was fit to the data using values of ε = 4 and q = 1.7.

Overall, these results suggest that Eq. [10] is a suitable approximation for describing the shrinkage curve of structured clay soils, while requiring only one variable (water content) and three other parameters (or five when presented in non-normalized terms).

As an example of how the proposed functions (Eq. [10], [12], and [14]) can accurately describe the various porosity domains of shrink–swell soils, we used them to model the results from a laboratory-based shrinking soil core experiment and the field-based crack study conducted by Hallaire (1984). The modeled curves were generated using values of ε = 4.5, q = 3.3, ϕ_max = 0.504, ϕ_min = 0.308, ϕ_pedon = 0.459, and u_max = 0.377. For the laboratory cores, Eq. [10] predicted Φ_aggr with a total RMSE of 0.024 (mean difference 0.018) compared with the observed values (Fig. 7a). Moreover, Eq. [12] predicted Φ_sub with a total RMSE of 0.030 (mean difference 0.026), and Eq. [14] predicted Φ_crack with a total RMSE of 0.020 (mean difference of 0.014) compared with the observed values. Thus, in this example we were able to match the observation of sample shrinkage, cracking, and subsidence within 3%, all while using the same set of five parameters. Using those same parameters, we were also able to model the field crack volumes that were measured in the same soil (Fig. 7b), albeit with an RMSE value of 0.063 and a mean difference of 5.2% compared with the measured crack volumes.

As another example, we compared predicted crack widths (using Eq. [21]) with actual crack widths (estimated using digital imaging) during an irrigation experiment at the Chile site (Fig. 8). Equation [21] was evaluated using values of ε = 4 and q = 1.7 (as predicted by the soil shrinkage curves quantified above). The model closely predicted crack widths as a function of soil water content, with some deviations observed for Cracks 2 and 6 toward the wet end of the curve (i.e., U < 0.6). This example again shows how soil shrinkage parameters estimated from laboratory samples can accurately predict field crack dynamics.

Although soil shrinkage is conveniently measured and modeled as a function of water content, in actuality the soil matric potential (or equivalently, soil suction) determines both water content and the degree of shrinkage. To demonstrate this relationship, we used the results from two previously published studies (Reeve and Hall, 1978; Salager et al., 2010) where these three variables were all measured (Fig. 9 and 10). The soil suction–water content relationship was modeled using Eq. [24], while the soil suction–shrinkage relationship was modeled using Eq. [25]. Specific parameters used to fit the curves are presented in Tables 1 and 2.

For the soil suction–water content relationship, Eq. [24] had a maximum RMSE of 0.038 compared with the Salager et al. (2010) data and 0.0095 compared with the Reeve and Hall (1978) data. For the soil suction–porosity relationship, Eq. [25] had a maximum RMSE of 0.0055 compared with the Salager et al. (2010) data and 0.0075 compared with the Reeve and Hall (1978) data. Thus, using only six total parameters (ε, q, α, n, ϕ_max, and ϕ_min), the two proposed equations fit both data sets with minimal error.

Note that we also predicted crack porosity, assuming isotropic shrinkage, as a function of soil suction for the Reeve and Hall (1978) data (Fig. 9). Interestingly, even though the three soil samples (Wrye Bw, Wrye Bg, and Fladbury) had distinct water retention curves and offset soil shrinkage curves, the predicted crack porosity was very similar for any given soil suction. This is probably an artifact of the limited range across which soil suction was measured (5–1500 kPa), as the Fladbury soil had the greatest overall amount of shrinkage between saturated and dry conditions.

In many clayey soils, porosity and other hydraulic properties vary with soil water content. While a number of approaches and solutions have been developed to account for nonequilibrium (i.e., preferential) flow in soils, there is still a need for accurate parsimonious models to predict how changes in soil water content affect the porosity distribution of the soil. In this study, we proposed a framework to describe the porosity distribution in shrink–swell clay soils, focusing on three porosity domains: (i) aggregate (ϕ_aggr); (ii) shrinkage cracks (ϕ_crack); and (iii) subsidence (ϕ_sub). The behavior of the aggregate domain can be understood through application of the soil shrinkage curve, for which we proposed a new expression that, when presented in normalized terms, requires only water content and two fitting parameters (or, alternatively, soil suction and four fitting parameters). This new soil shrinkage function is flexible and capable of describing many different shrinkage behaviors; it can also be readily integrated and differentiated, thus allowing estimations of shrinkage phase transitions.

More importantly, the soil shrinkage function can be used to predict soil cracking and subsidence using the same basic set of parameters, which represents a significant advancement from previous models. The shrinkage crack and subsidence domains can be collectively considered to form the “extra-aggregate” shrinkage porosity of the soil. By equating the decrease in aggregate porosity due to shrinkage with the increase in extra-aggregate porosity, our model can be used to predict the dynamic opening and closing of soil cracks and surface subsidence. This approach provides the framework for accurate and unified hydromechanical modeling of swelling soils.

Using data from previously published studies, we demonstrated the ability of our shrinkage model to accurately estimate changes in aggregate, crack, and subsidence volumes as well as crack widths. Given this additional capability and functionality, this new set of models should find broad application within the soil physics and hydrology communities.

This work was supported under National Science Foundation Award 0943682. Funding for this work was provided in part by the Virginia Agricultural Experiment Station and the Hatch Program of the National Institute of Food and Agriculture, USDA.

Extra-aggregate porosity—here defined as the porosity that forms outside of the soil matrix as a result of shrinkage—can be divided into subsidence (ϕ_sub) and cracking (ϕ_crack). The most common graphical depiction of soil shrinkage uses a cube that becomes more compact, resulting in an increase in extra-aggregate porosity (Bronswijk, 1990; Chertkov et al., 2004) (e.g., Fig. 3). It can be seen that subsidence is considered to be related to any change in the vertical direction (ΔH), while cracking is considered to be the remainder of the extra-aggregate volume, which may include cracks (of any orientation) within the soil matrix as well as those that are found outside of the soil matrix (also known as “inter-block” cracks).