# Vadose Zone Journal

- Soil Science Society of America

## Abstract

Macropores are important hydrologic features that result in preferential flow and transport even under partially saturated flow conditions. The objective of this study was to use numerical simulations to investigate the hydraulic representation of preferential flow in partially saturated macroporous soils. Tension infiltration experiments that exhibited varying degrees of preferential flow, primarily along worm burrows, provided the basis for the numerical simulations. Field measurements of infiltration, soil water content, and dye transport were used to calibrate the model and assess the results. A three-dimensional model was constructed such that the soil matrix contained discrete vertical macropores. The simulations were able to capture the relevant flow and transport characteristics during infiltration. Sensitivity analyses demonstrated the utility of cumulative infiltration and dye transport data for constraining numerical simulations of macroporous systems. Hydraulic conductivity estimates for both matrix and macropores were lower than expected, which may be due to an overly simplified description of macropore flow hydraulics. Simulated macropore discontinuities near the surface reduced the infiltration volume by >50% and the depth of dye transport by >80%. The simulations also showed that increasing macropore density was nearly linearly related to increases in preferential flow, and confirmed field observations that closer macropore spacing led to increased transport depths due to macropore–matrix interaction and conjoined wetting fronts between neighboring macropores. This discrete macropore approach provides a useful method for examining macropore flow and transport, and highlights gaps in our understanding of the unsaturated flow behavior of macropores.

Water flow and dye tracer transport from tension infiltration experiments were simulated with a three-dimensional discrete macropore model. Macropore connectivity and density were important controls on preferential flow. Unrealistic hydraulic parameter estimates suggested complex flow in macropores and a need to better describe macropore hydraulics.

**Preferential water flow along macropores** can have a major influence on water and contaminant movement in field soils (Beven and Germann, 1982; Flury et al., 1994). Macropores develop as a result of soil desiccation, animal burrows, and root channels and are common in many fine-grained structured soils. The importance of flow and transport within macropores has been recognized in recent decades and has resulted in a substantial increase in research on preferential flow dynamics (Jarvis, 2007). Macropore flow is generally associated with larger rainfall or irrigation events and can increase the rate and depth of infiltration (Akay and Fox, 2007; Eguchi and Hasegawa, 2008), leading to rapid vertical transport of solutes (Sander and Gerke, 2007; Cey and Rudolph, 2009) and biocolloids (Unc and Goss, 2003; Guzman et al., 2009) as well as influencing the overall runoff response of watersheds (Joerin et al., 2005).

Numerical models have been developed to evaluate water and solute flow in macroporous systems (Šimůnek et al., 2003; Köhne et al., 2009). Conceptually, models are broadly classified as either multicontinuum or discrete macropore models. Dual-permeability models are a subset of multicontinuum models that consist of two interacting domains superimposed in space. One domain represents the primary porosity of the soil matrix and the second domain represents secondary porosity (i.e., macropores and fractures). Hydraulic properties of the fast flow paths are thus separated from the surrounding, less permeable matrix material, resulting in an efficient and stable numerical solution for preferential flow under partially saturated conditions (Gerke and van Genuchten, 1993). In addition, the specification of mass transfer terms is critical for dual-permeability approaches to describe water and solute coupling between domains (Šimůnek et al., 2003).

Discrete macropore models provide an alternative approach whereby macropores are explicitly represented within subregions of a model domain (e.g., Cey et al., 2006; Rosenbom et al., 2009). Discrete macropore models typically require fewer parameters (e.g., no mass transfer terms), use fewer assumptions than dual-continuum models, and thus can better represent the important physical macropore features and small-scale flow processes (Samardzioska and Popov, 2005). The main drawbacks of the discrete macropore approach are the lack of data on the macropore network and the associated difficulties with numerical discretization and computational resources. As a result, discrete macropore models are best suited to smaller scale applications where sufficient information on macropore geometry is available, such as the field infiltration experiments simulated in this study.

The description of flow in the macropore domain is particularly important when simulating macroporous systems because it can have a major influence on water and solute fluxes during preferential flow (Cey et al., 2006; Köhne et al., 2009). Again, various conceptualizations exist, but flow within macropores is often characterized by one of the following approaches: (i) capillary and gravity-driven flow described by Richards equation (e.g., Gerke and van Genuchten, 1993); (ii) gravity-driven flow described in a manner such as the kinematic wave equation (e.g., Larsbo et al., 2005); or (iii) a generalized capacity or “tipping bucket” approach, where downward flow occurs only after the soil reaches a specified saturation level such as field capacity (e.g., Weiler, 2005). It is not yet clear which of these approaches (if any) is suitable for describing macropore flow phenomena or what the parameterization of macropore hydraulics suggests about the underlying dynamics of unsaturated flow.

A difficulty with the application of virtually any unsaturated preferential flow model is parameterization (Šimůnek et al., 2003). A large number of parameters is required, some of which are difficult, if not impossible, to physically measure (e.g., macropore hydraulic properties, mass transfer terms, and macropore geometry). The large number of parameters can give rise to non-unique solutions, which has been well documented in the literature (Šimůnek and van Genuchten, 1996; Šimůnek et al., 2003; Köhne et al., 2009). It is therefore important to evaluate which parameters have the greatest influence on preferential flow behavior. This typically requires simplifying the natural system, often with column studies (Schwartz et al., 2000; Jarvis et al., 2007) or artificial macropores (Allaire-Leung et al., 2000; Akay and Fox, 2007), to isolate certain aspects of the flow system and obtain more detailed data. Transient data from tension infiltrometer (TI) experiments have been particularly useful in estimating model hydraulic parameters near saturation (Šimůnek and van Genuchten, 1996; Šimůnek et al., 2003) and more recently have been combined with Guelph permeameter data to derive macropore domain and mass transfer parameters for a dual-permeability model (Kodešová et al., 2010). In a similar way, Larsbo and Jarvis (2005) and Köhne et al. (2006) have shown that tracer transport data are useful for improving predictions of flow parameters. Köhne et al. (2009) have suggested that tracer methods can therefore be combined with other flow or hydraulic measurements to obtain better parameter estimates in field soils.

Recently, inverse methods of parameter estimation have been increasingly used to deal with issues of equifinality and parameter uncertainty. To date, inverse methods dealing with preferential flow in partially saturated systems primarily have been limited to one- and occasionally two-dimensional simulations (see review by Köhne et al., 2009). This can lead to the omission of critical topological information about the three-dimensional macropore network. The three-dimensional geometry of macropores is an important factor in preferential flow, yet little research has been conducted to evaluate its importance either under field conditions or within numerical simulations. Important considerations include the lateral distances between adjacent (but not connected) macropores as well as the vertical connectivity of individual macropores with the soil surface (i.e., buried macropores). The areal density of macropores controls macropore spacing, which can strongly influence water flow and solute transport in the vadose zone. A study by Weiler and Naef (2003) concluded that macropore density had a strong influence on the macropore drainage area and infiltration volumes. A series of laboratory and modeling studies (Allaire et al., 2002; Allaire-Leung et al., 2000) showed the importance of vertical and lateral macropore connectivity on water flow, tracer breakthrough, and matrix–macropore mass exchange. The influence of surface-connected and buried macropores was also evaluated using three-dimensional simulations in the context of hydraulic connections between macropores and a subsurface drain (Akay et al., 2008). Akay et al. (2008) demonstrated that hydraulic nonequilibrium (i.e., preferential flow within the macropore) did not occur as rapidly in buried macropores as in surface-connected macropores, primarily because flow from the soil surface to the top of the macropore was buffered by soil matrix material.

The main objectives of this study were to use a three-dimensional, variably saturated flow model for simulating preferential flow and transport in macroporous soils and evaluate the influence of hydraulic properties and macropore geometry in controlling flow and transport dynamics within macropores. The HYDRUS 2D/3D model was used to simulate water flow and dye tracer transport in small-scale TI experiments presented by Cey and Rudolph (2009). The three-dimensional macropore geometry, which has not been used in previous TI modeling studies, was explicitly included by utilizing a discrete macropore approach based on detailed field measurements of soil macropore characteristics. Multistep TI tests run at different tensions were used to calibrate the model and characterize flow in both the matrix and macropore domains. Sensitivity analyses were then conducted to elucidate the effect of various matrix and macropore properties on flow and transport behavior. Emphasis was placed on estimating and evaluating macropore hydraulic properties because of their nonlinear control of the flow system near saturation. The goals throughout the model calibration and sensitivity analyses were to: (i) calibrate macropore and matrix hydraulic properties to a range of infiltration conditions with and without macropore flow; (ii) evaluate the suitability of the matrix and macropore hydraulic parameter estimates and the commonly used van Genuchten water retention model (van Genuchten, 1980) for representing macropore flow; (iii) assess the influence of different macropore geometries, including varying macropore density, proximity, and vertical connectivity, on water flow and dye transport; and (iv) determine the relative value of different types of field measurements for assessing model performance.

## Methods

### Tension Infiltration Experiments

A series of field TI experiments were used as the basis to construct the numerical model and analyze the results (Cey and Rudolph, 2009). The infiltration experiments were conducted in July 2006 in an agricultural field near Walkerton, ON, Canada. The conservation-tilled field was cropped in a corn (*Zea mays* L.)–soybean [*Glycine max* (L.) Merr.] rotation and was planted with soybean at the time of the tests. The silt loam soil was classified as an Orthic Luvisol. The water table ranged from 1.1 m below grade in the spring to 1.5 m at the end of the growing season. Soil macropores were examined and mapped at the test site, revealing the macropores to be comprised of worm burrows, root holes, and fractures. Although fractures and root holes were present in the top 20 cm of the soil profile, the vast majority of macropores (>85%) were vertical worm burrows with diameters of about 0.1 to 1.0 cm. These worm burrows were the primary pathways for preferential transport of both solute and colloids (Cey et al., 2009). The areal density of cylindrical macropores ranged from 240 to 391 m^{−2} in the top 10 cm of the soil profile and increased to a peak of 783 m^{−2} at 40 cm depth (Cey and Rudolph, 2009). Approximately 18% of cylindrical macropores had diameters >0.5 cm, and the remainder were <0.5 cm in diameter.

The benefit of utilizing the TI for this research was the ability to control the degree of macropore flow by varying the tension on the TI disk between tests. Under different applied tensions, temporal (i.e., infiltration rates and soil moisture content) and spatial (i.e., end-run dye patterns) data sets were collected to characterize unsaturated flow behavior. Three infiltration tests (WK-D1, WK-D2, and WK-D3) were examined in the numerical analyses. Dye tracer was infiltrated into the soil using a TI with a 20.9-cm-diameter disk. The soil surfaces were prepared by removing surface vegetation and a thin layer of soil to obtain a level surface. A 0.3-cm-thick layer of glass bead contact material, recommended by Reynolds and Zebchuk (1996), was used to improve the hydraulic contact between the soil and the TI disk. Each infiltration test had two levels of pressure head applied: a short initial portion at a higher tension to wet the soil, followed by a longer portion at a lower tension during which the majority of the 1.0-L infiltration volume was applied. The second applied tension (maximum pressure head) largely controlled the extent of macropore flow. The maximum pressure heads, *h*_{max}, used were −5.2, −2.5, and −0.4 cm for Tests WK-D1, WK-D2, and WK-D3, respectively. Details for all TI tests are provided in Table 1 . Modeling focused on TI Test WK-D3 because it exhibited the greatest degree of flow and transport within macropores. For each test, cumulative infiltration (CI) and soil water content (SWC) directly beneath the TI disk were recorded automatically. The water content data represent averages from three time domain reflectometry probes that were inserted into the soil at a 45° angle from the perimeter of the TI disk to a depth of 7 cm.

Brilliant Blue FCF dye (Acid Blue 9, C.I. 42090) was used as a tracer during the infiltration tests to assess subsurface flow patterns. Brilliant Blue FCF was applied in the infiltration water at a concentration of 4000 mg L^{−1}, which was sufficient to produce visible dye stain patterns in all soil horizons. Following infiltration, the infiltration sites were excavated to determine the extent of dye-stained soil, take detailed photographs, map macropore and soil features, and collect soil samples. Image analysis was completed on the soil photographs to convert the observed dye patterns into three dye concentration classes (200–1500, 1500–3000, and >3000 mg L^{−1}) based on the intensity of dye staining (Cey and Rudolph, 2009).

Soil matrix properties were obtained from 10 intact soil cores collected from the A, B, and C horizons (Cey and Rudolph, 2009) and formed the basis for hydraulic properties used in the initial numerical simulations (Table 2 ). Soil water retention data were measured on the intact soil cores according to the pressure plate method described by Topp et al. (1993) at pressure heads of 0, −0.1, −0.2, −0.5, −1.0, and −3.0 m. Average water retention functions for each soil horizon were determined by least-squares fitting to the van Genuchten (1980) soil hydraulic function. The saturated hydraulic conductivity, *K*_{s}, was also measured on the same soil cores using the falling-head test procedure (Reynolds, 1993).

### Numerical Simulations

The numerical simulations were conducted using version 1.07 of the HYDRUS 2D/3D code, a finite element model capable of simulating three-dimensional water flow and solute transport in variably saturated porous media (Šimůnek et al., 2007). A three-dimensional, single-continuum, discrete macropore approach was implemented within HYDRUS to represent both the soil matrix and a discrete macropore domain, similar to that used by Akay et al. (2008) to examine laboratory column studies of macropore flow to tile drains.

### Three-Dimensional Discrete Macropore Model

The model domain was 100 by 50 cm horizontally and extended from the ground surface to the 100-cm depth to accommodate the maximum lateral and vertical extent of dye transport observed during field experiments (Fig. 1 ). The model domain represented half the infiltration site (taken through the centerline of the TI disk) to reduce simulation times and produce dye stain results comparable to the field data. As a result, the field measurements of CI that are presented have been adjusted accordingly (i.e., reduced by 50%). Domain discretization was designed to minimize water and solute mass balance errors. In the region directly beneath the TI disk, model grid spacing was reduced to a minimum of 0.5 cm to maintain stability and was gradually increased to a maximum of 5 cm away from the disk. This resulted in 64,119 modeled nodes and 354,816 elements.

The domain was subdivided into three soil horizons, as identified at the field site, with vertical macropores cutting across soil layers. The size of the discrete macropores included in the model was assumed to represent the mean macropore diameter observed in the field. This was based on field observations that identified macropores with mean diameters of approximately 0.5 cm as the dominant active macropore element. Furthermore, the macropore diameter was assumed to be constant for the entire vertical length of the macropore. Macropores were explicitly represented as vertically aligned soil pores with hydraulic properties representative of highly conductive macropore features. The areal density of macropores was consistent with field observations at the research site and was further evaluated during the sensitivity analyses. All macropore locations were randomly assigned in the region below and surrounding the TI to prevent bias in assigning locations. Due to limitations on time and computational resources, only one realization was simulated for a given macropore density. Restrictions on the generality of results were not analyzed, but are likely to have been small relative to other sources of uncertainty.

Flow boundaries were based on conditions during the TI experiments (Cey and Rudolph, 2009). Lateral boundaries were assigned no-flow conditions. A time-dependent specified head boundary was used for the infiltration surface. The head was specified at the base of the model to approximate the measured water table depth. The initial pressure head distribution within the model was based on the water table depth and SWC measured in the soil profile before the infiltration tests.

The glass bead contact layer used in the field was not included in the simulations; rather, the upper TI boundary was applied directly to the soil surface. This raises the question of whether the glass beads in the simulations might restrict preferential flow initiation due to a lack of macropore connectivity with the TI disk. Additional numerical simulations were conducted to test this hypothesis. The upper 1 cm of the model domain beneath the TI was assigned literature-derived hydraulic parameters (grain size range 106–250 μm, *K*_{s} = 0.66 cm min^{−1}) to represent the glass bead contact material (Reynolds and Zebchuk, 1996). The results of this assessment showed that the influence of the contact layer was minor, with <10% change in the simulated values for SWC, CI, or the depth of dye transport. The hydraulic characteristics of the glass bead layer prevented it from impeding preferential flow, and it was concluded that the simulations presented here (without a contact layer) were sufficient to describe the flow behavior observed in the field.

A specified mass flux was imposed on the infiltration area to represent the application of dye from the TI (input concentration = 4000 mg L^{−1}). The solute transport boundary condition at the base of the model was a third-type (or Cauchy) boundary. The initial solute concentration within the model was set to zero. Solute transport parameters for Brilliant Blue FCF dye were obtained from literature values. A free water diffusion coefficient of 3 × 10^{−4} cm^{2} min^{−1} was used (Hu, 2000). It is known that Brilliant Blue FCF will undergo sorption, with the clay fraction playing a major role (Ketelsen and Meyer-Windel, 1999). Thus, Langmuir-type isotherm coefficients were selected (maximum absorption *Q*_{max} = 8.6 g kg^{−1}; Langmuir equilibrium constant *K* = 4.7 L g^{−1}) from data by Ketelsen and Meyer-Windel (1999) for a soil with a similar grain size distribution, notably the clay content. Longitudinal and transverse dispersivity values in each soil horizon and macropore domain were assumed to be 10 and 1%, respectively, of the observed transport length, which is within the range of values reported by Gelhar et al. (1992) for porous media and more recently by Vanderborght and Vereecken (2007) for unsaturated soils.

The van Genuchten–Mualem (VG) model (van Genuchten, 1980) was used to describe the constitutive relations between capillary pressure, water content, and hydraulic conductivity for both matrix and macropore domains following previous work (e.g., Akay et al., 2008; Allaire et al., 2002; Rosenbom et al., 2009). The water retention relation for the VG model is given by[1]where *S*_{e}(*h*) = [θ(*h*) − θ_{r}]/(θ_{s} − θ_{r}) is the effective saturation, *h* is the pressure head, α, *n*, and *m* are empirical fitting parameters, θ is the volumetric water content, θ_{s} is the saturated water content, and θ_{r} is the residual water content. The hydraulic conductivity function, *K*(*h*), is defined as[2]where *l* is an empirical pore-connectivity coefficient. Hereafter, the subscript *f* will be used to denote the macropore domain and subscript *m* for the matrix domain. The VG α_{m} and *n*_{m} parameters assigned to the matrix were derived from least-squares fitting to laboratory measurements on intact samples from the field site (Table 2). Note that *K*_{s,m} was initially taken from laboratory falling-head tests but was later modified during the model calibration process (see below). The empirical parameter *l* was taken as 0.5 for porous media after van Genuchten (1980) and was assumed to be 1.0 for a vertically continuous macropore.

The hydraulic properties for the macropore domain, which were a particular focus of this study, were based initially on theoretical and literature values. The equivalent *K*_{s} of a smooth-walled capillary tube assumed to represent a macropore, *K*_{s,f}, is given by Poiseuille's equation (Bear, 1972):[3]where ρ_{w} is the density of water, *g* is acceleration due to gravity, *r* is the macropore radius, and μ is the viscosity of water. A theoretical *K*_{s,f} was calculated using a 0.5-cm-diameter macropore, but it was quickly noted that the simulation results were highly sensitive to this value and the theoretical value provided a poor fit to the observations. The value of *K*_{s,f} was subsequently reduced to better fit the observed data. Likewise, an estimate for the VG α_{f} parameter was obtained assuming that it represents the inverse of the absolute pore air-entry pressure (*h*_{ae}) given by[4]where σ is the interfacial tension between air and water and β is the contact angle measured through the wetting phase. A VG α_{f} value of 1.7 cm^{−1} was calculated for a 0.5-cm-diameter macropore. Theoretical estimates of the VG *n*_{f} parameter could not be obtained, so VG *n*_{f} values for a range of coarse-grained soil types were selected. It was found that model results were also sensitive to the VG *n*_{f} parameter so it became another key fitting parameter during model calibration and a focus for the sensitivity analysis.

The three-dimensional model was initially calibrated as follows to obtain a Base Case or starting point for the sensitivity analysis. First, the matrix hydraulic parameters were assigned from laboratory-derived values and macropore parameters were assigned from theoretical and literature values as noted above. Second, the *K*_{s,f} and VG *n*_{f} parameters were optimized based on fitting to TI Test WK-D3 (*h*_{max} = −0.4 cm), which exhibited the greatest degree of macropore flow. The laboratory-derived *K*_{s,m} value was also revised to improve the fit during this stage. Model results were assessed by comparing the simulation results to the field-measured CI, the average SWC below the TI disk, and the maximum depth of dye transport. The selected model input parameters were revised in an iterative process to arrive at an initial set of best-fit parameters. Third, *K*_{s,m} was optimized by comparing the simulated and observed results from Test WK-D1, which was conducted at a much lower infiltration pressure head (*h*_{max} = −5.2 cm) and had almost no evidence of preferential flow. Because Test WK-D1 produced flow primarily in the matrix domain of the A horizon and no information was available to calibrate the deeper soil horizons, the *K*_{s,m} value determined in this step was applied to all soil horizons. Finally, Test WK-D2, which had an intermediate pressure head (*h*_{max} = −2.5 cm) yet still had considerable macropore flow observed, was simulated to further refine the *K*_{s,f} and VG *n*_{f} values. The calibrated and refined model input parameters determined from this procedure are collectively referred to as the calibrated Base Case parameters and are given in Table 2. Inverse modeling was not conducted to optimize parameters because HYDRUS currently does not support this function for three-dimensional simulations. Furthermore, the goal of model calibration was not to completely optimize the input parameters but rather to arrive at a set of values that could reasonably represent the observed flow and transport behavior and serve as the basis for more detailed examination during the sensitivity analysis.

### Hydraulic Parameter Sensitivity Analyses

Using the calibrated hydraulic parameters from the Base Case, sensitivity analyses were conducted for TI Test WK-D3 to evaluate the model response to changes in the matrix and macropore hydraulic parameters and to assess the reliability of parameter estimates obtained from the model calibration. Key hydraulic parameters were adjusted upward and downward to varying degrees to gain insights into unsaturated flow and transport behavior. The matrix parameters *K*_{s,m} and VG α_{m} were modified upward and downward by up to one order of magnitude. The VG *n*_{m} parameter is typically less variable and has a practical lower limit of 1.0, hence it was only increased by up to a factor of 4 and decreased by up to 8% in the matrix. The different sensitivity cases for the soil matrix properties and the resulting parameter adjustments are listed in Table 3 . Similarly, *K*_{s,f} and VG α_{f} for the macropore domain were increased and decreased by up to one order of magnitude, while VG *n*_{f} was increased by as much as a factor of 3.3 and decreased by up to 30%. The adjusted model parameters for the macropore sensitivity runs are listed in Table 4 .

Parameter sensitivity was evaluated by comparing the simulation results to the field-measured values of SWC, CI, and the maximum depth of dye (as measured at the end of the infiltration period). The time-integrated transient responses for SWC and CI were evaluated by calculating the RMSE between the measured and simulated results for the duration of the infiltration period. The RMSE was calculated as(5)where *x*_{m} is the measured parameter value, *x*_{s} is the simulated parameter value, and *n* is the number of paired parameter values. The results were normalized by dividing the RMSE for each sensitivity run by the RMSE for the Base Case.

### Macropore Connectivity and Density Analyses

Additional simulations were also completed to test the influence of macropore connectivity and macropore density. To examine the presence of buried macropores, a situation identical to the Base Case was simulated except that all macropores were closed from the soil surface to a depth of 1 cm below grade. In this case, soil matrix properties were used for the upper 1-cm portion of the macropore and macropore hydraulic properties were used at depths >1 cm. A final set of sensitivity runs was conducted with the areal macropore density doubled and halved in comparison to the field-measured values (denoted in Table 4 as Cases 4a and 4b). In the case of macropore doubling (Case 4a), the same macropore configuration as the Base Case model was used and then an equal number of additional macropore locations were randomly assigned. The assigned macropore properties were identical to those from the Base Case and all macropores extended vertically through the entire model domain. To reduce the number of macropores in half (Case 4b), one of every two macropores from the Base Case model was effectively removed by assigning those cells the properties of the soil matrix. The macropore locations to be removed were also selected at random.

## Results and Discussion

### Three-Dimensional Model Calibration

Using initial laboratory, literature, and theoretically derived model parameters, Test WK-D3 was simulated to assess model performance. These results showed that the model could replicate the complex flow patterns observed in the field, especially in terms of vertical transport along macropores and mass exchange with the surrounding matrix, yet initial estimates for certain hydraulic parameters, most notably *K*_{s,f}, were clearly not suitable. The model overestimated the degree of preferential flow and transport and hence the selected hydraulic parameter estimates (e.g., *K*_{s,f}, VG *n*_{f}, and *K*_{s,m}) were revised to improve the model fit.

Test WK-D1, which primarily had transport in the matrix, was then simulated to calibrate *K*_{s,m}. The value of *K*_{s,m} was reduced to 0.01 cm min^{−1} from an initial laboratory-derived value of 0.54 cm min^{−1} (A horizon) to better match the field measurements. Reduction of *K*_{s,m} was not unexpected because the laboratory methodology (falling-head permeameter) was conducted on intact samples that included matrix and macropore components. The model fit to Test WK-D1 provided confidence that the matrix hydraulic parameters were well defined.

Final calibration to TI Test WK-D2 resulted in further refinement of *K*_{s,f} and VG *n*_{f}. Test WK-D2 had considerable flow in both matrix and macropore domains, yielding macropore parameters applicable to a wider range of infiltration conditions. The calibration procedure was admittedly not unique but took advantage of measured data across a range of flow conditions and provided a firm basis for heuristic investigation. The final set of calibrated input parameters are referred to as the Base Case and were used as the basis for further simulations and the sensitivity analysis.

The Base Case input parameters given in Table 2 were used to simulate Test WK-D3 (Fig. 2 ). Generally, the simulated results agreed closely with field observations of SWC, CI, and the dye distributions. The modeled SWC slightly underestimated the magnitude of SWC changes (Fig. 2a), whereas CI was often above the field measurements (Fig. 2b). Nevertheless, the timing of both responses correlated closely with the field data and captured the increase in the infiltration rate when the applied tension was changed at time *t* = 15 min. The RMSE of the SWC and CI data were 2.0 × 10^{−2} (unitless) and 0.0599 L, respectively. Of note is the unusual concave-upward shape of the CI curve, which could not be matched by the simulations. An increase in the infiltration rate was seen in other TI tests, particularly tests with larger applied pressure heads and increased macropore flow (see, for example, Test WK-D2 in Cey et al., 2009), but the reason for the unusual curves is unclear. It may have been a result of increasing macropore density with depth (peak density at the 30-cm depth) that was not incorporated in the simulations. It could also be a temperature effect that was not accounted for. In the field, the TI was shaded to prevent solar heating, but unavoidable temperature increases during the course of the day could have reduced fluid viscosity and increased *K*_{s}.

The simulated dye patterns were also consistent with field observations. The model results agreed with the maximum depth of dye transport (70 cm) and reproduced the characteristic dye-stained halos that were observed around macropores by Cey et al. (2009), as shown in Fig. 3 . This suggests that the model captured the relevant mass transfer processes between the matrix and macropore domains that are an important control on water flow and solute transport (Šimůnek et al., 2003).

The multistep parameter estimation approach used here is similar to that used by Kodešová et al. (2010), who used a two-dimensional axisymmetric model to estimate dual-permeability parameters by numerical inversion of TI and Guelph permeameter measurements. Kodešová et al. (2010) first determined matrix parameters using TI results, then optimized the macropore and mass exchange parameters using Guelph permeameter data. Our parameter estimation method differs in that we used a three-dimensional discrete macropore model to explicitly represent the macropore geometry (no mass exchange and macropore shape-factor terms), and multiple TI tests at different tensions were used to characterize both the matrix and macropore parameters. The three-dimensional discrete macropore model shows promise for estimating macropore hydraulic parameters, provided that detailed information on the macropore geometry is known.

As mentioned above, it became obvious early in the calibration procedure that the *K*_{s,f} estimates obtained using Poiseuille's equation were too large to properly describe the infiltration rates and dye transport observed in the field. The extent of macropore flow and transport were significantly overestimated using the theoretical *K*_{s,f} value of 5 × 10^{4} cm min^{−1} calculated for a 0.5-cm-diameter capillary tube. The calibration exercise yielded a *K*_{s,f} estimate on the order of 10^{2} cm min^{−1}, which equates to a macropore diameter of 0.03 cm, or approximately one order of magnitude smaller than observed in the field. Potential explanations for the observed discrepancy between the theoretical and modeled results are: (i) film, rivulet, or droplet flow was potentially present within the macropores, (ii) there were localized macropore constrictions or discontinuities, or (iii) the VG model is not suitable for describing macropore hydraulic properties. Liquid film flow within fracture or macropore systems has been documented in previous literature, including Tokunaga and Wan (1997), Allaire-Leung et al. (2000), and Cey and Rudolph (2009), and would effectively decrease macropore relative *K* values. When flow rates into a macropore are less than the full capacity of the macropore, instability can develop that results in flow as droplets or rivulets (Dragila and Weisbrod, 2003). Even for steady inflow conditions, flow within the macropore in these situations is often intermittent or pulsating (Su et al., 1999; Tofteng et al., 2002). These complex flows cannot be described by Darcy's equation and are probably not in hydraulic equilibrium with the surrounding matrix. Similarly, constrictions or discontinuities along macropores could reduce the maximum potential flow rate through a macropore to a fraction of that predicted for a full capillary tube. Advances are being made in characterizing the structure of macropores and macropore networks, which would help address this issue, but these studies are primarily limited to pore-scale analysis (Perret et al., 1999; Sander et al., 2008).

An overarching uncertainty is the applicability of the VG model to macroporous soils, which is closely related to the complex and intermittent flows described above. Even though the VG model has commonly been applied to the macropore domain in unsaturated flow simulations, no studies have clearly shown the applicability of the VG model to describe soil–water relationships in macropores. As an alternative, the physically based liquid configuration model of Or and Tuller (Or and Tuller, 1999; Tuller and Or, 2002) shows promise for describing unsaturated hydraulic conductivity in structured soils and compares reasonably well with the more empirical VG model. Both these approaches rely on describing flow across a range of pore sizes, whereas the approach in this study was to simulate flow in macropores of a single size. Because the Or and Tuller (Tuller and Or, 2002) model has the capability of accounting for film flow in irregularly shaped macropores, it could be used to develop a soil–water relationship appropriate for a single macropore or narrow range of macropore sizes. The Or and Tuller model, however, is limited to laminar flow in relatively small soil pores and cannot account for pulsating non-Darcian flows that have been observed when fluids invade dry macropores. The Richards equation may be an invalid representation for macropore flow processes. New models for representing unsaturated flow in macroporous systems are needed along with carefully measured data sets for model testing.

### Sensitivity to Hydraulic Parameters

Sensitivity of the calibrated three-dimensional model was analyzed by varying the hydraulic properties for both matrix (Table 3) and macropores (Table 4). The results showing the influence of hydraulic properties on transient SWC and CI responses are presented in Fig. 4 . Hydraulic parameters were increased or decreased by a given multiplication factor, which is plotted on the *x* axis. The resulting RMSE values relative to the Base Case are plotted on the *y* axis, such that RMSE ratios <1 signify a better fit to the observed data than the Base Case and ratios >1 signify a poorer fit.

The results in Fig. 4 support the selection of the calibrated parameter values. Adjusting any one of the hydraulic parameters can improve the fit to a particular measured data set but at the expense of a poorer fit overall. For example, increasing VG α_{m} or reducing *K*_{s,m} and VG *n*_{m} can reduce the RMSE for CI by as much as 50% (Fig. 4b); however, this would result in a larger RMSE for the SWC in all instances. The transient SWC response was generally not highly sensitive to changes in macropore parameters (Fig. 4a). Although small improvements in model performance (<10%) were seen by increasing the VG α_{f} values and decreasing the *K*_{s,f} and VG *n*_{f} values for the macropores, these parameters did not provide an improved match to the CI data (Fig. 4b and 5c) or the depth of dye transport (Fig. 5a). None of the modified macropore hydraulic parameters used in the sensitivity simulations reduced the RMSE for CI relative to the Base Case (Fig. 4b). The largest RMSE values were associated with the CI response to changes in macropore parameters, indicating that CI was a useful measure for refining the model parameters, especially the difficult-to-measure macropore hydraulic parameters.

The influence of hydraulic properties on end-of-infiltration values is plotted in Fig. 5 . Note that the CI values shown are half the actual measured values because the simulations were conducted using only half the TI disk. There are two key findings in these data. First, the system response was generally more sensitive to macropore parameters than matrix parameters. Changing the macropore hydraulic properties had precisely the opposite effect on dye transport as changes to matrix properties, although the magnitude of the response was considerably greater in the case of the macropores (see Fig. 5a). On the other hand, the influence of the hydraulic properties on the total CI was very similar for both matrix and macropores (Fig. 5c), suggesting a similar dynamic response to changes in hydraulic properties, but the magnitude of the response was again much larger for the macropores. For instance, an increase of *K*_{s} and VG *n* values resulted in an increase in the simulated CI in both the matrix and macropore cases. Increasing *K*_{s,m} by a factor of 10 times, however, resulted in a modest increase in the simulated infiltration volume to 0.83 L, whereas the same increase in *K*_{s,f} resulted in a simulated CI of >4.0 L. Assuming that our conceptual approach is valid, the results imply that macropore parameters can be more readily constrained than matrix parameters given a suitable measured data set.

The second key finding is that the hydraulic parameters had more influence on the depth of dye transport (Fig. 5a) than on the final SWC or CI (Fig. 5b and 5c). Notice that the data in Fig. 5a are on a logarithmic scale that covers four orders of magnitude. The depth of dye transport is a critical parameter because it reflects the likelihood of contaminant breakthrough to a receptor, such as a tile drain or underlying aquifer. Studies by Köhne et al. (2006) and Larsbo and Jarvis (2005) have shown that incorporation of flow and transport data together significantly improves hydraulic parameter estimation in soils with preferential flow. The fact that dye depth was most sensitive to changes in the macropore properties is consistent with field observations of extensive preferential flow within vertical macropores and the recognition that dye transport was controlled by macropore flow. Increasing VG *n*_{f} and *K*_{s,f} by factors of 3.3 and 10 times, respectively, resulted in a simulated dye transport distance of about 1400 cm (or 20 times the Base Case). The sensitivity of dye depth to the macropore hydraulic parameters indicates that the solute transport distance was the most useful measure for refining these parameters. Our results suggest that solute transport measurements should be included in field studies and associated simulations to better constrain hydraulic parameters in macroporous soils.

An instructive pattern appears in the comparison of Fig. 5a and 5c. For matrix parameters, the sensitivity plots for the depth of dye and final CI are nearly mirror images of each other about a horizontal plane. In other words, increasing *K*_{s,m} increases CI, but the depth of dye transport decreases (and vice versa). The same holds true for the other matrix hydraulic parameters considered during the sensitivity analyses and suggests a strong negative correlation of the matrix parameters. These findings point to non-unique model parameterization even with the use of solute transport information. The macropore hydraulic parameters do not behave in the same manner. In fact, the sensitivity curves for dye depth and CI show a strikingly similar pattern when the macropore hydraulic parameters are varied (Fig. 5a and 5c), indicating fewer issues with parameter uniqueness for the macropores.

Field measurements of dye transport patterns and infiltration volume (rate) were most useful for estimating the hydraulic properties of the macropores. Data on infiltration volumes and dye depth were necessary for the model to capture the macropore-dominated flow behavior observed in the field. Water or solute flux data in the subsurface would also have been beneficial for macropore parameterization, as has been demonstrated in previous studies of preferential flow to tile drains or drainage lysimeters (e.g., Akay et al., 2008; Gerke and Köhne, 2004; Kung et al., 2005). On the other hand, SWC data had the least value for constraining model input parameters in the current study. Because macroporosity is a small fraction of the total porosity (<1%), SWC measurements had little influence on macropore parameters. Temporal SWC data to monitor the wetting front advance was beneficial for refining the matrix hydraulic properties, but the final SWC values were insensitive to matrix and macropore parameters (Fig. 5b) because the soil beneath the TI disk was at or near saturation at the end of the infiltration experiment.

In spite of the fit of the simulation results to field observations and the appearance of well-constrained hydraulic parameters, several macropore parameter values fell well outside the expected range. In addition to the *K*_{s,f} values discussed above, the final VG *n*_{f} value was also much smaller than expected based on extension of the values for coarse-textured media. The calibration and sensitivity results, however, showed that larger VG *n*_{f} values would not have improved the model results. As discussed above, this result may also stem from non-Darcian flow within macropores. Treating macropores using standard methods for porous media (i.e., Darcy's law, the Richards equation, and VG soil hydraulic functions) clearly has limitations for describing complex flows in macropores.

### Influence of Macropore Connectivity

Field and laboratory studies have shown that the number and size of surface macropores alone are not sufficient to predict water and solute transport through structured soils but that macropore length and continuity can also dictate the extent of preferential flow (Akay and Fox, 2007; Allaire-Leung et al., 2000; Logsdon, 1995; Munyankusi et al., 1994). Macropore connectivity was evaluated in this study by disconnecting the simulated macropores from the ground surface. The model was run using macropores that were closed at the surface (simulated as soil matrix) but became open at a depth of 1 cm below grade. All other model domain and input parameters were the same as the calibrated Base Case.

As shown in Fig. 6, the 1-cm disconnection at the top of the macropores had a major effect on the simulation results. The model was able to match early infiltration into the soil matrix during the initial applied tension (up to *t* = 15 min) but could not reproduce the SWC or CI responses at later times when macropore flow was prevalent (Fig. 6a and 6b ). The total CI was reduced by >50% compared with the field-measured values or the Base Case simulation that included surface-connected macropores. Solute transport was reduced, with dye reaching depths of <10 cm at all points in the model (Fig. 6c). While initiation of macropore flow can be observed at the end of the simulation, all infiltration occurred via the matrix until sufficient pore pressures developed in the subsurface to initiate flow within the buried macropores. Macropore connectivity was difficult to evaluate in the field (Cey and Rudolph, 2009), but the degree of vertical dye transport and observations of worm burrows near surface suggested that some fraction of the macropores were connected with the soil surface; however, not all macropores were open to the surface.

The model results suggest that a small discontinuity may be sufficient to significantly reduce the likelihood of preferential flow and transport. As described in Cey and Rudolph (2009), “bulbs” of dye-stained soil developed near the base of vertical macropores, indicating that the lack of macropore continuity restricted preferential flow and resulted in increased mass exchange with the soil matrix. These observations are entirely consistent with the simulation results and indicate a need to further examine more realistic macropore networks, with macropores opening and terminating at different depths. The numerical results from this study also agree with other studies (Allaire-Leung et al., 2000; Akay and Fox, 2007; Akay et al., 2008) that have identified reduced macropore flow as a result of buried macropores.

Macropore discontinuities may also explain the difference between the theoretical and fitted *K*_{s,f} values. The absence of surface-connected macropores significantly reduced the vertical migration of water and solute in the simulations. An increase in *K*_{s,f} would counteract this in part, as identified in the sensitivity analyses (Fig. 5a and 5c), and bring the value closer to the theoretical prediction. Nearly all preferential flow and transport simulations are based on the use of a fully connected macropore network (Gerke and van Genuchten, 1993; Larsbo et al., 2005; Köhne et al., 2006). A small disconnection along a macropore could alter flow, as shown here, and thereby significantly change the parameter values obtained from such models. Macropore connectivity should thus be a high-priority area of research, with more work needed to evaluate macropore connectivity in field settings and properly incorporate this information into preferential flow models.

### Influence of Macropore Density

The total number of macropores mapped in the TI experiments was considerably greater than the number of active dye-stained macropores (Cey and Rudolph, 2009). Given the importance of macropore flow processes, the model was used to examine the influence of macropore density on flow and transport. Based on field observations of active macropores, the Base Case model included four macropores below the half TI disk. Macropore density was varied by doubling the number of macropores to eight below the TI disk, as well as halving it to two under the disk. The results of these simulations are shown in Fig. 7 . Changing the macropore density had relatively little influence on the SWC (Fig. 7a), although the time-integrated SWC response was improved in comparison to the Base Case (Fig. 2a). The influence was much more pronounced on the CI response and the depth of dye transport. Infiltration volumes appeared to be almost linearly related to the number of macropores (Fig. 7b). Dye transport behavior, in contrast, was highly nonlinear. Doubling the macropore density caused the dye to reach the bottom of the model domain (1.0 m) within the simulation period, but a reduction to two macropores resulted in the simulated dye reaching only the 60-cm depth (Fig. 7c).

Because macropore density exerts a strong influence on the overall model response, the calibrated simulation parameters are a function of macropore density. Although it is possible to recalibrate the model hydraulic parameters to account for different macropore densities, this is not a trivial task. Dual-permeability approaches have an advantage in this regard because the model input parameters can be modified to account for varying macropore density (or spacing) without the need to refine the model domain or reassign properties within discrete macropores.

### Influence of Macropore Spacing

Changes in the flow volume and transport depth are due not only to changes in the number of macropores but also to the changing proximity of macropores as the density increases or decreases. Using the results from Case 4b (see Table 4; Fig. 7), we can compare the solute transport results for macropores with different spacings. Figure 8 shows the result of dye transport on a vertical section through two macropores: one situated within about 1.5 cm of another active macropore (located behind the section) and the other a more distal macropore that was not closely neighboring another active macropore (about 4 cm from the nearest macropore). The model results demonstrated that vertical dye transport was approximately 40% deeper for the macropores that were in closer proximity. This is due in large part to the significant macropore–matrix interaction that occurred in these systems. Flow in closely neighboring macropores led to pressure buildup in the soil matrix and interconnection of wetting fronts surrounding the macropores, which subsequently enhanced preferential flow within the group of macropores.

Our simulation results compare well with field observations of increased macropore flow and deeper dye transport around clusters of macropores. Cey and Rudolph (2009) demonstrated that the presence of neighboring macropores enhanced vertical dye transport, and they observed significant lateral flow between adjacent macropores. Allaire-Leung et al. (2000) concluded that interaction between neighboring macropores can enhance the vertical transport of solutes. They further concluded that buried macropores can become important conduits for subsurface flow and transport in the presence of closely spaced macropores. The influence of neighboring macropores is not unlike the merger of gravity-driven flow fingers (Glass and Nicholl, 1996) or the occurrence of funneled flow (Kung, 1990), both of which result in localized fluid concentration and increased vertical flow velocities. The proximity of macropores can be accounted for directly in discrete macropore simulations as shown here. The more common dual-permeability approach, however, generally assumes a uniform distribution of macropores, potentially underestimating the degree of preferential flow in regions with increased macropore density or proximity.

## Summary and Conclusions

Tension infiltration experiments previously described by Cey and Rudolph (2009) were simulated to gain insights into preferential flow behavior within macroporous soil under near-saturated conditions. The explicit representation of a simplified three-dimensional macropore geometry in the simulations provided an alternative to more frequently used dual-permeability models for examining macropore properties and flow processes. The three-dimensional discrete macropore model was iteratively calibrated to a series of three TI tests, each run at different tensions, yielding manually optimized matrix and macropore hydraulic parameters. The calibrated Base Case simulation of Test WK-D3 appeared to capture the relevant components of flow and transport during infiltration, including the temporal response to changes in the infiltration rate, total depths of dye transport, and the significant matrix–macropore interaction that occurred. Subsequent sensitivity analyses of the matrix and macropore hydraulic properties revealed that dye transport was controlled primarily by the macropore parameters, thereby constraining macropore hydraulic properties that are difficult to estimate independently. Field measurements of CI and (especially) the depth of dye staining were most useful for refining the macropore hydraulic parameters. Monitoring solute transport depths and breakthroughs is recommended, in addition to more traditional flow and water content data, to improve parameter estimates in field studies of preferential flow.

Although sensitivity analysis confirmed the suitability of the hydraulic parameter estimates in terms of matching the field observations, questions remain about the applicability of the results. The model produced unexpected estimates for several hydraulic parameters. The calibrated three-dimensional model had *K*_{s,m} more than an order of magnitude lower than laboratory estimates, suggesting that laboratory estimates of *K*_{s,m} may not properly represent matrix properties unless measures are taken to exclude macropore flow. More significantly, *K*_{s,f} and VG *n*_{f} determined from the calibrated model were much lower than theoretical estimates would dictate, indicating limitations in the use of conventional porous medium representations for flow in macropores. These findings indicate that a Darcy-based model may not adequately represent unsaturated flow in macroporous systems. Research is needed to develop physically based models of the hydraulic functioning of macropores near saturation. Careful measurements of flow characteristics in larger gravitational pores are also needed to aid in the development and testing of these models.

The number of surface-connected macropores played an important role in the determination of *K*_{s,f}. Simulation results showed that a small (1-cm) break in vertical macropore connectivity at surface severely restricted flow and transport along macropores at greater depth. Similarly, increasing the macropore density and proximity strengthened the lateral connections between different macropores, producing greater infiltration volumes and depths of solute transport. The strong link between the macropore parameter estimates and the dye transport results implies that macropore geometry and connectivity appreciably influence *K*_{s,f} (and other parameter) estimates. Currently, the most common preferential flow models do not account for vertical macropore connectivity and variable macropore spacing, which may lead to concerns about the accuracy of macropore hydraulic parameter and mass exchange parameter estimates, particularly at larger model scales or during more extreme preferential flow events. Finally, the results suggest that land use practices that can alter the number and connectivity of macropores will have major implications for preferential flow in macroporous soils. More research is needed to improve measurements of macropore geometry and connectivity in the field and incorporate that information into flow and transport models.

## Acknowledgments

We gratefully acknowledge financial support from the University of Calgary. We also thank Anneli Schöniger for assistance with the initial numerical simulations.

## Footnotes

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.

- Received September 3, 2010.