# Vadose Zone Journal

- Copyright © by the Soil Science Society of America, Inc.

## Abstract

Many clayey soils shrink as they dry, causing a shift of porosity from inside to outside the soil aggregates and leading to the formation of shrinkage cracks and/or surface subsidence. During swelling, shrinkage cracks begin to seal and/or the soil surface rises. Previous models have focused on describing shrinkage at the aggregate level, with little success in predicting soil cracking and subsidence. To remedy this shortcoming, we provide a unified, physically based set of governing equations for these three pore domains (aggregates, cracks, and subsidence) and predict the porosity distribution among domains as a function of soil water content and minimal (up to six) additional parameters. Examples collected from a variety of soils show how these functions describe shrinkage of soil samples in the laboratory; quantify the relationships among soil suction, soil shrinkage, and water content using the same set of parameters; and predict sealing of soil cracks in the field. This approach provides the framework for accurate and unified hydromechanical modeling of swelling soils.

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.

## Theory

### Physical Considerations and Assumptions

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.

One of the challenges of working with shrinking soils is that several descriptions of the relationship between volume and water content are used. For example, the volume of a soil aggregate may be reported in terms of the void ratio *e* (volume of aggregate voids over volume of solids), specific volume *v* (volume of aggregate over mass of solids), or porosity ϕ (volume of aggregate voids over total volume). Likewise, water content can be expressed as a moisture ratio ϑ (volume of water over volume of solids) or gravimetric water content *u* (mass of water over mass of solids). However, if the absolute maximum and minimum of each variable are known, then all of these variables can be normalized to a single curve of Φ vs. *U*:where the subscripts *min* and *max* refer to the respective minimum and maximum values for each parameter. The following constraints are often assumed:where ρ_{w} is the density of water and ρ_{s} is the density of the solid particles. To allow easy comparison among different data sets and models, we have opted to present our model using both absolute and normalized terms.

### Derivation of Governing Equations

The presence or absence of a void at any point *p* within our control volume, ∀, is indicated byand the total porosity, ϕ_{total}, is defined as

Note that we are assuming that the total mass of solid particles (*m*_{s}) within this volume is constant with time and can be determined as

Using the definitions and assumptions outlined above, we apportion the porosity of a swelling soil as

In structured clay soils, repeated shrinking and swelling cycles can produce an equilibrium shrinkage curve with well-defined endpoints (Tripathy et al., 2002). The difference between the minimum point (ϕ_{min}) and maximum point (ϕ_{max}) of the equilibrium shrinkage curve can be considered to be the reversibly exchangeable porosity, ϕ_{rev}:or in normalized terms,where

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.

### Aggregate Porosity

Several approaches exist for modeling the soil shrinkage curve, including quantifying only certain phases of the overall curve (Giráldez et al., 1983; McGarry and Malafant, 1987), quantifying all phases with assumptions on the shape of each (Tariq and Durnford, 1993; Braudeau et al., 2004), or developing expressions based on the shape of the entire curve (Groenevelt and Grant, 2001; Peng and Horn, 2005; Leong and Wijaya, 2015). For our derivation, we utilize the latter approach and approximate the normalized soil shrinkage curve with a single continuous function:orwhere ε and *q* are fitting parameters. Note that the function proposed in Eq. [10] has only four free parameters (ε, *q*, ϕ_{max}, and ϕ_{min}), one fewer than the Groenevelt and Grant (2001) and Peng and Horn (2005) models. In theory, this makes parameter estimation easier, with reduced uncertainty, compared with competing models.

### Extra-aggregate Shrinkage Porosity

The combination of crack and subsidence porosity represents the extra-aggregate shrinkage porosity, Φ_{ext}. As the interparticle and interaggregate porosities decrease due to shrinkage, the extra-aggregate porosity increases. This relationship can be seen by substituting Eq. [10] into Eq. [9]:

### Subsidence Porosity

For purposes of this derivation, we approximate the soil subsidence as being proportional to the total extra-aggregate porosity:where ϕ_{pedon} is the total porosity attributable to the combination of cracks and aggregate porosities at *U* = 0. In other words, ϕ_{pedon} represents all porosity contained below the soil surface at the point of maximum subsidence.

Note that we can also express ϕ_{pedon} in terms of the shrinkage geometry factor χ (Bronswijk, 1990) as

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.

### Crack Porosity

Substituting Eq. [12] into Eq. [10] and rearranging allows us to describe the crack porosity as

### Required Parameters

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.

### Pore and Crack Geometries

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 *j*th 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:

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.

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 (2*L*_{j} *− w*_{j}) by (*H* − Δ*H*) (Fig. 3e). We assume that all shrinkage cracks have rectangular geometry.

We start with the aggregate porosity ϕ_{aggr}(*U*), which is calculated aswhere *V*_{v,aggr} is the volume of aggregate voids and *V*_{t} is the total volume (i.e., *L*^{2}*H*). Note that Δ*H* can be estimated from the subsidence porosity as Δ*H* = Δ*H*_{max}[(1 − *U*^{q})/(1 + ε*U*^{q})], where Δ*H*_{max} is the maximum subsidence (at *U* = 0).

As the soil wets and dries, the aggregate porosity will change as a function of soil water content. We can quantify this variation by substituting Eq. [10] into Eq. [15]:

Note that in the case where the height of the soil block is much greater than the change in height due to subsidence (i.e., *H* >> Δ*H*_{max}), Eq. [16] will simplify toThe shrinkage crack porosity ϕ_{crack}(*U*) can be written aswhere *w*_{j,max} is the maximum crack width corresponding to the *j*th soil block (when *U* = 0).

Equation [18] can be rearranged asSo long as 2*L*_{j} >> *w*_{j,max}, the crack width can be expressed in terms of *w*_{j,max} as

Moreover, if Δ*H*_{max} << *H*, Eq. [20] reduces to

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.

### Relation with Water Retention

Soil shrinkage and cracking can also be described in terms of the soil suction, *h*. Indeed, the relationships between shrinkage and suction, shrinkage and water content, and water content and suction have been combined to create three-dimensional soil water retention surfaces for various soils and/or compaction conditions (e.g., Salager et al., 2010). Due to the limited number of parameters and the overall simplicity of the proposed shrinkage model (i.e., Eq. [10]), it is straightforward to incorporate classic water retention models such as those proposed by van Genuchten (1980) and Brooks and Corey (1964). For example, the van Genuchten (1980) equation can be written aswhere α, *n*, and *m* are fitting parameters. Then, substituting Eq. [22] into Eq.[10] and using the condition that *m* = 1 − 1/*n*, we obtain a compaction model in terms of suction pressure:

Likewise, we can use the Brooks and Corey (1964) function:to express the soil shrinkage as

where *h*_{b} is known as the bubbling pressure and λ is a pore size distribution parameter.

Finally, we can also estimate the relationship between soil suction and crack porosity, e.g., by combining Eq. [14] and [22]:

## Methods

### Study Sites

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.

### Soil Shrinkage Parameters

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.

### Crack Widths

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.

### Shrinkage Measurements from Previous Studies

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*_{0} = 1.01, *e*_{0} = 0.86, and *e*_{0} = 0.68, where *e*_{0} is equivalent to *e*_{max}).

## Results and Discussion

### Soil Shrinkage Curve

Equation [10] is a highly flexible function capable of producing a variety of curve shapes (Fig. 5). Even though the fitting parameters *q* and ε lack true physical meaning, they each control certain aspects of the soil shrinkage curve. For example, the fitting parameter *q* can be considered to be related to the amount of residual or zero shrinkage, with small values of *q* (i.e., *q* < 1) reflecting soils that have little or no zero or residual shrinkage. The value of ε, on the other hand, is more reflective of the amount of structural shrinkage, with larger values of ε associated with soils that have greater amounts of structural shrinkage. Furthermore, as discussed by Groenevelt and Grant (2001), the second derivative of the soil shrinkage curve, κ(*U*), can be used to identify the approximate locations of the transition points between structural and proportional shrinkage, and again between proportional and residual shrinkage. In the case of Eq. [10], the structural–proportional shrinkage transition can be estimated as the maximum value of κ(*U*) [and the proportional–residual shrinkage transition can be estimated as the minimum value of κ(*U*)], where κ(*U*) is found by

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).

### Predicting Field Crack Volumes from Soil Shrinkage Measurements

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.

### Soil Suction–Water Content–Porosity Relationships

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.

## Conclusions

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.

## Acknowledgments

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.

## Appendix

### Incorporation of Shrinkage Geometry Factor into Model

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).

The shrinkage geometry factor, χ, was devised as a way to relate the relative change in total layer volume Δ*V*/*V* to the relative change in layer thickness Δ*H*/*H*, i.e.,where *V*_{t} is the total (initial) volume of the cube (i.e., *L*^{2}*H*).

The left-hand side of Eq. [A1] can be rewritten aswhere *V*_{v,aggr} is the volume of voids within the aggregates (matrix), *V*_{s} is the volume of solids within the matrix, ϕ_{aggr} is the porosity of the aggregate, and ϕ_{max} is the total porosity of the cube.

In the context of this cubic soil model, the subsidence porosity can be written aswhere *V*_{sub} is the volume of the void associated with subsidence (i.e., *L*_{j}^{2}Δ*H*).

Substituting Eq. [A2] and [A3] into [A1], we obtain

Evaluation of Eq. [A4] at *U* = 0 gives

We define ϕ_{pedon} as all porosity contained below the soil surface at the point of maximum subsidence:

## Footnotes

Reproducible Research—This article includes supplemental material online to permit readers to analyze the data in a manner similar to that presented and reproduce results.

- Received November 11, 2015.
- Accepted January 6, 2016.

This is an open access article distributed under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).