New Analytical Model for Cumulative Infiltration into Dual-Permeability Soils

Water quality must be protected because water is used for many different purposes (e.g., drinking, irrigation, and industry) and provides a living environment for marine, aquatic, and continental ecosystems, yet climate change is predicted to drastically impact the water cycle and reduce water resources in many terrestrial environments, requiring proper management of their quantity and quality (McDonald et al., 1996; Gascoin et al., 2009). Additionally, there is a large spectrum of chemical substances that are released into the environment by human activities and that constitute a risk to groundwater quality. In particular, in urban areas, runoff water carries significant loads of pollutants, including heavy metals, hydrocarbons, pesticides, bacteria, and nutrients (Mikkelsen et al., 1997).

In recent decades, field and laboratory studies have demonstrated that water flow and solute transport in the unsaturated zone often occur through a small fraction of the soil along preferential flow paths (Flury et al., 1994; Buttle and Leigh, 1997; Allaire-Leung et al., 2000a, 2000b). Preferential flow reduces the access of pollutants to the soil matrix and thus their opportunity to adsorb onto soil particles (Lassabatere et al., 2004; Lassabatere et al., 2007; Lamy et al., 2009). Therefore, there is an urgent need for scientific advances in developing new measurement techniques and theoretical approaches to better understand preferential water flow and physical nonequilibrium solute transport, particularly for the case of water infiltration into the soil and related pollutant transport.

Physically based modeling concepts for water flow and solute transport in soils can be broadly classified into continuum, dual- and multi-continuum, and network approaches (Köhne et al., 2009a, 2009b). In a single-continuum approach, the traditional conceptualization of soil is that of a single-permeability porous medium, usually with unimodal water retention and hydraulic conductivity curves (Šimůnek and van Genuchten, 2008). The dual-permeability (DP) approach assumes that the soil has structural pores representing a fast-flow region and the soil matrix, with a much lower or even zero saturated hydraulic conductivity (Šimůnek and van Genuchten, 2008; Köhne et al., 2009a, 2009b). The transfer term representing the lateral exchange of water and/or solute between the two regions is usually described using the first-order process (van Genuchten and Dalton, 1986; Gerke and van Genuchten, 1993). The dual-permeability approach of Gerke and van Genuchten (1993) assumes that flow in both regions obeys the Buckingham–Darcy law. Such an assumption can be invalid in situations where soils have large macropores, cracks, or fractures, and when specific approaches are required to account for nonlaminar and turbulent flow in these structural features (Coppola et al., 2009b). The focus of this study, however, was dual-permeability soil, for which the same conceptual approach can be used to describe flow and transport in both fast and slow domains.

With regard to water infiltration into single-permeability soil, many analytical models have been obtained by analytically solving the Richards equation, using the concept of soil sorptivity and flux concentration (for a detailed review, see Haverkamp et al., 2006). In particular, Haverkamp et al. (1994) proposed a set of analytical equations for water infiltration from a disk source, assuming a constant pressure head at the soil surface and a uniform initial pressure head profile. This model is valid for any combination of zero or negative water pressure head values in the soil profile and at the surface when the water pressure head at the soil surface is higher than in the soil profile. Their analytical model was validated against both experimental infiltration data (Clausnitzer et al., 1998) and numerically generated data (Lassabatere et al., 2009). This analytical model led to several methods for the determination of the soil hydraulic properties of single-permeability soils from water infiltration experimental data, e.g., the BEST method (Lassabatere et al., 2004; Yilmaz et al., 2010).

In this study, we developed an extension of the model of Haverkamp et al. (1994) for cumulative water infiltration into dual-permeability soils. This extended model was based on several assumptions. First, we assumed that the total cumulative infiltration at the soil surface is not affected by exchange between the two pore domains. This means that cumulative infiltration into a dual-permeability soil with no water exchange between the matrix and fast-flow regions can be used to represent infiltration into a dual-permeability soil with both slow and fast exchange between the two regions. In the case of no water exchange, water infiltration into a dual-permeability soil implies independent infiltration into two separate single-permeability (SP) systems. The total water infiltration at the soil surface is then evaluated as a simple linear combination of cumulative infiltration in the fast-flow and matrix regions, considering the volumetric fractions occupied by each region. We then assumed that each infiltration into two single-permeability domains can be properly described using the analytical model of Haverkamp et al. (1994), provided its shape parameters are adjusted according to the soil type as suggested by Lassabatere et al. (2009). Finally, the proposed model was tested against its capability to well describe infiltration into dual-permeability soils with significant water exchange between the matrix and fast-flow regions and to simulate water infiltration experiments (direct modeling) and invert synthetic and real experimental data.

We first present the theoretical background of the proposed model. Using numerical simulations, we validate the assumption that water exchange between the soil matrix and fast-flow regions has no significant effect on the overall surface water infiltration. The proposed analytical model is then used to model single- and multi-tension water infiltration experiments for several synthetic dual-permeability media. This step helps in assessing the effect of the fast-flow region on water infiltration. Finally, we assess the capability of the proposed model to use either generated or experimental data in the inverse mode to estimate relevant soil hydraulic parameters. In this study, we first use analytically generated data using both SP and DP approaches to compare the accuracy of model predictions and the precision of the obtained parameter estimates. Finally, we assess the capability of the proposed model to invert the experimental data obtained for two coarse materials for which dual-permeability behavior was expected.

For this study, we considered several synthetic DP soils, which were formed as a combination of several types of matrices and fast-flow regions. Fast-flow regions represent several types of porous media with mesopores or small macropores. As a first step, we calculated wetting fronts and cumulative infiltration as a function of interfacial saturated hydraulic conductivity, which is the parameter that governs water exchange between the matrix and fast-flow regions (Step 1). Once the hypothesis that inter-region water exchange has a negligible effect on infiltration was validated, we applied the proposed analytical model to the same synthetic DP soils to calculate cumulative infiltration as a function of the initial and surface pressure heads (Step 2). Analytically generated data were then inverted, using the proposed DP or SP models, and the accuracy of model descriptions and estimated parameters was determined (Step 3). Finally, we inverted two experimental data sets using the proposed model (Step 4). Below, we first describe the synthetic and experimental soils and then the numerical and analytical methods used for both direct and inverse modeling.

The pore sizes were chosen to be in the range of small macropores and mesopores (Luxmoore, 1981). We assumed that this choice does not contradict the validity of the Buckingham–Darcy law. All hydraulic parameters are given in Table 1, and corresponding DP soil hydraulic functions (water retention and hydraulic conductivity functions) are displayed in Fig. 1.

Two experimental soils were selected due to their coarse texture and expected bimodal pore-size distribution. The first material was a bimodal porous material used extensively for commercial green roofs (Société des Terreaux et Amendements Rochelais, Forges d’Aunis, France). The growth medium was composed of 70% mineral fraction and 30% organic matter. The mineral fraction was a mixture of 3- to 20-mm pumice stones and 3- to 7- and 7- to 15-mm pozzolanic ash. The organic matter consisted of composted maritime pine bark, and 3- to 20-mm black and blond peat. The material specific density, ρ_s, was 1.88 ± 0.04 g cm⁻³ measured by the pycnometer method. The dry bulk density, ρ_b, was 0.60 ± 0.06 g cm⁻³, leading to a porosity on the order of 68.1%. The field capacity was evaluated to be 49.4%, indicating that a significant amount of water was stored in the porosity inside the particles of the material.

The second material was a basic oxygen furnace (BOF) slag, which was studied by Yilmaz et al. (2010, 2013) with respect to the evolution of its hydraulic properties as a result of weathering and clogging processes. The material was produced in March 2003 at the site of Fos-sur-Mer in the south of France. It has a narrow and coarse particle size distribution, with a specific density of 2.97 g cm⁻³. Additional information about the geochemical aspects of the material characteristics were provided by Yilmaz et al. (2010) and Deneele et al. (2008).

The influence of interfacial hydraulic conductivity, K_s,a, on cumulative infiltration was first investigated for the DP soils evaluated in this study. The effect of capillary and gravity forces was investigated step by step by solving Eq. [1–4], first without gravity to model horizontal infiltration due only to capillary forces, and then with gravity to model vertical infiltration due to both capillary and gravity forces.

One-dimensional water flow for horizontal and vertical infiltration was modeled using HYDRUS-1D (Šimůnek et al., 2008). Water infiltration was modeled for initially dry soils with an initial pressure head of −10 m and ponding conditions at the soil boundary (h_surf = 0 and θ_surf = θ_s). The initial water content and pressure head profiles were assumed to be uniform (θ₀, h₀), as required by the model of Haverkamp et al. (1994), i.e., Eq. [5–7]. The zero pressure head at the boundary was chosen to maximize the role played by macroporosity. Calculations with different values of the initial pressure heads led to similar conclusions.

The finite-element mesh and other numerical options used in HYDRUS-1D simulations were similar to those used by Lassabatere et al. (2009). The flow domain was assumed to be 1 m deep for loam and silt and 40 cm deep for silty clay. Finite elements of 2 mm were used to discretize the transport domain. The maximum time was 90 min, which corresponds to the average duration of water infiltration experiments. The minimum time step was 10⁻³ min for loam and silt and 10⁻⁴ min for silty clay. Precision tolerances for pressure heads and water contents were fixed for all soils at HYDRUS default values of 10 mm and 0.001 m³ m⁻³, respectively.

The impact of water exchange at the interface between the matrix and fast-flow regions on the water infiltration rates, cumulative infiltration, and wetting fronts in the system is evaluated below. For each synthetic DP soil, the water exchange was evaluated using Eq. [2–4] and assuming the HYDRUS default values, i.e., β_a = 15, d = 1 mm, and γ_w = 0.4. The interfacial conductivity, K_s,a, was first assumed to be zero (no exchange) and then gradually increased from 10⁻¹¹ (low exchange) to 10⁻³ mm min⁻¹ (almost instantaneous water exchange between the matrix and fast-flow regions).

Soil sorptivities were calculated using Eq. [6] for selected hydraulic conditions and soil hydraulic parameters (Table 1). Sorptivity values were then used to calculate cumulative one-dimensional infiltration into the two regions using Eq. [8]. The shape parameter, β, was derived according to Lassabatere et al. (2009), who optimized β by fitting Eq. [5] to numerically generated data for four types of soils (silty clay, silt, loam, and sand). They showed that optimized values of the shape parameter β vary significantly among different soil types, with a minimum value of 0.344 for sand (coarser texture), a maximum value of 1.65 for silty clay (finer texture), and a value of 1.56 for silt. The values of β considered for the synthetic soils are given in Table 1.

Calculations were performed for conditions representing single- and multiple-tension water infiltration experiments. For the single-tension infiltration (STI) experiments, while several values were considered for initial pressure heads, only a zero pressure head was considered at the soil surface, as in the widely used Beerkan experiment (Braud et al., 2005; Lassabatere et al., 2006, 2010). Multi-tension infiltration (MTI) experiments were modeled as successive infiltration steps with the following boundary pressure heads: −150, −60, −30, and 0 mm. These values were proposed in several previous studies (Watson and Luxmoore, 1986; Timlin et al., 1994; Malone et al., 2004) and are often considered as a guide for water infiltration experiments using tension disk infiltrometers (Mallants et al., 1997; Jacques et al., 2002).

The cumulative MTIs were calculated as follows. The first infiltration was calculated using Eq. [8] with the initial pressure head h₀ and the lowest value of the boundary pressure head. At the end of the first infiltration, the second infiltration was calculated, considering as the initial pressure head the value previously imposed at the boundary and the second boundary pressure head. The same procedure was used for the third, fourth, and other infiltration steps. Total cumulative infiltration was then derived as a combination of all cumulative infiltration steps. Time and/or infiltration volumes used at each of the infiltration steps were selected according to several types of infiltration protocols (see below).

We used the analytically generated data to assess the capability of the proposed model in inverse mode and to compare the accuracy of model descriptions and precision of obtained parameter estimates for both SP and DP approaches. Afterward, both SP and DP approaches were fitted to the experimental data obtained for the two soils described above. For these soils, all water infiltration experiments were conducted using the SW-080 infiltrometer developed by SDEC (Deurer et al., 2008). A series of successive tensions was applied during the infiltration experiment: −155, −130, −70, −40, −20, −10, and −1 and −160, −130, −100, −70, −40, −10, −1 mm for two replicates of the growth material and −110, −50, −20, and 0 mm for the BOF slag studied by Yilmaz et al. (2010, 2013). Small increments in applied tensions were used to maximize the information obtained during the infiltration experiment. The data extracted from Yilmaz et al. (2013) were obtained 4 mo after embedding of the BOF slag. The data were initially unsuccessfully analyzed using the HYDRUS-3D inverse modeling procedure while considering a SP approach.

The impact of water exchange at the matrix–fast-flow interface on infiltration rates, cumulative infiltration, and wetting fronts was evaluated for all combinations of synthetic soil matrices and fast-flow regions listed in Table 1 and for both horizontal and vertical infiltration. Because the flow enhancement due to the presence of the fast-flow region is larger when flow is in the vertical direction, the effect of the interfacial conductivity, K_s,a, on the wetting fronts and cumulative infiltration is also larger for these flow conditions.

The effect of K_s,a on the wetting fronts is depicted in Fig. 2 for vertical flow. When the interfacial conductivity is zero, there is no exchange of water between the matrix and the fast-flow region (Fig. 2a). Water separately and independently infiltrates into the matrix and the fast-flow region without any interaction between the two regions. This produces a much deeper moisture front in the fast-flow region (Fig. 2a). When the interfacial conductivity is intermediate (Fig. 2b), water infiltration is faster into the fast-flow region, but a fraction of infiltrated water moves laterally from the fast-flow region to the matrix due to a gradient in the water pressure head. Such exchange is reflected in the second, more dispersed, wetting front ahead of the main wetting front in the soil matrix (Fig. 2b). Water infiltration in both regions is now linked. When the interfacial hydraulic conductivity is very high, water exchange seems to instantaneously balance the water pressure heads between the matrix and the fast-flow regions (Fig. 2c). Consequently, wetting fronts in the two regions move downward at almost the same rate. In this case, cumulative total infiltration at the soil surface corresponds to that of the single-permeability model with bimodal water retention and hydraulic conductivity functions, as depicted in Fig. 1. Clearly, the interfacial hydraulic conductivity between the matrix and fast-flow regions significantly affects the movement of water in the DP porous medium.

The effect of the interfacial hydraulic conductivity on cumulative infiltration is illustrated in Fig. 3. Although this figure shows results only for horizontal infiltration, the same conclusions can also be obtained from an analysis of vertical infiltration. The interfacial hydraulic conductivity impacts cumulative infiltration in both the matrix and the fast-flow region (Fig. 3). Any increase in K_s,a causes an increase in cumulative infiltration in the fast-flow region (Fig. 3g–3i) and a decrease in cumulative infiltration in the matrix (Fig. 3d–3f). Water flowing through the fast-flow region laterally infiltrates into the matrix, and the fast-flow region thus acts as a leaking tank. This loss of water lowers the pressure head in the fast-flow region, increases the water pressure gradient, and correspondingly also increases infiltration at the surface. The opposite occurs in the matrix. Water from the fast-flow region increases water contents in the matrix, reduces the pressure head gradient, and correspondingly also reduces the infiltration rate. The total infiltration is the result of the summation of an increase in the fast-flow region and a decrease in the matrix. These two opposing effects compensate each other and produce only small differences in the total infiltration (Fig. 3a–3c). Total cumulative infiltration is thus relatively independent of the interfacial hydraulic conductivity.

Cumulative infiltration into DP soils with K_s,a larger than zero are compared with reference to the total cumulative infiltration into DP soils with K_s,a equal to zero in Table 2, which lists time-averaged relative differences for all matrix and fast-flow region combinations. These differences rarely exceed 15%, with most of them being <5% (in absolute values). Cumulative infiltration for soils with different interfacial conductivities are clearly similar to those with zero interfacial conductivity. As a result, total cumulative infiltration can be approximated by a linear combination of cumulative infiltration into the matrix and the fast-flow regions (Eq. [8a]). Assuming that infiltration in each region can be described using the analytical model developed by Haverkamp et al. (1994), the proposed model (Eq. [8]) can be considered to be valid.

The analytical model, i.e., Eq. [8], was used to describe cumulative infiltration for single- and multi-tension infiltration experiments into synthetic SP and DP soils and to analyze real experimental data. Because the analytical results were similar for all combinations of synthetic DP soils, they are depicted only for silt. Cumulative STIs are depicted for different combinations of initial and surface water pressure heads (Fig. 4). The initial water pressure head had a similar effect on infiltration into both DP and SP soils. An increase in the initial water pressure head resulted in a decrease in the hydraulic gradient [difference (h_surf − h₀)] and a small decrease in cumulative infiltration, yet the differences remained relatively small (Fig. 4a, 4c, and 4e). Although the water pressure head at the soil surface (h_surf) impacted water infiltration into both SP and DP soils, it impacted flow enhancement by the fast-flow region much more. When the water pressure head at the surface was zero (Fig. 4a, 4b, 4c, and 4e), water infiltration into DP soils was higher than infiltration into SP soils. When the surface water pressure head was sufficiently high to activate flow into the fast-flow region, it significantly increased total cumulative infiltration into the entire DP system. When the water pressure head at the soil surface was reduced, infiltration into the fast-flow region decreased until it became even lower than infiltration into the matrix (Fig. 4d and 4f). This resulted from the fact that the fast-flow region is characterized by low water retention and a very low hydraulic conductivity for unsaturated conditions. These results are in full agreement with previous field observations and laboratory experiments (Durner, 1994; Rousseau et al., 2004; Majdalani et al., 2007; Lamy et al., 2009).

For MTI experiments, two different procedures were investigated: (i) the same infiltration volumes for all boundary water pressure heads (CV for constant volume) and (ii) the same infiltration duration (CT for constant time). For both CV and CT procedures, the same total amount of water was supposed to infiltrate, corresponding to the reservoir capacity of the SDEC infiltration disk SD-080 infiltrometer (Deurer et al., 2008). The set of imposed boundary pressure heads was equal to −150, −60, −30, and 0 mm, as suggested by Watson and Luxmoore (1986). When MTI experiments were performed following the CT procedure, infiltration curves for the synthetic SP and DP soils were almost indistinguishable for the first three suctions and strongly deviated for the last suction, i.e., for zero water pressure head (Fig. 5a). The slope of the last part of the cumulative infiltration curve reflects the saturated hydraulic conductivity of a particular system.

For the CV procedure, visible differences appeared even for the first boundary suctions (Fig. 5b). In this case, infiltration of identical volumes requires longer times for lower boundary pressure heads, giving more weight to these parts of the infiltration curves. For lower pressure heads, all synthetic DP soils had similar cumulative infiltration because only the soil matrices were activated. When the surface water pressure head was increased to −60 mm, the pores in the fast-flow region of the DP-r_min soil were also activated, and the infiltration curve for this system separated from the other curves (Fig. 5b, DP-r_min). When the surface water pressure head was increased to −30 mm, the pores of the fast-flow region in the DP-r_mean soil were activated, and the corresponding infiltration curve (Fig. 5b, DP-r_mean) started deviating from that of the DP soil with the largest pores. Finally, zero pressure head was needed to activate the fast-flow region of the DP-r_max soil and to increase the infiltration rate into this system (Fig. 5b, DP-r_max). These results show that the CV procedure is preferable for identifying the successive activation of pores associated with the fast-flow region. The differences between infiltration curves in Fig. 5a and 5b clearly demonstrate that the shapes of cumulative infiltration are deeply affected by the choice of the experimental protocol.

Infiltration rates were computed as well, and their relative values at the end of each infiltration period are reported against the surface pressure head (Fig. 5c and 5d). For the SP soils, the increase in the surface pressure head led to a continuous increase in final infiltration rates. For the DP soils, the same increase led to a larger increase in infiltration rates due to the activation of pores in the fast-flow region. Such a large increase in the infiltration rate for DP soils close to saturation has already been reported (Deurer et al., 2008). It is interesting to note that the graphical representations exhibit a similar pattern for both protocols when the infiltration rate is considered (Fig. 5c and 5d), whereas the shape greatly differs when cumulative infiltration is considered (Fig. 5a and 5b). In addition, such a plot allows us to clearly separate the SP and DP systems and to sort the DP systems as a function of the size of pores associated with the fast-flow region. These results indicate that one of the ways to deal with infiltration data could be to work with infiltration rates instead of cumulative infiltration, yet the estimation of experimental rates that is based on the mathematical derivation of experimental cumulative infiltration requires an extreme precision, potentially unreachable when using usual devices. Because of such a difficulty, most studies use and invert cumulative infiltration (Angulo-Jaramillo et al., 2000).

Next we address the feasibility of detecting dual-permeability behavior using inversion procedures. The cumulative infiltration results obtained for the DP systems were inverted using both SP and DP models. The analysis of errors allows assessments of the capability of the SP and DP models to reproduce the DP data.

Similarly to the results above, similar results were obtained for all DP systems. Figure 6 illustrates the results for the DP system with a silty matrix for (i) STI (with h_surf = 0) and (ii) MTI experiments. For the STI experiments, the SP model could accurately fit the infiltration data (Fig. 6a) with small errors (Fig. 6c and 6e) and with the NSE close to 1 and a CVRMSE of <1% (Table 3). These fits were obtained with high values of the saturated hydraulic conductivity and high values of the parameter α (Table 3). This means that the inversion of infiltration data for the DP soil, assuming the SP model, will provide an accurate fit of the data, with no visible differences compared with the DP approach, but having wrong estimates of the soil hydraulic parameters. Given the fact that both SP and DP approaches provide a similar goodness of fit, it is impossible to discriminate which model is the correct model, and the dual-permeability behavior of the soil can be overlooked.

When using MTI data, the fits remained quite accurate for the synthetic DP-r_min system but clearly less accurate for the DP-r_mean and DP-r_max systems, i.e., systems with larger pores in the fast-flow region (Fig. 6b and 6d). Lower values of NSE and larger values of RMSE and CVRMSE confirmed these observations (Table 3). Clearly, the fits underestimated cumulative infiltration for the first and last suctions and overestimated it for the two intermediate suctions (Fig. 6f). This reflects the inadequacy of the SP model to reproduce the large increase in cumulative infiltration close to saturation that resulted from the activation of the fast-flow region. Similar conclusions could be derived when inverting the CT data for the DP-r_min and DP-r_mean systems. No fits could be obtained when inverting data for the DP-r_max system. Meanwhile, the DP approach provided an accurate description of the data, with NSE close to unity and very low values for RMSE and CVRMSE. This means that the SP and DP approaches can be distinguished in terms of their accuracy in describing the data, and a proper choice can be made.

These results indicate that inverting MTI experimental data using both SP and DP approaches allows the detection of dual-permeability behavior, whereas this is not the case for STI experiments. These results reflect the fact that the MTI procedure contains more information than the STI procedure and must be considered for the detection of dual-permeability behavior. The lack of information in the STI data was already pointed out by several studies, even for the determination of unimodal hydraulic curves (e.g., Šimůnek and van Genuchten, 1996, 1997). In addition, inverting MTI data with the DP approach led to accurate estimates of parameters of both the matrix and fast-flow regions.

We investigated above the adequacy of the proposed model to invert analytically generated data, yet experimental data contain additional sources of error (due to measurements, treatment of data, etc.), which may lead to the failure of the parameter estimation process. The proposed DP approach was thus assessed with regard to its capability to invert real experimental data. Figures 7 and 8 show the results of the analysis for the growth medium and the BOF slag, respectively.

The two replicates of the growth material exhibit similar cumulative infiltration (Fig. 7a and 7b). For the first replicate, both SP and DP approaches provided accurate fits (Fig. 7a, 7c, and 7e). Similarly, the values of the NSE, RMSE, and CVRMSE were on the same order, indicating no significant increase in the quality of the fit when considering the DP approach rather than the SP approach (Table 4, Growth Medium Replicate no. 1). In addition, the hydraulic properties did not differ significantly between the matrix and fast-flow regions (Table 4, Growth Medium Replicate no. 1, DP fit), with K_s in the fast-flow region being only ≈20 times larger than K_s in the matrix. Moreover, the matrix of the DP system resembled the SP model in terms of hydraulic properties (Table 4, Growth Medium Replicate no. 1, DP fit vs. SP fit). From these results, it can be concluded that the experimental system may correspond to a single-permeability medium.

For the second replicate, larger differences between the SP and DP approaches appeared in terms of errors and average errors (Fig. 7d and 7f). The values of NSE were slightly higher and RMSE and CVRMSE slightly lower for the DP model (Table 4, Growth Medium Replicate no. 2). The optimized values seem to indicate a significant increase in the saturated hydraulic conductivity and an increase in the parameter α with a fast-flow region content around 2.7%, yet it is difficult to make conclusions on the basis of the analysis of errors and model precision. These results resemble the inversion results of data analytically generated using the DP model with small pores in the fast-flow region using the SP model (Table 3; Fig. 6, MTI, DP-r_min). It seems difficult to distinguish between (i) a SP system with a relatively high saturated hydraulic conductivity and a high value of α and (ii) a DP system with a matrix exhibiting usual values of hydraulic parameters and a fast-flow region with a high value of saturated hydraulic conductivity and a high value of parameter α. As a result, one cannot conclude that this material exhibits a dual-permeability behavior.

For the BOF slag, the increase in the model adequacy (from the SP to DP approach) seems more significant, with a significant increase in the NSE and a decrease in the RMSE and CVRMSE. Fits were more accurate using the DP approach (Fig. 8b vs. 8a), with randomly distributed small errors (Fig. 8c) and suction-averaged errors close to zero (Fig. 8d). The SP approach exhibited larger errors and suction-averaged errors, similarly to when analytically generated data for the DP model with large pores in the fast-flow region were inverted using the SP approach (Fig. 6, DP-r_max). These errors highlight the inadequacy of the SP approach and clearly indicate a dual-permeability behavior. The estimated hydraulic parameters for the BOF slag indicate that macroporosity occupied approximately 10% of the DP system. The fast-flow region had a saturated hydraulic conductivity about 10 times larger than the matrix, the α parameter was about seven times larger, and the shape parameter n was almost three times larger, indicating the behavior of a very coarse material (Schaap et al., 2001).

This last example illustrates how the dual-permeability approach can be needed to accurately describe experimental data. Yilmaz et al. (2013) studied the evolution of the BOF slag hydraulic properties due to weathering and concluded that its porosity was clogged and that the saturated hydraulic conductivity was reduced due to clogging. They also discussed the lack of information in the ponded infiltration data in comparison with multiple tension disk infiltrations, however, and could not obtain accurate fits for some of the multiple tension disk infiltration experiments (Yilmaz et al., 2013, Fig. 2a—data related to 4 mo). In this study, the dual-permeability approach provided a perfect fit of the data, whereas the SP approach failed. This may indicate the appearance of a dual-porosity system in the BOF slag with weathering. Such results are consistent with previous studies that proved that weathering processes, and in particular the dissolution of portlandite followed by the precipitation of calcite, could clog intermediate pores, leading to a bimodal porosity with the concomitancy of large and small pores (Yilmaz et al., 2013). These results indicate the need to properly model the pore size distribution and to account for dual porosity (or dual permeability) when studying these types of materials.

We have presented an analytical model for cumulative infiltration into dual-permeability soils. The total water infiltration was evaluated as a simple linear combination of cumulative infiltration into the fast-flow and matrix regions. Each infiltration was described using the analytical model of Haverkamp et al. (1994), with adjustment of its shape parameters as a function of the soil type. We now (i) synthesize and discuss the main analytical and numerical results with regard to the physics of water flow in soils, (ii) propose a methodology to detect the dual-permeability behavior in the context of the hydraulic characterization of soils, and finally, (iii) address the shortcomings of the proposed model and discuss ongoing research.

The numerical and analytical results showed that the fast-flow region enhances water infiltration when surface pressure heads are high enough to activate pores of the fast-flow region, as already mentioned in the literature (Allaire et al., 2009; Köhne et al., 2009a, 2009b). We also showed that this enhancement of infiltration does not depend on water exchange at the interface between the matrix and the fast-flow regions but only on the global hydraulic conductivity of the system. Water exchange affects only the wetting fronts in the soil profile.

These results are important with respect to modeling flow and pollutant transport. When only water flow is of interest, then the single-permeability model with global soil hydraulic properties can be used, as already suggested by Durner (1994) and Zurmühl and Durner (1998). However, the dual-permeability approach is needed for an accurate description of pollutant transport (Dousset et al., 2007), because preferential flow not only restricts the access of solutes to sorption sites located in the pockets of stagnant water (Fesch et al., 1998a, 1998b; Padilla et al., 1999; Lassabatere et al., 2004; Lamy et al., 2009) but also increases the pore water velocity of the mobile water fraction, thus reducing the contact time between solutes and sorption sites and limiting the effect of kinetic processes (Lassabatere et al., 2007; Hanna et al., 2010, 2012). These are some of the reasons why a dual-permeability behavior must be detected and fully identified when the fate of pollutants is of interest.

Water infiltration experiments constitute one of the main experimental means to characterize soil hydraulic properties (Angulo-Jaramillo et al., 2000). However, there is a need to fill in the gaps in our knowledge regarding the detection and characterization of dual-permeability soils (Coppola et al., 2009a; Bagarello et al., 2013). The numerical and analytical results of this study provide guidance on how to detect such behavior.

First, it was concluded that the interfacial hydraulic conductivity has a negligible impact on cumulative infiltration. This means that this parameter cannot be estimated by fitting cumulative infiltration data. For that, specific information about the location of wetting fronts and water exchange between the matrix and fast-flow regions would be needed. Tracer or dye experiments could be used to quantify solute and water exchange, as indicated by Buttle and Leigh (1997), Angulo-Jaramillo et al. (2000), or Castiglione et al. (2003, 2005).

Second, experimental protocols must be adequately designed to avoid problems of non-uniqueness and ill-posed estimates. Single-tension infiltration was shown to be insufficient for detecting dual-permeability behavior, as already reported by Olioso et al. (2002) and Šimůnek and van Genuchten (1997). In our study, we clearly proved that experimental protocols should be based on multi-tension infiltration experiments, with values of the applied tensions being a function of the expected pore size associated with the fast-flow region. In addition, the application times should be longer for higher tensions to facilitate the detection of a dual-permeability behavior.

Third, specific mathematical developments are needed to improve the characterization process. In this study, we used the quasi-exact formulation for dual-permeability soils as defined by Eq. [8], yet these equations can be complicated to compute and their inversion can lead to problems with nonconvergence for experimental data with substantial errors. For single-permeability media, the complex implicit formulation, i.e., Eq. [5–7], was simplified to obtain approximations for short and long times (Lassabatere et al., 2009). These approximations were incorporated into the BEST method for derivation of the whole set of hydraulic parameters from single-tension experiments with a zero surface water pressure head (Lassabatere et al., 2006; Yilmaz et al., 2010, 2013). The same approach could be considered for the derivation of an adapted BEST method for dual-permeability soils provided that: (i) approximations of Eq. [8] are proposed for short and long times and then implemented into the updated BEST algorithm, (ii) the updated BEST algorithm is adapted for the analysis of infiltration tests with negative surface pressure heads and for multi-tension infiltration, and (iii) complementary information is obtained and used together with cumulative infiltration data to derive the interfacial hydraulic conductivity. This is the subject of ongoing research.

The proposed model can be questioned with regard to its validity and conditions of use. Indeed, the proposed model inherits restrictions of its application from the Haverkamp model for single-permeability media (Haverkamp et al., 1994; Smettem et al., 1994). This means that it can be used for infiltration experiments with a constant initial water pressure head and a constant surface water pressure head equal to zero or negative. The condition of a uniform initial pressure head profile is problematic, especially for successive infiltration. Clearly, such conditions may not be fulfilled in the case of a constant time procedure, when a small volume of water may infiltrate during the application of the first suction, followed by larger volumes during the applications of subsequent suctions. Finally, Haverkamp et al. (1994) recommended that the initial water content should be less than a quarter of the final water content to ensure the validity of their equation, which constitutes an additional restriction.

The use of the proposed model requires that the fast-flow region obey the Buckingham–Darcy law. This hypothesis seems plausible for soils with a fast-flow region consisting of mesopores or small macropores. On the other hand, specific approaches are required to account for non-laminar and turbulent flow in large macropores (Hincapié and Germann, 2009). Coppola et al. (2009b, 2009c) suggested distinguishing real macropores, i.e., non-capillary pores, from interaggregate pores, which can still be considered to be capillary pores. They recommended a combination of the Richards equation for capillary pores and the kinematic wave equation (Germann, 1990; Di Pietro et al., 2003) or the Hagen–Poiseuille equation (Ahuja et al., 2000; Malone et al., 2004) for non-capillary pores. Some authors, however, observed substantial discrepancies between the results of all these models and the experimental data obtained for solutes. For instance, Castiglione et al. (2003) and Lamy et al. (2009) revealed experimentally that large macropores not only funnel preferential flow but also enhance flow in the surrounding matrix. Water flow is much more complex than purely bimodal, with high fluxes in macropores and slow fluxes in the matrix. These observed flow patterns are inconsistent with a dual-permeability approach or any other approach discussed above. Further research is thus required in the area of modeling preferential flow and its impact on water infiltration and pollutant transfer.

We would like to acknowledge Usha Pillai (Sustainable Mining Institute, Centre for Mined Land Rehabilitation) for helpful discussion regarding soil science. We would like to thank the editor and the anonymous reviewers for their time and their constructive suggestions and helpful comments.