Effects of dynamic changes of desiccation cracks on preferential flow: experimental investigation and numerical modeling

Preferential flow induced by desiccation cracks (PF-DC) has been proven to be an important hydrological effect that could cause various geotechnical engineering and ecological environment problems. Investigation on the PF-DC remains a great challenge due to the soil shrinking–swelling behavior. This work presents an experimental and numerical study of the PF-DC considering the dynamic changes of desiccation cracks. A soil column test was conducted under wetting–drying cycles to investigate the dynamic changes of desiccation cracks and their hydrological response. The ratios between the crack area and soil matrix area (crack ratio), crack aperture and depth were measured. The soil water content, matrix suction and water drainage were monitored. A new dynamic dual-permeability preferential flow model (DPMDy) was developed, which includes physically consistent functions in describing the variation of both porosity and hydraulic conductivity in crack and matrix domains. Its performance was compared to the single-domain model (SDM) and rigid dual-permeability model (DPM) with fixed crack ratio and hydraulic conductivity. The experimental results showed that the maximum crack ratio and aperture decreased when the evaporation intensity was excessively raised. The self-closure phenomenon of cracks and increased surficial water content was observed during low-evaporation periods. The simulation results showed that the matrix evaporation modeled by the DPMDy is lower than that of the SDM and DPM, but its crack evaporation is the highest. Compared to the DPM, the DPMDy simulated a faster pressure head building-up process in the crack domain and higher water exchange rates from the crack to the matrix domain during rainfall. Using a fixed crack ratio in the DPM, whether it is the maximum or the average value from the experiment data, will overestimate the infiltration fluxes of PF-DC but underestimate its contribution to the matrix domain. In conclusion, the DPMDy better described the underlying physics involving crack evolution and hydrological response with respect to the SDM and DPM. Further improvement of the DPMDy should focus on the hysteresis effect of the soil water retention curve and soil deformation during wetting–drying cycles.

1 Introduction

Desiccation cracks are prevalent in clay-dominated soils due to water loss and often lead water to bypass the surface soil matrix and rapidly infiltrate into subsoil as preferential flow (Davidson, 1984; Weiler, 2005). Positively, the preferential flow induced by desiccation cracks (PF-DC) can promote the migration of farmland organic matter (Vervoort et al., 2003) and reduce surface runoff (Pei et al., 2020; Zhang et al., 2021a). Negatively, it also has proven to be an important hydrological mechanism that could lead to geotechnical engineering and ecological environment problems, such as dike and slope instability (Jamalinia et al., 2020; Zhang et al., 2021b), shallow landslides (Bogaard and Greco, 2015; Caris and Van Asch, 1991; Luo et al., 2021), groundwater pollution (Chaduvula et al., 2022; Chen et al., 2002; Mooney and Morris, 2008; Schlögl et al., 2022) and reduction of irrigation efficiency (Greve et al., 2010; Smith et al., 2005; Wang et al., 2018; Wang et al., 2022). Under the current background of frequent extreme flood–drought climate events, its negative effects will be more prominent (Tichavsky et al., 2019). Investigations on the PF-DC are of great significance in guiding scientific research and practical design in the above disciplines.

A unique characteristic of the desiccation cracks is their dynamic features, often causing instantaneous variation of crack proportion, depth and connectivity with moisture content. Previous efforts have attempted to reveal the effects of crack dynamics on the PF-DC through experimental studies, but most of them focused on the short-term wetting process and obtained only qualitative results, and debates remained. For instance, Favre et al. (1997) and Liu et al. (2003) stated that crack closure due to wetting can cause a significant reduction or even disappearances in the preferential flow. However, other studies found that the PF-DC also leads water to rapidly infiltrate into deep soil, even when desiccation cracks are nearly closed (Baram et al., 2012; Greve et al., 2010; Luo et al., 2021; Tuong et al., 1996; Sander and Gerke, 2007). Cheng et al. (2021) conducted a series of constant-head permeability tests with the hydraulic head gradient of 15 kPa. They stated that 4 % of surface crack ratio could be a critical value for determining whether desiccation cracks cause a significant increase in the infiltration rate or not. However, this value may vary with different soils, rainfall patterns and sample scales and thus lacks general applicability. Indeed, PF-DC has long-term and complex spatiotemporal variability due to crack dynamics during wetting–drying cycles. Therefore, short-term and small-scale infiltration tests (i.e., laboratory permeability tests) are not enough to reveal the complex hydrological process induced by PF-DC. Meanwhile, it is also difficult to quantitatively study PF-DC only through experiments. An improved understanding of the PF-DC combined with theoretical methods is also needed.

Regarding the theoretical methods, explicit crack models (EMs) (Hendrickx and Flury, 2001; Khan et al., 2017; Xie et al., 2020), dual-porosity models (DPoMs) (van Genuchten, 1980; van Genuchten and Wierenga, 1976) and dual-permeability models (DPMs) (Aguilar-López et al., 2020; Gerke and van Genuchten, 1993b, a) were developed to simulate preferential flow in cracked clay soils. EMs were constructed based on the single-domain (or single-permeability) framework, which requires definition of the details involving the geometry, spatial distribution and hydrological properties of each crack. Such a requirement may be conceptually correct but makes it difficult to simulate network-distributed desiccation cracks due to considerable computational burden (Aguilar-López et al., 2020). The DPoM and DPM concepts belong to the dual-domain framework that assumes the soil pore system can be represented as two overlapping interacting regions, one which represents the matrix domain with micropores and the other one which represents the crack domain with mesopores–macropores (Šimùnek et al., 2003). Those models represent the cracks in the soil as an implicit form which need not prescribe geometrical and spatial features of the desiccation cracks. The DPoM concept holds the simplifying stipulation that water only flows through the shrinkage cracks rather than the soil matrix, which is unrealistic in many cases. To remedy this shortcoming, a classical DPM was developed, where the water flow in soil matrix and crack domain was simulated using the Richards' equation (Aguilar-López et al., 2020; Coppola et al., 2012; Gerke and Maximilian Köhne, 2004; Gerke and van Genuchten, 1993a) or Green–Ampt model (Davidson, 1984; Stewart, 2019; Weiler, 2005) building on Darcy's law. However, some critics emerged that the Richards' equation building on the capillarity, not existing in large PF paths (e.g., tensile cracks and biological holes), is not suitable to simulate the PF (Larsbo and Jarvis, 2003; Nimmo, 2010, Nimmo et al., 2021). Consequently, some improved DPMs were developed, where water flow in the crack domain was simulated by the Navier–Stokes equation (Germann and Karlen, 2016; Nimmo, 2010), kinematic wave equation (Greco, 2002; Larsbo and Jarvis, 2003) and Poiseuille model (Lepore et al., 2009). Although these improved DPMs better captured the characteristics of the water flow in the crack domain, the classical DPM concept has still been widely accepted and used in simulating preferential flow in soils due to its easily available parameters, reasonably satisfactory prediction to the measurements and high computation efficiency (Jarvis et al., 2016). Most importantly, a recent numerical study conducted by Aguilar-López et al. (2020) proved that effective parameter selection in the DPMs can achieve similar modeling results to the EMs.

Nevertheless, classical DPMs often adopt the assumption that crack volume and hydrological properties remain constant in both time and space, which is unfeasible to capture the full dynamics of PF-DC. Some attempts have been made to incorporate the dynamic nature of desiccation cracks into DPM including the SWAP family of models, i.e., LEACHM, which simulates PF-DC using a shrinkage characteristic and water loss (Kroes et al., 2000) but neglects the water exchange process occurring at the interface between two domains. Such a process has widely been confirmed to be significant in cracked soils (Greve et al., 2010; Krisnanto et al., 2016; Tuong et al., 1996). A later modification of SWAP incorporated the aforementioned process but at the cost of neglecting shrink–swell behavior of soil. The VIMAC model developed by Greco (2002) solved previous problems but against the cost of inducing many parameters which are difficult to determine from experiments or measurements. Coppola et al. (2012, 2015) took another step forward to allow crack volume and/or hydrological properties to vary as a function of soil shrinkage. However, the relationship proposed in the model, an empirical natural logarithm function involving the suction head and crack proportion, is not directly transferable to other types of soil. Stewart et al. (2016b) deduced a shrinking–swelling model, with relatively clear physical meaning and high consistency, and recently incorporated it into a Green–Ampt-based DPM (Stewart, 2018). While an analytical solution was obtained, the intrinsic limitation of the Green–Ampt approach (i.e., hypothesis of the wetting front and request for a constant boundary condition) hindered the further application of this model in complicated scenarios.

The objective of this research was to investigate the PF-DC from the experimental perspective in combination with an effective modeling approach. Hence, a soil column test was conducted to investigate the dynamic changes of desiccation cracks and hydrological response. The variation of crack geometry, including crack ratio, width and depth, was measured. The soil moisture content, matrix suction and water drainage were also monitored. Meanwhile, we developed a dynamic dual-permeability preferential flow model by incorporating the shrinking–swelling model proposed by Stewart et al. (2016b). The performance of the model was evaluated by comparing the simulated results with measured data.

2 Experimental study

2.1 Testing apparatus

To investigate the effects of dynamic changes of desiccation cracks on preferential flow, a soil column infiltration test was conducted under wetting–drying cycles (abbreviated as WD cycles hereafter). The testing apparatus consisted of a rainfall-evaporation system, an environment monitoring device, a plexiglass column, a HD camera, hydrological sensors and a drainage measurement device (Fig. 1).

Figure 1 Schematic design and photos of the soil column test.

The rainfall-evaporation system included a rainfall simulator and two warm lamps as well as a small fan. The rainfall simulator was 0.5 m above the soil surface, which can produce rainfall with the intensity of 24–120 mm h⁻¹. The warm lamps and a small fan were put near the soil surface to accelerate water evaporation. The environment monitoring device consisted of a thermo-hygrometer that connected a probe above the soil surface to detect the environmental temperature and humidity and a water container to measure the potential evaporation.

The plexiglass column was composed of a column (with a height of 60 cm and a diameter of 50 cm) placed on a catchment hopper which was used to collect and drain out water from the soil column.

The HD camera (TTQ-J2, constant focal length: 35 mm) was fixed on the slope above the soil surface to take photos at regular intervals during the drying periods.

Hydrological sensors, including five soil moisture content–temperature sensors (Acclima, TDR-310s, with a measurement moisture content range of 0 %–100 % and an accuracy of ±2 % and a temperature range of $-\mathrm{40}\pm \mathrm{60}$^∘C and an accuracy of ±0.2^∘C) and five water potential sensors (Campbell, WP-257, with a measurement range of −200–0 kPa and an accuracy of ±0.5 kPa), were used to monitor the hydrological response during WD cycles. Five TDR-310s and five WP-257s were inserted into the soil column from the two opposite sides of the plexiglass column, respectively, with the same height spacing of 10 cm from top to bottom.

The drainage measurement device, including two electronic balances, was used to record the cumulative water drainage from the soil column.

2.2 Materials

The soil used in the test was taken from Zongyang County, Anhui, China. Table 1 shows the basic physical parameters and main mineral composition of the soil samples. The soil found in this study is classified as weak expansive soil. The saturated hydraulic conductivity was measured on reconstituted soil cores with a dry density of 1.55 g cm⁻³ (the same as the soil column). In addition, the shrinkage curve of the saturated soil core was also obtained using a similar method proposed in Wen et al. (2021). The difference is that we measured the vertical deformation in regular time intervals instead of continuous monitoring. Figure 2 shows the variation of soil porosity with the volumetric water content.

Table 1 Basic physical parameters of the soil sample.

G_s – specific gravity (–); ω_opt – optimal moisture content (g g⁻¹), which refers to the water content corresponding to the maximum dry density; ρ_d,max – the maximum dry density (g cm⁻³); L_l – liquid limit (%); P_l – plastic limit (%); δ_ef – free swelling ratio (%); C_illite, C_kaolinite, C_quartz and C_albite – content of illite, kaolinite, quartz and albite, respectively (%); and K_s – saturated hydraulic conductivity (m s⁻¹).

Figure 2 Shrinkage curve of the test soil.

To ensure the homogeneity of the soil column, soil samples were compacted in 10 layers, and each layer was 5 cm thick. Prior to filling soil into the plexiglass column, the soil samples with the total weight required for each layer were prepared according to the designed density (dry density of 1.55 g cm⁻³) and gravimetric water content (10 %). Then, the soil samples were compacted in the plexiglass column using a rubber hammer. The soil column was constructed within 1 d. After that, the soil column was allowed to stand for 3 d to obtain stable records of the hydrological sensors.

2.3 Data collection

In the soil column test, the following data were collected:

1. Boundary conditions were collected, including rainfall intensity (r; mm h⁻¹) and potential evaporation (PE; mm h⁻¹) at a 1 h time interval and temperature (T; ^∘C) and relative humidity (RH; %) at a 5 min time interval.

2. Hydrological data, including volume water content (θ_exp; %) and soil matrix suction (S_exp; kPa), as well as cumulative drainage from the top (D_top; g) and bottom (D_bottom; g) of the soil column, were collected at different depths at a 5 min time interval.

3. Crack geometric data, including the crack ratio (w_c,exp), the crack aperture (w_j,exp) and the maximum crack depth (d_max; mm), were collected. The w_c,exp and w_j,exp were obtained via processing of the crack photos, which were taken at 20 min intervals during drying periods. The image processing method mainly includes two steps as shown in Fig. 3. The d_max was measured by thin wire before each rainfall event.

Figure 3 Process of crack image processing. (a) A photo obtained from the HD camera, 800 pixels × 1400 pixels. (b) A crack image after cropping and pixel enhancement, 1044 pixels × 1005 pixels. (c) A crack image after binarization and denoising; the crack ratio was calculated as the crack area divided by the overall area of interest, and the crack aperture was calculated as the average value of crack aperture from three different positions.

2.4 Test procedure

The overall experimental process included two stages of WD cycles. The purpose of the first stage was to generate a relatively stable surface pattern of the desiccation cracks. The first stage was from 5 January 2022 at 15:00 to 28 February 2022 at 09:00 and included 13 WD cycles.

The second stage was from 28 February 2022 at 09:00 to 28 March 2022 at 22:30 and included seven WD cycles. Figure 4 presents the variation of rainfall, evaporation, temperature and relative humidity in the entire experiment process. Because the two warm lamps and fan were closed during the night, two kinds of evaporation intensity can be observed during the drying periods. In addition, the average environment temperature in the fifth WD cycle was higher because we turned up the power of the two warm lamps. In this current study, we mainly focus on the second stage of WD cycles.

Figure 4 Environmental conditions of the experiment. (a) Time series of temperature and relative humidity. (b) Rainfall intensity and potential evaporation.

3 Model description

3.1 Dual-permeability model (DPM)

The DPM concept used in this study corresponds to the one developed by Gerke and van Genuchten (1993a). The model divides the flow domain into two overlapping and interacting continua according to the volumetric ratios of each domain, where two coupled 2-D Richards' equations are used to describe the matrix flow and preferential flow as

$$\begin{array}{}\text{(1)}& {\displaystyle }{C}_{\mathrm{c}}\left(h\right){\displaystyle \frac{\partial h_{\mathrm{c}}}{\partial t}}=\mathrm{\nabla }\left[K_{\mathrm{c}}\right(h\left)\mathrm{\nabla }\right(h_{\mathrm{c}}+z\left)\right]-{\displaystyle \frac{{\mathrm{\Gamma }}_{\mathrm{w}}}{w_{\mathrm{c}}}}\text{(2)}& {\displaystyle }{C}_{\mathrm{m}}\left(h\right){\displaystyle \frac{\partial h_{\mathrm{m}}}{\partial t}}=\mathrm{\nabla }\left[K_{\mathrm{m}}\right(h\left)\mathrm{\nabla }\right(h_{\mathrm{m}}+z\left)\right]+{\displaystyle \frac{{\mathrm{\Gamma }}_{\mathrm{w}}}{w_{\mathrm{m}}}}\text{(3)}& {\displaystyle }{\mathrm{\Gamma }}_{\mathrm{w}}={\mathit{\alpha }}_{\mathrm{w}}{K}_{\mathrm{a}}(h_{\mathrm{c}}-h_{\mathrm{m}})\text{(4)}& {\displaystyle }{w}_{\mathrm{c}}+w_{\mathrm{m}}=\mathrm{1},\end{array}$$

where the subscripts “c” and “m” indicate the crack and matrix domains, respectively; h (m) is the pressure head; C represents the specific water capacity, $\mathrm{d}\mathit{\theta }/\mathrm{d}h$ (1 m⁻¹); θ (–) is the volumetric water content; K (m s⁻¹) is the isotropic hydraulic conductivity; z (m) is the elevation head; w (–) is the volumetric ratio of the crack domain or matrix domain over the bulk soil volume; Γ_w is the water exchange term (1 s⁻¹) between the two domains; α_w (1 m⁻²) is the effective water transfer coefficient; and K_a (m s⁻¹) is the interface hydraulic conductivity.

The hydraulic properties of the two domains are parameterized based on the Mualem–van Genuchten soil water retention curves (SWRCs) (Mualem, 1976; van Genuchten, 1980) as

$$\begin{array}{}\text{(5)}& {\displaystyle }{S}_{\mathrm{e}}\left(h\right)={\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}={\left[\mathrm{1}+(\left|\mathit{\alpha }h\right|{)}^n\right]}^{-m}\text{(6)}& {\displaystyle }K\left(S_{\mathrm{e}}\right)=K_{\mathrm{s}}{K}_{\mathrm{r}}\left(S_{\mathrm{e}}\right)=K_{\mathrm{s}}{S}_{\mathrm{e}}^{\mathrm{0.5}}[\mathrm{1}-(\mathrm{1}-S_{\mathrm{e}}^{\mathrm{1}/m}{)}^m{]}^{\mathrm{2}},\end{array}$$

where S_e (–) is the effective saturation; θ_s (–) and θ_r (–) are the saturated and residual volumetric water content, respectively; α (1 m⁻¹), n (–) and m (–) are fitting parameters; K_s (m s⁻¹) is the saturated hydraulic conductivity; and K_r (–) is the relative hydraulic conductivity.

According to Gerke and van Genuchten (1993a), the total porosity ε (–), total volume water content θ (–), total hydraulic conductivity K (m s⁻¹) and total volumetric flux (m s⁻¹) in terms of the volume ratio of each domain can be expressed as

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathit{\epsilon }=w_{\mathrm{c}}{\mathit{\epsilon }}_{\mathrm{c}}+w_{\mathrm{m}}{\mathit{\epsilon }}_{\mathrm{m}}\text{(8)}& {\displaystyle }\mathit{\theta }=w_{\mathrm{c}}{\mathit{\theta }}_{\mathrm{c}}+w_{\mathrm{m}}{\mathit{\theta }}_{\mathrm{m}}\text{(9)}& {\displaystyle }K=w_{\mathrm{c}}{K}_{\mathrm{c}}+w_{\mathrm{m}}{K}_{\mathrm{m}}.\end{array}$$

Note that the total porosity ε is defined as the total pore volume (V_p) divided by total soil volume (V), while ε_m (or ε_c) is defined as the pore volume in the matrix (V_p,m) (or crack, V_p,c) domain divided by the volume of that domain (V_m or V_c). The total volume water content has the same definition.

In the case of a DPM, a specified flux i is divided between the matrix and crack domains as

$$\begin{array}{}\text{(10)}& i=w_{\mathrm{c}}{i}_{\mathrm{c}}+w_{\mathrm{m}}{i}_{\mathrm{m}},\end{array}$$

where i_c and i_m are the effective boundary fluxes into each domain (m s⁻¹).

Considering a rainfall condition, the effective boundary fluxes of the two domains are initially equal to rainfall intensity (r) due to the infiltration capacity of each domain being larger than r (Dusek et al., 2008), and therefore the boundary fluxes of each domain can be written as

$$\begin{array}{}\text{(11)}& {\displaystyle }{i}_{\mathrm{c}}=r\text{(12)}& {\displaystyle }{i}_{\mathrm{m}}=r.\end{array}$$

As the soil keeps wetting, the decrease of the pressure head gradient may firstly lead to the infiltration capacity of the matrix domain dropping to a value less than r. Then, ponding occurs on the surface of the soil matrix, and the boundary condition changes to a specified pressure head boundary. This transformation can be achieved in COMSOL, a multiphysics solver and simulation software package building on the finite element method, using a combined type of boundary (Dirichlet and Neumann) proposed by Chui and Freyberg (2009). Once ponding occurs on the matrix domain, the surplus water from that domain infiltrates into the crack domain, and its effective flux increases to

$$\begin{array}{}\text{(13)}& i_{\mathrm{c}}=(r-w_{\mathrm{m}}{i}_{\mathrm{m}})/w_{\mathrm{c}}.\end{array}$$

When the retained water volume in the cracks exceeds its storage capacity, water will pond on the surface of the crack domain.

Considering an evaporation condition, the Wilson–Fredlund–Barbour–Penman experimental function model (Wilson et al., 1997) was used to calculate the actual evaporation of each domain:

$$\begin{array}{}\text{(14)}& \mathrm{AE}/\mathrm{PE}=\exp \left({\displaystyle \frac{-Sg{\mathit{\omega }}_v}{\mathit{\xi }(\mathrm{1}-h_{\mathrm{a}}){\mathit{\gamma }}_{\mathrm{w}}R(T_{\mathrm{s}}+\mathrm{273.15})}}\right),\end{array}$$

where AE is the actual evaporation; PE is the potential evaporation measured in the experiment; S (kPa) is total matric suction at the soil surface; g (m s⁻²) is the gravitational acceleration constant; ω_v is molecular mass of water, 0.018 kg mol⁻¹; ξ is a dimensional empirical parameter with a suggested value of 0.7; h_a is relative humidity of overlying air; γ_w is unit mass of water, 9.807 kN m⁻³; R is the universal gas constant, 8.314 J (mol K)⁻¹; and T_s (^∘C) is the soil surface temperature.

3.2 Dynamic dual-permeability model (DPMDy)

3.2.1 Porosity description

In Stewart et al. (2016a, b) and Stewart (2018), the total porosity (ϕ_max) of a cracked soil was divided into three domains: aggregates (or soil matrix), cracks (voids from horizontal deformation induced by desiccation cracks) and subsidence (voids from vertical deformation induced by desiccation cracks). In Stewart et al. (2016a), the distributions of these domains change as a function of a unified water content, U:

$$\begin{array}{}\text{(15)}& {\mathit{\varphi }}_{max}={\mathit{\varphi }}_{\mathrm{matrix}}\left(U\right)+{\mathit{\varphi }}_{\mathrm{crack}}\left(U\right)+{\mathit{\varphi }}_{\mathrm{sub}}\left(U\right),\end{array}$$

where the subscripts “matrix”, “crack” and “sub” refer to the aforementioned three domains. In this study, we assume that the horizontal deformation dominates the formation of desiccation cracks; thus ϕ_sub(U) can be neglected.

Stewart et al. (2016a) then deduced the porosities of each domain as

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathit{\varphi }}_{\mathrm{matrix}}\left(U\right)=({\mathit{\varphi }}_{max}-{\mathit{\varphi }}_{min})\left({\displaystyle \frac{p+\mathrm{1}}{p+U^{-q}}}\right)+{\mathit{\varphi }}_{min}\text{(17)}& {\displaystyle }{\mathit{\varphi }}_{\mathrm{crack}}\left(U\right)=({\mathit{\varphi }}_{max}-{\mathit{\varphi }}_{min})\left({\displaystyle \frac{\mathrm{1}-U^q}{\mathrm{1}+p{U}^q}}\right),\end{array}$$

where p and q are functional shape parameters; ϕ_max is the maximum porosity of a soil core prior to shrinkage and thus also represents the total porosity; ϕ_min is the minimum porosity of the matrix domain; and U is a unified water content (defined as water content u divided by its saturated value u_max), which can be approximately estimated to be the saturation degree (S_e,m) in an SWRC function of the soil matrix (Stewart et al., 2016a). Indeed, Eq. (16) represents a shrinkage curve function in which four parameters can be obtained through a shrinkage test.

Substituting S_e,m as U and incorporating Eq. (5) into Eqs. (16) and (17), we can obtain the porosity of the two domains as a function of pressure head h:

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{rl}{\mathit{\varphi }}_{\mathrm{matrix}}\left(h\right)& =({\mathit{\varphi }}_{max}-{\mathit{\varphi }}_{min})\left({\displaystyle \frac{p+\mathrm{1}}{p+S_{\mathrm{e},\mathrm{m}}^{-q}}}\right)+{\mathit{\varphi }}_{min}\\ & =({\mathit{\varphi }}_{max}-{\mathit{\varphi }}_{min})\left({\displaystyle \frac{p+\mathrm{1}}{p+{\left({\left[\mathrm{1}+(\left|{\mathit{\alpha }}_{\mathrm{m}}{h}_{\mathrm{m}}\right|{)}^{n_{\mathrm{m}}}\right]}^{-m_{\mathrm{m}}}\right)}^{-q}}}\right)\\ & +{\mathit{\varphi }}_{min}\end{array}\text{(19)}& {\displaystyle }\begin{array}{rl}{\mathit{\varphi }}_{\mathrm{crack}}\left(h\right)& =({\mathit{\varphi }}_{max}-{\mathit{\varphi }}_{min})\left({\displaystyle \frac{\mathrm{1}-S_{\mathrm{e},\mathrm{m}}^q}{\mathrm{1}+p{S}_{\mathrm{e},\mathrm{m}}^q}}\right)\\ & =({\mathit{\varphi }}_{max}-{\mathit{\varphi }}_{min})\left({\displaystyle \frac{\mathrm{1}-{\left({\left[\mathrm{1}+(\left|{\mathit{\alpha }}_{\mathrm{m}}{h}_{\mathrm{m}}\right|{)}^{n_{\mathrm{m}}}\right]}^{-m_{\mathrm{m}}}\right)}^q}{\mathrm{1}+p{\left({\left[\mathrm{1}+(\left|{\mathit{\alpha }}_{\mathrm{m}}{h}_{\mathrm{m}}\right|{)}^{n_{\mathrm{m}}}\right]}^{-m_{\mathrm{m}}}\right)}^q}}\right).\end{array}\end{array}$$

With these porosity equations in mind, we can rewrite Eqs. (4) and (7) as

$$\begin{array}{}\text{(20)}& {\mathit{\varphi }}_{max}=w_{\mathrm{c}}{\mathit{\epsilon }}_{\mathrm{c}}+(\mathrm{1}-w_{\mathrm{c}}){\mathit{\epsilon }}_{\mathrm{m}}.\end{array}$$

Because the crack domain is mainly composed of voids, we assume here that V_p,c is equal to V_c, and thus ε_c=1. Through this assumption, we obtained a physically consistent definition of how the porosity and crack volume vary as functions of saturation degree as follows:

$$\begin{array}{}\text{(21)}& {\displaystyle }{w}_{\mathrm{c}}{\mathit{\epsilon }}_{\mathrm{c}}=w_{\mathrm{c}}={\mathit{\varphi }}_{\mathrm{crack}}\left(S_{\mathrm{e},\mathrm{m}}\right)=({\mathit{\varphi }}_{max}-{\mathit{\varphi }}_{min})\left({\displaystyle \frac{\mathrm{1}-S_{\mathrm{e},\mathrm{m}}^q}{\mathrm{1}+p{S}_{\mathrm{e},\mathrm{m}}^q}}\right)\text{(22)}& {\displaystyle }\begin{array}{rl}{\mathit{\epsilon }}_{\mathrm{m}}& ={\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{matrix}}\left(S_{\mathrm{e},\mathrm{m}}\right)}{\mathrm{1}-w_{\mathrm{c}}}}=\left[({\mathit{\varphi }}_{max}-{\mathit{\varphi }}_{min})\left({\displaystyle \frac{p+\mathrm{1}}{p+S_{\mathrm{e},\mathrm{m}}^{-q}}}\right)+{\mathit{\varphi }}_{min}\right]\\ & /\left[\mathrm{1}-({\mathit{\varphi }}_{max}-{\mathit{\varphi }}_{min})\left({\displaystyle \frac{\mathrm{1}-S_{\mathrm{e},\mathrm{m}}^q}{\mathrm{1}+p{S}_{\mathrm{e},\mathrm{m}}^q}}\right)\right].\end{array}\end{array}$$

3.2.2 Water content and hydraulic conductivity

In terms of Eq. (8), the total water content of the soil volume can be expressed as

$$\begin{array}{}\text{(23)}& \mathit{\theta }={\mathit{\varphi }}_{\mathrm{crack}}\left(h\right){\mathit{\theta }}_{\mathrm{c}}+(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{crack}}(h\left)\right){\mathit{\theta }}_{\mathrm{m}}.\end{array}$$

Regarding the hydraulic conductivity of each domain, the classical DPM often assumed it to be equal to the product of a fixed K_s and the relative hydraulic conductivity K_r of the corresponding domain. In the DPMDy, the van Genuchten–Mualem equation is adopted in the matrix domain, while in the cracks, where capillarity has little effect on water flow, we set K_r=1. Hence, the following equations are obtained:

$$\begin{array}{}\text{(24)}& {\displaystyle }\left\{\begin{array}{l}{K}_{\mathrm{m}}=K_{\mathrm{m},\mathrm{s}}{K}_{\mathrm{r}}\left(S_{\mathrm{e},\mathrm{m}}\right)=K_{\mathrm{m},\mathrm{s}}{S}_{\mathrm{e},\mathrm{m}}^{\mathrm{0.5}}[\mathrm{1}-(\mathrm{1}-S_{\mathrm{e},\mathrm{m}}^{\mathrm{1}/m_{\mathrm{m}}}{)}^{m_{\mathrm{m}}}{]}^{\mathrm{2}}\\ K_{\mathrm{c}}=K_{\mathrm{c},\mathrm{s}}{K}_{\mathrm{r}}\left(S_{\mathrm{e},\mathrm{c}}\right)=K_{\mathrm{c},\mathrm{s}}{S}_{\mathrm{e},\mathrm{c}}^{\mathrm{0.5}}[\mathrm{1}-(\mathrm{1}-S_{\mathrm{e},\mathrm{c}}^{\mathrm{1}/m_{\mathrm{c}}}{)}^{m_{\mathrm{c}}}{]}^{\mathrm{2}}\end{array}\right.\text{(25)}& {\displaystyle }{K}_{\mathrm{c}}=K_{\mathrm{c},\mathrm{s}},\end{array}$$

where K_c,s and K_m,s refer to the saturated hydraulic conductivity in crack and matrix domains, respectively.

However, the K_c,s and K_m,s are transient variables that change with the crack geometries in the crack domain and porosity in the matrix domain, which should be taken into consideration in a shrinking–swelling soil. To solve this issue, Stewart et al. (2016b) further deduced models that describe the relationships between K_m,s, K_c,s and S_e,m.

$$\begin{array}{}\text{(26)}& {\displaystyle }{K}_{\mathrm{m},\mathrm{s}}\left(S_{\mathrm{e},\mathrm{m}}\right)=K_{\mathrm{m},\max }\left({\displaystyle \frac{p+\mathrm{1}}{p+S_{\mathrm{e},\mathrm{m}}^{-q}}}\right)\text{(27)}& {\displaystyle }{K}_{\mathrm{c},\mathrm{s}}\left(S_{\mathrm{e},\mathrm{m}}\right)=K_{\mathrm{c},max}{\left({\displaystyle \frac{\mathrm{1}-S_{\mathrm{e},\mathrm{m}}^q}{\mathrm{1}+p{S}_{\mathrm{e},\mathrm{m}}^q}}\right)}^{\mathrm{2}},\end{array}$$

where $K_{\mathrm{c},max}$ is the maximum saturated hydraulic conductivity of the crack domain (at $S_{\mathrm{e},\mathrm{m}}=$ 0) when the crack aperture achieves the maximum value, and $K_{\mathrm{m},max}$ is the maximum saturated hydraulic conductivity of the matrix domain (at $S_{\mathrm{e},\mathrm{m}}=$ 1) when the radius of cylindrical pores in that domain achieves the maximum value (see Eqs. 25 and 27 in Stewart et al., 2016b). The assumption $K_{\mathrm{r}}\left(S_{\mathrm{e},\mathrm{c}}\right)=$ 1 means that the magnitude of K_c only depends on the crack area or the saturated degree of the soil matrix domain.

$$\begin{array}{}\text{(28)}& K_{\mathrm{c}}=K_{\mathrm{c},\mathrm{s}}=K_{\mathrm{c},\mathrm{s}}\left(S_{\mathrm{e},\mathrm{m}}\right)=K_{\mathrm{c},max}{\left({\displaystyle \frac{\mathrm{1}-S_{\mathrm{e},\mathrm{m}}^q}{\mathrm{1}+p{S}_{\mathrm{e},\mathrm{m}}^q}}\right)}^{\mathrm{2}}\end{array}$$

Incorporating Eqs. (26) and (27) into Eq. (9) obtains

$$\begin{array}{}\text{(29)}& \begin{array}{rl}{K}_{\mathrm{s}}& ={\mathit{\varphi }}_{\mathrm{crack}}\left(h\right)K_{\mathrm{c},max}{\left({\displaystyle \frac{\mathrm{1}-S_{\mathrm{e},\mathrm{m}}^q}{\mathrm{1}+p{S}_{\mathrm{e},\mathrm{m}}^q}}\right)}^{\mathrm{2}}\\ & +(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{crack}}(h\left)\right)K_{\mathrm{m},max}\left({\displaystyle \frac{p+\mathrm{1}}{p+S_{\mathrm{e},\mathrm{m}}^{-q}}}\right).\end{array}\end{array}$$

Note that $K_{\mathrm{m},max}$ can be obtained by a laboratory-based infiltration test through a saturated soil core prior to shrinkage. Then, Eq. (29) can be used to fit the $K_{\mathrm{c},max}$ through the overall saturated hydraulic conductivity (measured K_s) under different crack volume or ratio. Alternatively, $K_{\mathrm{c},max}$ can also be approximately calculated as

$$\begin{array}{}\text{(30)}& K_{\mathrm{c},max}={\displaystyle \frac{w_{j,max}^{\mathrm{2}}g}{\mathrm{12}v}},\end{array}$$

where $w_{j,max}$ stands for the maximum crack aperture measured in the experiment (m), g is the gravity acceleration constant (m s⁻²) and v is the water kinematic viscosity (m² s⁻¹). This equation is a relation to the cubes of the aperture of a crack with respect to the crack inner flux, which is based on the derivation of laminar flow between parallel plates for Hagen–Poiseuille type of flow (Snow, 1965).

Eventually, we can simulate the hydrological process considering the dynamic changes of desiccation cracks by incorporating Eqs. (19), (21), (26), (27) and (28) into the DPM.

4 Experimental results

4.1 Crack dynamic changes

Figure 5 presents typical images of crack evolution during each WD cycle. Intuitively, it seems that the crack area and width did not show an obvious increasing trend with the WD cycles as expected. Conversely, during the first to fourth WD cycles, the cracks at the same moment after rainfall (Fig. 5b2–4) and the final state (Fig. 5c2–4) decreased significantly, even though the environmental temperature (T) and the potential evaporation (PE) increased in these periods. The cracks increased significantly from the fifth WD cycle on, but most of them were finer than before. Overall, cracks in the first WD cycle are wider than those formed in other cycles.

Figure 5 Typical images of crack evolution in seven wetting–drying cycles. (a1–7) Water ponds on the soil surface after rainfall, (b1–7) crack images at the 2135th min after each rainfall event and (c1–7) crack images at the end of the final high-evaporation period during each wetting–drying cycle.

Figure 6 quantitatively shows the variation of crack ratio (w_c,exp) and crack aperture (w_j,exp) in the experiment. Overall, the variation curves corresponded to the intuitive descriptions mentioned above. In particular, an unexpected result was that the T and PE in the fifth and sixth WD cycles were higher than in previous cycles, but their maximum w_c,exp and w_j,exp became smaller. During a single WD cycle, the w_c,exp and w_j,exp have a similar trend, which shows a dramatic decrease during rainfall, rapid increase in high-evaporation periods, and slow increase or even decrease in low-evaporation periods. More specifically, during the rainfall periods, the crack closure process was not significant until the water ponded on the soil matrix, then ponded water flowed into the cracks, leading to acceleration of the crack closure. Note that cracks were not completely closed, even when they were full of water (Fig. 5a1–7). The minimum crack ratio under such conditions is approximately 0.1 %. In the evaporation periods, the maximum crack ratio reaches 2.87 %, and the maximum crack aperture reaches 2.6 mm. In addition, Fig. 7 shows the maximum crack depth (d_max) measured after each cycle. It can be seen that d_max increased substantially after the first WD cycle and then slightly increased in the last six cycles, with a maximum value of 23.8 cm.

Figure 6 Time series of crack geometries. (a) Crack ratio. (b) Crack aperture.

Figure 7 The maximum crack depth measured after each wetting–drying cycle.

4.2 Hydrological response

Table 2 presents the manually recorded results of external hydrological responses involving ponding and drainage during each WD cycle. It can be seen that the ponding occurred on the soil surface within 5 min after each rainfall event. The ponding duration in each rainfall event mainly decreased with WD cycles. Note that the ponding depth in each rainfall event was below the upper drainage outlet. Regarding the water drainage, approximately 1.4 kg of water (the total water mass was 8 kg) was leaked during the first rainfall event due to the interspace between the soil and the plexiglass column and the hydrological sensors. Then, we sealed the interspace using clay powder and polyurethane cement (soft materials without constraining effects on the soil swelling) after each drying process, and subsequently, no water drainage was observed at the bottom outlet.

Table 2 Manual readings of external hydrological responses.

t_p (min) – beginning of ponding after each rainfall event.

Figure 8 shows the internal hydrological responses recorded by the soil moisture and water potential sensors. Because the M2 and M4 were damaged during soil compaction, no matric suction data were obtained at their depths. Overall, water content at all depths increased during rainfall and decreased during evaporation, where T1 showed the most sensitive responses to the WD cycles. During rainfall, the time for water content to respond to each rainfall event increased with depth, but the time difference among all depths decreased significantly from the second WD cycle on. During the drying periods, an interesting phenomenon was that the water content at 5 cm depth showed an overall decline trend, but transient increases of water content frequently appeared during low-evaporation periods. Such transient increases seem to be related to the slow decrease of crack ratio as mentioned in Sect. 4.1. Regarding the matric suction, its variation trend was similar to the water content but showed more delayed responses to the environmental conditions, especially in the last three WD cycles. Additionally, Fig. 8b also implies that soil at 5 cm depth reached saturation during each rainfall event, while soil below the 25 cm depth was in the unsaturated state in the whole experiment process.

Figure 8 Time series of volume water content (a) and matric suction (b) at different depths.

5 Numerical simulation

5.1 Setup of numerical model

The single-domain model (SDM), dual-permeability model (DPM) and dynamic DPM (DPMDy) were implemented in a finite element solver for Richards' equation as part of the COMSOL Multiphysics^® software (COMSOL 5.6). As shown in Fig. 9, they have the same 2-D size, boundary conditions, mesh structure and initial condition. The model domain is 0.5 m by 0.5 m, the same as the soil column. Because the measured maximum crack depth was 23.8 cm, we specified the crack domain existing within the upper 25 cm depth of the soil column.

Figure 9 Setup of the 2-D numerical model for the SDM, DPM and DPMDy.

The boundary conditions at the top were set as a combined type of boundary conditions (as mentioned in Sect. 3.1) for representing the rainfall, ponding and evaporation process recorded in the experiment; the bottom side is a seepage boundary condition; and the left and right sides of the model are no-flux boundaries.

Because the pressure head in the surface area may change frequently and drastically during WD cycles, a refined mesh structure with dense boundary layers was used to capture the transient hydrological conditions. The boundary layers included 15 layers of rectangular grid, with a minimum and maximum thickness of approximately 0.04 and 0.3 cm, respectively. A coarser free-triangle mesh (average length of 1.8 cm) was defined below the boundary layers. The initial condition both in matrix and crack domains was set as the distribution of porewater pressure measured from the experiment prior to the first WD cycle.

5.2 Parameters

5.2.1 Shrinkage parameters

As shown in Fig. 10, using Eq. (18) to fit the measured shrinkage curve in Fig. 2, we obtained the four shrinkage parameters as ${\mathit{\varphi }}_{min}=\mathrm{0.22}$, ${\mathit{\varphi }}_{max}=\mathrm{0.30}$, $p=\mathrm{8.8}\pm \mathrm{4.84}$ and $q=\mathrm{2.71}\pm \mathrm{0.85}$. Then, the variation of porosity in the crack domain (or crack ratio w_c) and the matrix domain (ε_m) could be obtained using Eqs. (21) and (22), respectively. Note that the minimum w_c calculated by Eq. (21) was set as 0.001 considering the incomplete closure of cracks during rainfall.

Figure 10 Fitted shrinkage curve (solid line) and modeled porosity variation of matrix (dashed line) and crack domains (dashed–dotted line).

5.2.2 Soil water retention parameters

Figure 11 shows the measured matric suction versus volume water content at different depths. The SWRCs were estimated using a best fit of the van Genuchten–Mualem equation to measured soil water retention data. It can be seen that the WD cycles lead to hysteretic curves in the SWRC at 5 cm and 25 cm depths, while those at the 45 cm depth rarely show hysteretic curves. This result may also indicate that most of the cracks exist within the upper 25 cm depth of the soil column. In this study, we simply estimated an approximate single SWRC of the soil matrix through experiment data instead of incorporating the hysteretic curves into the model. For instance, the estimated SWRC in Fig. 11a lies between the wetting SWRC and drying SWRC to capture the overall characteristics of wetting–drying SWRCs as far as possible. Note that the shape parameter n in the upper matrix domain is slightly smaller than the lower one considering the upper soil matrix may become denser after long-term WD cycles (13 times, 54 d). Regarding the SWRC of the crack domain, as we assume the crack domain does not contain any solids, the saturated water content (θ_c,s) and the residual water content (θ_c,s) of that domain were set to be 0.99 and 0.01, respectively. Meanwhile, because the SWRC of the crack domain cannot be experimentally determined, we assigned the other two SWRC parameters (a= 1.5 and n_c= 2) to mimic coarse-textured soil like behavior and to be consistent with Poiseuille law, which implies little effect of capillarity on flow through the cracks.

Figure 11 Measured and estimated SWRC at different depths: (a) 5 cm, (b) 25 cm and (c) 45 cm.

5.2.3 Hydraulic conductivity

As mentioned in Eq. (29), the maximum saturated hydraulic conductivity of matrix domain ($K_{\mathrm{m},max}$) equals the saturated hydraulic conductivity (K_s) measured in the laboratory. Here, we set $K_{\mathrm{m},max}=\mathrm{1.16}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹. Regarding the$K_{\mathrm{c},max}$, it was calculated using Eq. (30), where the $w_{j,max}$ was set to 2.6 mm obtained from Fig. 6b. Then, the variation curve of transient saturated hydraulic conductivity of the matrix domain (K_m,s) and the crack domain (K_c,s) could be obtained using Eqs. (27) and (28), respectively. Note that here we slightly modified Eq. (28) as follows:

$$\begin{array}{}\text{(28-b)}& K_{\mathrm{c},\mathrm{s}}\left(S_{\mathrm{e},\mathrm{m}}\right)=K_{\mathrm{c},max}{\left({\displaystyle \frac{\mathrm{1}-S_{\mathrm{e},\mathrm{m}}^q}{\mathrm{1}+p{S}_{\mathrm{e},\mathrm{m}}^q}}\right)}^{\mathrm{2}}+K_{\mathrm{c},min}.\end{array}$$

This modification not only avoided the K_c,s dropping to zero, thus benefitting the numerical convergence, but also was reasonable when considering the incomplete closure of cracks during rainfall. The $K_{\mathrm{c},min}$ was also estimated using Eq. (30) with a suggested $w_{j,max}=$ 0.01 mm. Further, the variation of K_m and K_c with the pressure head (h) in the DPMDy could be calculated by combining Eqs. (24), (26) and (28). Figure 12 presents K_m and K_c in the three models. Note that the pressure head in K_c(h_m) of the DPMDy refers to that of the matrix domain (h_m), while h in K_c(h_c) of the DPM refers to that of the crack domain (h_c).

Figure 12 Modeled hydraulic conductivity of each domain in the three models. (a) Saturated hydraulic conductivity of the matrix domain, (b) saturated hydraulic conductivity of the crack domain, (c) transient hydraulic conductivity of the matrix domain and (d) transient hydraulic conductivity of the crack domain.

5.2.4 Water exchange between and pore domains

In the dual-permeability concept, another important parameter is the hydraulic conductivity of the interface between matrix and crack domains (K_a). Generally, K_a was often estimated as the arithmetic mean of hydraulic conductivity of the two domains (Arora et al., 2011; Coppola et al., 2012, 2015; Gerke and van Genuchten, 1993b; Laine-Kaulio et al., 2014; Shao et al., 2015). However, this approximation may overestimate the K_a when the hydraulic conductivity of the crack domain is much higher than that of the matrix domain, especially in cracked clays. In our current study, a K_a function reformulated by Gerke et al. (2013) was adopted.

$$\begin{array}{}\text{(31)}& K_{{\mathrm{a}}_{min}}=\left\{\begin{array}{ll}\min \left\{K_{\mathrm{m}}\left(h_{\mathrm{c}}\right),K_{\mathrm{c}}\left(h_{\mathrm{c}}\right)\right\}& h_{\mathrm{c}}\ge h_{\mathrm{m}}\\ & \\ \min \left\{K_{\mathrm{m}}\left(h_{\mathrm{m}}\right),K_{\mathrm{c}}\left(h_{\mathrm{m}}\right)\right\}& h_{\mathrm{c}}<h_{\mathrm{m}}\end{array}\right.\end{array}$$

This formulation represents the flow occurring from the highest head toward the lowest head but being regulated by the less permeable of the two subsystems in that instant of time (Aguilar-López et al., 2020).

Regarding the α_w, experimental results presented by Song et al. (2018) showed that the saturated K_a may be 1 order of magnitude larger than the K_m,s, which will represent an enlarging coefficient ranging from 10 to 18. Hence, the α_w was set as 10 m⁻² considering the saturated $K_{{\mathrm{a}}_{min}}$ determined by Eq. (31) is equal to the K_m,s.

All parameters and parametric methods for the SDM, DPM and DPMDy are listed in Table 3.

Table 3 Summary of parameters and parametric methods for the SDM, DPM and DPMDy.

SDM: single-domain model. DPM: dual-permeability model neglecting crack dynamic changes. DPMDy: dynamic DPM. w_c= 0.01 and 0.03 refer to the average and the maximum value of the measured crack ratio, respectively.

5.3 Simulation results

5.3.1 Boundary flow

Figure 13 shows the temporal evolution of the boundary flow velocity simulated by the SDM, DPM and DPMDy. As shown in Fig. 13a1 and a3, during drying periods, the matrix domain dominates the soil evaporation process and was responsible for 97 %–99 % of the total evaporation in all the dual-permeability models. The matrix evaporation rate (e_m) simulated by the DPMDy was overall lower than that of the SDM and DPM during high-intensity evaporation periods, but the crack evaporation rate (e_c) simulated by the DPMDy, especially during the last three drying periods, was approximately 1 to 2 orders of magnitude larger than that of the DPM (see the enlarged image in Fig. A1).

Figure 13 Boundary flow simulated by the SDM, DPM and DPMDy. (a) Flow velocity of the boundary conditions and simulated results; (a1) and (a2) are the enlarged images of the flow velocity during the second drying and third wetting process, respectively; (a3) and (a4) are the enlarged images of the flow velocity during the fifth drying and sixth wetting process, respectively. The positive value is for infiltration and negative for evaporation.

With regard to the wetting process, Fig. 13a2 and a4 represent two typical infiltration patterns before and after the fifth drying period (with significantly increased evaporation intensity). Overall, matrix flow still dominated the infiltration process in all the dual-permeability models due to the relatively small crack ratio and depth. For the SDM, all the rainfall infiltrates into the soil during the beginning of rainfall events. When the soil surface gets saturated, water ponding occurred, and the soil infiltration rate gradually decreased. In the DPM and DPMDy, the surplus water after matrix ponding infiltrates into the crack domain as preferential flow, and water will pond on the overall soil surface when the crack domain reached its storage capacity. Recall that the crack volume in the DPMDy decreases with the matrix getting moist, while that in the DPM keeps constant. Consequently, the ponding time of the crack domain simulated by the DPMDy in the third rainfall event (inflection point of the dashed red line in Fig. 13a2) was 1.6 and 4.8 min earlier than that of the DPM-0.01 and DPM-0.03 (rigid dual-permeability model with crack ratio of 0.01 and 0.03, respectively), respectively. The cumulative preferential flow simulated by the DPMDy was 87.4 % and 95.2 % less than that of the DPM-0.01 and DPM-0.03, respectively. A similar rainfall pattern was obtained during the sixth rainfall event.

5.3.2 Water balance

By integrating the boundary flow velocity in Fig. 13a, the total cumulative flux for the experiment and the three models was obtained (Fig. 14a). In the experiment, the variation of water flux was estimated by calculating the sum of the difference between θ_ini (initial volume water content) and θ_t=i (volume water content at any time) in the five monitoring depths. Meanwhile, the water evaporation during water ponding was also estimated and added to the total flux volume. Regarding the numerical model, the water balance was obtained by integrating all flow components along the upper and lower boundaries. The steep increase stage of each curve represents cumulative input water flux during wetting periods, and the gradual decrease stage represents cumulative output water flux during drying periods. To evaluate the performance of each model on the water balance, the measured cumulative input and output water fluxes in each wetting and drying stage were compared to the simulated ones (Fig. 14b).

Figure 14 Water balance for the measured and simulated results. (a) Temporal evolution of total water flux calculated from the measured water content, SDM, DPM and DPMDy, and (b) measured versus modeled cumulative flux during each drying and wetting stage.

In Fig. 14a, the results show that the total infiltration (I_t,inf) and evaporation flux (E_t,eva) estimated from measured θ_exp were 171 and 138.95 mm, respectively. The I_t,inf was 5.86 % less than the supplied water (183.44 mm) due to the water leakage. The E_t,eva was 16.48 % less than the cumulative PE (166.36 mm) because of the limit of the soil actual evaporation. Regarding the simulation results, the coefficient of determination (R²) and intercept were used to evaluate the errors made by the three models. As shown in Fig. 14b, the slope of each fitting curve was fixed as 1. The SDM and DPMDy have relatively smaller intercepts and slightly higher R² than that of the DPM-0.01 and DPM-0.03, indicative of a better coincidence to the measured data. Overall, the errors in water balance caused by the three models were acceptable in this study.

5.3.3 Crack dynamic changes and hydrological response

Figure 15 shows part of the comparison results between the measured data and the three models. Detailed descriptions of all the comparison results are presented in Appendix A. Overall, all models show similar response trends with the measured data. Divergences among the three models mainly appeared during drying.

Figure 15 Temporal evolution of the measured and simulated crack ratio, matric suction and volumetric water content. (a) Measured and simulated crack ratio on soil surface, (b) measured and simulated matric suction at 25 cm depth and (c) measured and simulated total water content at 25 cm depth.

In Fig. 15a, the simulated surficial w_c,sim was not only generally close to the w_c,exp in value and trend, but also it captured the transient slow decrease in w_c,exp during low-evaporation periods. Notably, significant overprediction appeared in the sixth and seventh wetting–drying cycles.

In Fig. 15b, the matric suction (S_sim) at the 25 cm depth simulated by SDM and DPMDy was close to each other and had an average divergence 2.26 kPa to the measured data. The S_sim simulated by DPM had a greater average divergence of 3.4 kPa to the measured data. These values showed systematic underprediction compared to the S_sim simulated by SDM and DPMDy, but their differences became smaller with the increasing WD cycles.

In Fig. 15c, the total volumetric water content θ_sim simulated by SDM was much lower with respect to the DPMDy and DPM. The θ_sim simulated by DPM-0.01 and DPM-0.03 overpredicted the volumetric water content. The DPMDy provided better prediction results but also showed slight underprediction to the measured data in the last two WD cycles.

6 Discussions

6.1 Crack dynamic changes

Our experimental results demonstrated that the crack evolution is not always positively correlated to the increase of the WD cycles, T and PE. For instance, the 5 cm θ_exp at the end of the final three WD cycles was lower than that in the first WD cycle due to the increased T and PE, but the maximum w_c,exp measured during the final three WD cycles was much less than that in the first WD cycle. From the energy-driven perspective, soil cracking and propagation can be regarded as a process that the shrinkage energy (or stress), built up from the evaporation and thermal radiation, was released until a critical moment when the tensile strength of soil is reached (Peron et al., 2009; Xu et al., 2022; Tian et al., 2022). If the environmental condition changes in a stable range, the desiccation cracks will vary within the crack pattern and the maximum w_j,exp that were formed under the maximum shrinkage energy. In this case, new desiccation cracks will not appear in the remained soil matrix during WD cycles (Fig. 5b1–b4). One reason is that the shrinkage energy can be fully released via previous cracks. The other reason is that the shrinkage energy is not high enough to split the soil matrix that has a denser structure (or higher tensile strength) than its initial state prior to shrinkage (Luo et al., 2021). However, once the evaporation rate and thermal radiation increase to exceed the stable range, higher shrinkage energy will lead to new cracks appearing in the soil matrix that will concurrently restrain the width increase of the previous cracks (Wang et al., 2018; Xu et al., 2021). This is the reason that cracks in the final three WD cycles are finer than the first four WD cycles. Our model describes the crack evolution mainly from the hydrological-driven perspective that assumed the surface crack pattern has become stable after undergoing 13 WD cycles and has a constant function relationship with the water content. Indeed, this assumption is reasonable for natural soils under atmospheric environmental conditions. However, our experiment not only used reconstituted soil, but also intensely changed the environmental conditions since the fifth WD cycle. Therefore, the model overpredicted w_c,exp at the end of the sixth and seventh WD cycles.

Figure 16 Crack images at t= 3702 and 4327 min (photo at the top of the panels) as well as the vertical distribution of water content in the numerical model (lower part of the panels) during the low-evaporation process. Panel (a) shows the beginning of the final low-evaporation stage during the first drying periods, while panel (b) shows the end of the final low-evaporation stage during the first drying periods.

In addition, another interesting phenomenon is the transient decrease in w_c,exp and increase in θ_exp measured at 5 cm depth during low-evaporation periods, which we called a “self-closure” process. In light of Figs. 6 and 8, the self-closure process appeared to be always accompanied by relatively high RH. From the insight of the experiment, it is natural and common to infer that the moist air wetted the surface soil from top to bottom, resulting in the self-closure phenomenon. Interestingly, our model does not incorporate the vapor flow into the boundary conditions, and also the evaporation boundary only involves the outflow of water, but it still managed to capture the self-closure process. Figure 16 shows the crack images at t= 3702 and 4327 min as well as the corresponding cloud chart of θ_sim. It can be seen that the soil surface became moist during the low-evaporation period, which is a typical external phenomenon reflecting the self-closure process. The simulation results show that θ_sim near the surface soil increased during evaporation, while θ_sim at deep soils decreased, indicative of evaporation causing the deep water to move up and wet the surface soil from bottom to top. We further found that the process occurred because the water flow driven by the soil water potential gradients, existing between the wet and dry soil layers, overcame the gravity. Indeed, this kind of “hydraulic lift” process frequently occurs in planted soils where root zone soil can force water flow from moist, deep soil layers to dry, shallow soil layers (Richards and Caldwell, 1987; Bauerle et al., 2008) but was rarely reported in homogeneous bare soil. We infer that the evaporation boundary conditions using Eq. (14) might play a positive role in causing water to move up and constraining it within the surficial soil depths when the evaporation intensity decreased. In any case, our results provide an additional possible explanation to the self-closure phenomenon. Further quantitative analysis based on a gas–liquid two-phase flow model is needed to compare the contribution of hydraulic lift and moist air to the self-closure process of cracks.

Figure 17 Variation of boundary K_m and K_c in each model during the fifth drying periods. (a) K_m and (b) K_c.

6.2 Water flow with dynamic changes of desiccation cracks

6.2.1 Water fluxes

As mentioned in Sect. 5.3.1, during the drying process, the matrix and crack evaporation simulated by the DPMDy are overall lower and higher than other models, respectively. It can be explained by looking at the variation of boundary K_m and K_c in each model. Take the time span in Fig. 13a2 as an example; because the DPMDy considers the effects of matrix shrinkage on the K_m using Eq. (26), the K_m,DPMDy is always approximately 20 % and 30 % lower than that of the SDM and DPM, respectively (Fig. 17a). Conversely, because the DPM links the K_c with the saturation degree of the crack domain (see Eq. 25), the K_c,DPM is destined to decrease with the decreased saturation degree of the crack domain induced by drying, while the K_c,DPMDy increases with the crack development induced by drying in light of Eq. (28-b). The ultimate K_c,DPMDy is 80 % higher than the K_c,DPM (Fig. 17b). Indeed, the decrease in K_c with the drying process is an unrealistic and physically unreasonable result. We can imagine that after long-term drought, the K_c,DPM will decline to nearly zero according to Fig. 12d, which will greatly underestimate the propagation of the PF-DC in the subsequent storm event. However, many laboratory and field experiments have observed that heavy rainfall following a long-term drought facilitated PF-DC (Baram et al., 2012, 2013; Greve et al., 2010; Kurtzman and Scanlon, 2011; Schlögl et al., 2022). By contrast, the DPMDy has the potential to capture this process for its increasing K_c with the enlarging desiccation crack during the long-term drought. In this study, because the experiment scale (or crack volume) is small, the increment of PF-DC simulated by the DPMDy after high-intensity evaporation is not significant (despite increment = 25 %), but we believe the DPMDy will have a better performance in a larger scale (i.e slope scale).

6.2.2 Water exchange and distribution

For the dual-permeability model, the two domains are coupled by the water exchange term (Eq. 3) that is governed by the pressure head difference between the two domains ($\mathrm{\Delta }h=h_{\mathrm{c}}-h_{\mathrm{m}}$), water exchange coefficient (α_w) and the hydraulic conductivity between the two domains (K_a). The higher the Γ_w, the quicker the two domains equilibrate. Generally, the higher Γ_w leads to faster water exchange from the crack domain into the matrix domain and thus boosts the contribution of preferential flow to the water distribution in the soil matrix. According to previous studies, the commonly used magnitude of the product of saturated α_wK_a in clay soils ranges from 10⁻⁵ m⁻¹ s⁻¹ (Aguilar-López et al., 2020) to 10⁻⁶ m⁻¹s⁻¹ (Coppola et al., 2012, 2015; Gerke and Maximilian Köhne, 2004; Vogel et al., 2000). In this study, the saturated α_wK_a is $\mathrm{1.16}\times {\mathrm{10}}^{-\mathrm{5}}$ m⁻¹ s⁻¹, which falls in the reasonable range. Building on the above statement, the Δh and water exchange rates (${\mathrm{\Gamma }}_{\mathrm{w}}/w_{\mathrm{m}}$) for both the DPM and DPMDy at the 5, 15 and 25 cm depths during the sixth rainfall event are graphed in Fig. 18.

Figure 18 Pressure head difference (a1–a3) and water exchange rate (b1–b3) between the two domains at the 5, 15 and 25 cm depths during the sixth rainfall event. The positive value of water exchange rate is for the water flowing from the crack to the matrix domain, while the negative value is for the opposite direction.

As shown in Fig. 18a1–a3, Δh at all depths simulated by both the DPM and DPMDy rapidly reaches a positive peak value and gradually decreases with the rainfall process. The rapidly increasing positive value is because the crack domain reaches saturation earlier than the surrounding soil matrix due to the influx of preferential flow and the small crack storage space in this study. The decrease in the Δh is ascribed to the increase in h_m with water exchanging from the crack to the matrix domain. Notably, the crack closure process during the rainfall process leads to a decrease in crack volume (or crack water storage space); the “water table” (saturated zone) in the shrinking cracks elevates faster than that in the constant larger crack volume, which means the h_c simulated by DPMDy is higher than the DPM-0.01 and DPM-0.03. Consequently, the time for Δh reaching the peak value simulated by the DPMDy is the earliest at all the three depths, followed by the DPM-0.01 and DPM-0.03. The ${\mathrm{\Gamma }}_{\mathrm{w}}/w_{\mathrm{m}}$ simulated by the DPMDy shows the similar trend to the Δh (Fig. 18b1–b3). During the sixth rainfall event, its cumulative ${\mathrm{\Gamma }}_{\mathrm{w}}/w_{\mathrm{m}}$ at the 5, 15 and 25 cm depths is (26 %, 50 %), (10 %, 26 %) and (3 %, 14 %) larger than that of the DPM-0.01 and DPM-0.03, respectively.

This result means that the crack closure during wetting benefits the building-up process of the pressure head in the crack domain and thus can promote water exchange from the crack into the matrix domain. It corresponds to some experimental results that the PF-DC also exists and leads water to rapidly infiltrate into soils even if desiccation cracks are nearly closed (Baram et al., 2012; Greve et al., 2010; Luo et al., 2021; Sander and Gerke, 2007; Tuong et al., 1996). It also means using the DPM may overestimate the flux of PF-DC but underestimate the water exchange coming from the PF-DC. Because the experimental scale, crack ratio and depth in this study are small, the difference of simulation result involving the matric suction and water content between the DPM and DPMDy is not very significant. However, we can imagine that the deviation caused by the DPM at a larger scale will be more significant, especially in a typical shrinking–swelling soil slope under long-term WD cycles.

6.3 Model performance

We evaluated the prediction errors of different models to the measured matric suction, water content and crack ratio using a fixed slope line as the same in Sect. 5.3.2 (see Fig. A3 and Table 4). Overall, the DPMDy, which incorporates the dynamic changes of desiccation cracks and hydraulic conductivity into the dual-permeability model, has an overall better performance than the SDM and DPM, as indicated by the small intercept and high R². With regard to the water flux, while the three models all give a good fit with the measured data, the DPM overpredicted the water flux of PF-DC but underestimated the water exchange from cracks to the soil matrix compared to other models. It implies that adopting a constant crack volume in the DPM, whether it is an average or a maximum value of the measured crack ratio, will overestimate the PF-DC, which may be unsuitable to evaluate the irrigation efficiency. With regard to the matric suction (or porewater pressure), although the SDM has good performance as the DPMDy does, it significantly underpredicted the volume water content and thus may overestimate landslide stability in a moisture-content-dependent threshold method. Further, we expect that the SDM may show much poorer performance if one applies it to scenarios where the cracks are deeper and the soil has a higher swelling–shrinking ability than that of our experiment. A comprehensive model sensitivity analysis will be conducted in our future work.

Table 4 Summary of fitting performance of different models to measured data.

Compared to other dynamic preferential flow models, the DPMDy developed in this study also has its unique advantages. Firstly, the variation of crack volume (or crack ratio) in our model is deduced from the changes of matrix porosity due to shrinkage and thus has a universal definition. Instead, Coppola et al. (2012, 2015) linked the crack ratio to the suction head with an empirical natural logarithm function, which is not transferable to other types of soils. Secondly, the results support the suitability, in the crack domain, where capillarity has little effect, of V–G SWRC with effective parameters and a constant relative hydraulic conductivity (K_r= 1). In fact, a common defect in classical DPMs is that they often set the hydraulic conductivity of the crack domain (K_c) to vary as a function of the saturated degree calculated from the SWRC of the crack domain (i.e Eq. 25). This will lead to an unreasonable extremely low K_c in drying initial conditions (Aguilar-López et al., 2020). Setting K_r= 1 ensures that the magnitude of K_c only depends on the crack area or the saturated degree of the soil matrix domain, which provides a potential solution for remedying the shortcoming mentioned above. Thirdly, compared to some dynamic preferential flow models neglecting the water exchange between the two domains (Jamalinia et al., 2020; Kroes et al., 2000; Luo et al., 2021; Stewart, 2018) or ones adopting an arithmetic mean of hydraulic conductivity of the two domains (Coppola et al., 2012, 2015; Laine-Kaulio et al., 2014; Shao et al., 2018) that tends to overestimate the water exchange, our model tentatively adopts an improved exchange term proposed by Gerke et al. (2013), which we showed to be a logical and satisfactory improvement in simulating water exchange in our experiment.

However, in the current study, the hysteresis effect was neglected in both the soil deformation and SWRC because we assumed the soil shrinking–swelling behavior has less influence on the pore-size distribution (or SWRC shape) but more influence on the porosity (or hydraulic conductivity). This assumption inevitably caused some errors when compared to the measured water content, especially for the surficial soil layer that has been significantly affected by the WD cycles. Our future work will try to incorporate the hysteresis effect into the current model to further improve the prediction strength. In addition, we have to remind the reader again that because the shrinking–swelling model in our method is developed based on the hydrological-driven perspective, it may be more suitable in the natural soil layer where the crack pattern already has a stable state after long-term WD cycles.

7 Conclusions

This study combined an experimental study and a numerical simulation to quantify the preferential flow induced by dynamic changes of desiccation cracks (PF-DC). A soil column infiltration test under wetting–drying conditions was conducted to investigate dynamic changes of desiccation cracks and the accompanying water infiltration process. The variation of crack geometry, including crack ratio, width and depth, was measured. The soil volumetric water content, matric suction and water drainage were also monitored. A new dynamic dual-permeability model (DPMDy) was developed to account for the PF-DC, which includes physically consistent functions in describing the variation of both porosity and hydraulic conductivity in crack and matrix domains. The performance of the single-domain model (SDM), rigid dual-permeability model (DPM) and DPMDy was evaluated by comparing their simulation results to the monitoring data.

Overall, the DPMDy performed not only better in its prediction of the crack evolution and hydrological response with respect to the SDM and DPM, but it also provided much better descriptions of the underlying physics involving the PF-DC. During the drying periods, the matrix evaporation modeled by the DPMDy is lower than that of the SDM and DPM due to considering the permeability decay induced by soil shrinkage. But the crack evaporation modeled in the DPMDy approach is the highest because it managed to capture the raised crack permeability induced by drying–enlarging desiccation cracks. Compared to the DPM with a fixed crack volume, the DPMDy revealed that the crack closure process during wetting will lead to a faster pressure head building-up process in the crack domain and higher water exchange rates from the crack to the matrix domain. Additionally, using a fixed crack ratio in the DPM, whether it is the maximum or the average value from the experiment data, will overestimate the infiltration fluxes of PF-DC but underestimate its contribution to the matrix domain.

The DPMDy developed here has a physically consistent definition. It remedies the shortcomings DPM and other dynamic preferential flow models in defining the dynamic changes of desiccation cracks and hydraulic properties of the crack domain and interface. Future works should focus on considering the hysteresis effect of the SWRC during wetting–drying cycles in the model and its application to complex field situations.

Appendix A

Figures A1 and A2 show the temporal evolution of the measured and simulated crack ratio on the soil surface, matric suction (negative porewater pressure) and volumetric water contents at the five monitoring depths (5, 15, 25, 35 and 45 cm).

In Fig. A1a, the simulated w_c,sim was not only generally close to the w_c,exp in value and trend, but also it captured the transient slow decrease in w_c,exp during low-evaporation periods.

In Fig. A1b–f, the matric suction (S_sim) simulated by SDM and DPMDy is close to each other and has average divergence of 2.75, 2.26 and 5.02 kPa to the measured data at the 5, 25 and 45 cm depths, respectively. The S_sim simulated by DPM has a greater average divergence of 2.78, 3.4 and 7.43 kPa to the measured data at the three corresponding depths.

In Fig. A2a–e, the volumetric water content θ_sim simulated by SDM was much lower than that simulated by DPMDy and DPM. In most depths (except the 5 and 45 cm depth), SDM systematically underpredicted the volumetric water content during both wetting and drying periods. By contrast, the θ_sim simulated by DPM-0.01 and DPM-0.03 overpredicted the volumetric water content. The DPMDy gave overall better prediction results in most depths but has significant divergences to the measured data at the depth of 5 cm, and so are the other two models.

Figure A1 Temporal evolution of the measured and simulated crack ratio and matric suction at different depths. (a) Measured and simulated crack ratio (dynamic DPM) on soil surface and (b–f) measured and simulated matric suction (single-domain model, DPM and dynamic DPM) at depths of 5, 15, 25, 35 and 45 cm.

Figure A2 Temporal evolution of the measured and simulated volumetric water content at depths of 5, 15, 25, 35 and 45 cm. Note that the simulated volumetric water content demonstrated here is the total volumetric water content that was combined with the combined matrix and crack domains using Eq. (8).

Figure A3 Scatter plots of modeled vs. measured data. Panels (a), (b), (c) and (d) refer to the matric suction simulated by SDM, DPMDy, DPM-0.03 and DPM-0.01, respectively; panels (e), (f), (g) and (h) refer to the volumetric water content simulated by SDM, DPMDy, DPM-0.03 and DPM-0.01, respectively; and panel (i) refers to the crack ratio simulated by DPMDy. “Sim.” means “simulated”, and “Meas.” means “measured”.

Appendix B: Notation

PF-DC	Preferential flow induced by desiccation cracks
SDM	Single-domain model
EMs	Explicit crack models
DPoM	Dual-porosity model
DPM	Rigid dual-permeability model with fixed crack ratio and hydraulic conductivity
DPM-0.01	Rigid dual-permeability model with crack ratio of 0.01
DPM-0.03	Rigid dual-permeability model with crack ratio of 0.03
DPMDy	Dynamic DPM with changing crack ratio and hydraulic conductivity
WD cycles	Wetting–drying cycles
θ	Total water content (combined matrix and crack domains) (m3 m−3)
θexp	Volumetric water content measured in the experiment (m3 m−3)
θm	Volumetric water content of the matrix domain (m3 m−3)
θc	Volumetric water content of the crack domain (m3 m−3)
θm,s	Saturated volumetric water content of the matrix domain (m3 m−3)
θm,r	Residual volumetric water content of the matrix domain (m3 m−3)
θc,s	Saturated volumetric water content of the crack domain (m3 m−3)
θc,r	Residual volumetric water content of the crack domain (m3 m−3)
Se,m	Saturation degree of the matrix domain (m3 m−3)
Se,c	Saturation degree of the crack domain (m3 m−3)
αm	Parameter for the van Genuchten water retention curve of the matrix domain (1 m−1)
nm	Parameter for the van Genuchten water retention curve of the matrix domain (1 m−1)
mm	Parameter for the van Genuchten water retention curve of the matrix domain (1 m−1)
αc	Parameter for the van Genuchten water retention curve of the crack domain (1 m−1)
nc	Parameter for the van Genuchten water retention curve of the crack domain (1 m−1)

mc	Parameter for the van Genuchten water retention curve of the crack domain (1 m−1)
hm	Pressure head of the matrix domain (m)
hc	Pressure head of the crack domain (m)
Cc	Specific water capacity of the crack domain, which is defined as $d{\mathit{\theta }}_{\mathrm{c}}/d{h}_{\mathrm{c}}$ (1 m−1)
Cm	Specific water capacity of the matrix domain, which is defined as $\mathrm{d}{\mathit{\theta }}_{\mathrm{m}}/\mathrm{d}{h}_{\mathrm{m}}$ (1 m−1)
Ks	Total transient saturated hydraulic conductivity of the soil (combined matrix and crack domains) (m s−1)
Kc	Transient hydraulic conductivity of the crack domain (m s−1)
Kc,s	Saturated hydraulic conductivity of the crack domain (m s−1)
$K_{\mathrm{c},max}$	The maximum crack hydraulic conductivity when the crack reaches its maximum crack aperture (m s−1)
$K_{\mathrm{c},min}$	The minimum crack hydraulic conductivity when the crack reaches its minimum crack aperture (m s−1)
Kc,r	Relative hydraulic conductivity of the crack domain (m3 m−3)
Km	Transient hydraulic conductivity of the matrix domain (m s−1)
Km,s	Saturated hydraulic conductivity of the matrix domain (m s−1)
$K_{\mathrm{m},max}$	The maximum matrix hydraulic conductivity prior to soil shrinkage (m s−1)
Km,r	Relative hydraulic conductivity of the matrix domain (m3 m−3)
Ka	Hydraulic conductivity between the matrix and crack domains (m s−1)
$K_{{\mathrm{a}}_{min}}$	An improved hydraulic conductivity between the matrix and crack domains reformulated by Gerke et al. (2013) (m s−1)
Γw	Water exchange term between the crack and matrix domains (1 s−1)
wc	Crack ratio, which is defined as volumetric ratio between the crack domain and the overall soil volume (m3 m−3)
wc,exp	Surface crack ratio measured in the experiment (m2 m−2)
wj,exp	Average crack aperture (or crack width) measured in the experiment (m)
$w_{j,max}$	The maximum average crack aperture measured in the experiment (m)
dmax	The maximum crack depth measured in the experiment (m)
wm	Volumetric ratio between the matrix domain and the overall soil volume (m3 m−3)
αw	Effective water transfer coefficient (1 m−2)
V	Total soil volume (combined matrix and crack domains) (m3)
Vm	Volume of the soil matrix domain (m3)
Vc	Volume of the crack domain (m3)
Vp	Total pore volume (m3)
Vp,m	Pore volume in the matrix domain (m3)
Vp,c	Pore volume in the crack domain (m3)
ε	Total soil porosity (combined matrix and crack domains), which is defined as $V_{\mathrm{p}}/V$ (m3 m−3)
εm	Effective porosity of the matrix domain, which is defined as $V_{\mathrm{p},\mathrm{m}}/V_{\mathrm{m}}$
εc	Effective porosity of the crack domain, which is defined as $V_{\mathrm{p},\mathrm{c}}/V_{\mathrm{c}}$
i	Total effective infiltration rate (combined matrix and crack domains) (m s−1)
im	Effective infiltration rate of the matrix domain (m s−1)
ic	Effective infiltration rate of the crack domain (m s−1)
em	Effective evaporation rate of the matrix domain (m s−1)
ec	Effective evaporation rate of the crack domain (m s−1)
r	Rainfall intensity (m s−1)
AE	Actual evaporation rate (m s−1)
PE	Potential evaporation rate (m s−1)
S	Total matric suction at the soil surface (kPa)
Sexp	Soil matric suction measured in the experiment (kPa)
g	Gravitational acceleration constant (m s−2)
ωv	Molecular mass of water (kg mol−1)
ξ	Dimensional empirical parameter with a suggested value of 0.7
ha	Relative humidity of soil overlying air
γw	Unit mass of water (kN m−3)
R	Universal gas constant J (mol K)−1
Ts	Soil surface temperature (∘C)
ϕmax	Total porosity (or the maximum porosity) of a soil core prior to soil shrinkage, which is defined as $V_{\mathrm{p}}/V$ and thus is equal to the ε (m3 m−3)
ϕmin	The minimum porosity of the matrix domain (m3 m−3)
ϕmatrix	Porosity of the matrix domain, which is defined as $V_{\mathrm{p},\mathrm{m}}/V$ (m3 m−3)
ϕcrack	Porosity of the crack domain, which is defined as $V_{\mathrm{p},\mathrm{c}}/(V_{\mathrm{m}}+V_{\mathrm{c}})$ (m3 m−3)
ϕsub	Porosity of the subsidence zone, which is defined as voids induced by soil subsidence divided by the total soil volume (m3 m−3)
U	A unified water content, which is defined as the gravimetric water content u divided by its saturated value umax
p	Functional shape parameters of the soil shrinkage curve
q	Functional shape parameters of the soil shrinkage curve
v	Water kinematic viscosity (m2 s−1)
tp	Beginning of ponding time after each rainfall event (min)
Δh	Pressure difference between the crack and matrix domains, which is defined as hc−hm