Performance of preferential flow models in predicting infiltration through a remolded soil with artificial macropores

Preferential flow in soils contributes to water and contaminant leaching to groundwater and subsurface drains. At practical scales, preferential flow is often neglected due to the lack of consensus or guidance on how to simulate this process. Dual‐permeability (DP) approaches are a standard, but their applicability in practice is limited by the large number of parameters involved. The source‐responsive (SR) model simulates preferential flow as flow films along the walls of macropores that interact with the soil matrix. The objective of this study was to evaluate the ability of the SR model and a DP model with calibration limited to macropore parameters to simulate hydrologic variables of interest for practical applications. Infiltration experiments were conducted in a 150‐cm‐long by 50‐cm‐wide by 40‐cm‐deep soil column with vertical artificial macropores, with preferential flow being the most prevalent flow process. Two macropore configurations and various rainfall rates were tested. The DP model provided relatively accurate predictions of the bottom drainage rates, but inaccurate soil parameter estimates resulted in poor predictions of matrix flow. Critical limitations in the SR model formulation were highlighted, and the model was revised. The refined SR model provided accurate predictions of bottom drainage rates and water contents under steady‐state rainfall. The DP model with a reduced calibration appeared valid to estimate subsurface fluxes with negligible matrix flow. The SR model applicability is limited to steady‐state scenarios with negligible matrix flow. Critical developments are needed for the SR model to be applicable in practice.

and hillslope scales is critical to improving water pollution management (Djabelkhir, Lauvernet, Kraft, & Carluer, 2017;Weiler, 2005). In particular, riparian buffers implemented to reduce contaminant loading from runoff to surface waters are uniquely susceptible to preferential flow because of their geomorphic configuration, abundant root networks and other biological activity, and frequent drying-wetting cycles (Fox, 2019;Orozco-López, Muñoz-Carpena, Gao, & Fox, 2018). Also, their proximity to water bodies often implies the presence of a shallow water table, which is reached more easily in case of preferential flow (Fox, Muñoz-Carpena, & Purvis, 2018).
Quasiphysical modeling approaches based on simplified representations for macropore flow have been developed as an alternative to simulate preferential water flow processes at larger scales relevant to management (Gerke, 2006;Jarvis, 2007;Köhne, Köhne, & Šimůnek, 2009;Šimůnek, Jarvis, van Genuchten, & Gärdenäs, 2003). Most of them are one-dimensional and can be classified into two types. Single-domain approaches are typically based on the Richards equation for variably saturated flow and use modified water retention curves to account for preferential flow (Ross & Smettem, 2000). Dual-domain approaches divide the porous medium into two domains corresponding to the soil matrix and the macropore system where water flows at distinct rates. The domains are not geometrically separated but are conceptually superimposed over the same soil volume (Köhne et al., 2009). Among dualdomain approaches, mobile-immobile approaches assume that water in the matrix domain is immobile and that water flows in the macropore domain only. Dual-porosity approaches are similar to mobile-immobile approaches but account for water exchanges between the matrix and macropore domains. Dual-permeability (DP) approaches are standard approaches that simulate water flow in both domains, as well as exchanges. Flow in the matrix domain is assumed to be controlled by capillary and gravity forces and is modeled by the Richards equation. Different approaches have been proposed to simulate macropore flow in DP models. Some models assume that macropore flow is governed by capillary and gravity as well and

Core Ideas
• A numerical study of water flow in a remolded soil with artificial macropores was performed. • Source-responsive and dual-permeability models were compared and evaluated. • Dual permeability model with reduced calibration predicts bottom drainage rates. • An empirical, refined source-responsive model performs well for simplified scenarios. • A source-responsive model cannot be applied to dynamic rainfall scenarios.
use a second Richards equation (Gerke & van Genuchten, 1993), whereas other models assume that macropore flow is dominated by gravity effects and use a kinematic wave equation (Beven & Germann, 1981;Jarvis, 1994). Dualporosity approaches based on the kinematic wave equation require slightly fewer parameters, which can be seen as an advantage. Several studies have been conducted to assess the performance of DP models at different scales. Köhne and Mohanty (2005) evaluated a DP model's ability to simulate intradomain and interdomain water flow processes measured in a sand column with an artificial macropore. The DP model was also compared with the three-dimensional axisymmetric Richards equation used as a reference model. The DP model was almost equally capable of simulating intradomain and interdomain water flows as the reference model. Gärdenäs, Šimůnek, Jarvis, and van Genuchten (2006) compared the performance of a single-domain model using modified hydraulic properties near saturation, a mobile-immobile model, and a DP model based on two Richards equations in simulating drainage and pesticide concentrations measured at a tiledrained field. Their results showed that the DP approach most accurately simulated preferential drainage flow and that the single-domain and mobile-immobile approaches failed to capture this process. Other studies confirmed that DP models provide suitable results when calibrated for a specific scenario at the field and catchment scale (Christiansen, Thorsen, Clausen, Hansen, & Refsgaard, 2004;Larsson & Jarvis, 1999). However, the large number of parameters involved in DP models has impeded their application to practical scales because most parameters are difficult or impossible to quantify due to the spatial heterogeneity in vadose zone properties (Larsbo et al., 2005). Also, the discrepancy between the spatial and temporal scales of observations and models limits the application of inverse calibration techniques (Weiler, 2005).

F I G U R E 1 Experimental setup:
(1) topsoil, (2) gravel layer, (3) rainfall simulator, (4) drainage collector, (5) runoff collector, and (6) soil moisture sensors Another approach proposed by Nimmo (2010) is to simulate the macropore domain as a source-responsive (SR) domain that is not always active and in which flow is conceptualized as free-surface films along the walls of macropores. The concept of film-flow used in the SR model has been recently put forward as a more realistic way of simulating preferential flow (Germann, 2020). In this approach, model parameters are representative of material properties and can easily be interpreted physically. Only one formulation that describes the SR model has been proposed to date. It is based on simplified conditions: constant rainfall rates and negligible contribution of surface infiltration in the matrix domain (Nimmo & Mitchell, 2013). This version of the SR model has only been evaluated based on a couple of laboratory experiments and point-scale measurements of water content without consideration for water fluxes (Nimmo & Mitchell, 2013).
The objective of this research was to evaluate the performance and applicability of the SR model for practical scales in comparison with a DP model with calibration focused on the macropore domain parameters. Infiltration experiments were conducted in a 150-cm-long, 50-cmwide, and 40-cm-deep soil column with artificial macropores to obtain comprehensive, detailed data at a scale more representative of practical scales. "Plot" boxes that are a few meters in length and width and more or less half a meter in depth have been used to simulate infiltration and runoff processes representative of field plots (Kleinman, Sharpley, Veith, Maguire, & Vadas, 2004) and vegetative filter strips . Such "mesoscale" experiments allow greater control of confounding variables including soil properties and heterogeneity and comprehensive monitoring at a scale more representative of field conditions.

INFILTRATION EXPERIMENTS
The experimental setup consisted of a 150-cm long, 50-cmwide, and 40-cm-deep box with a 2% slope ( Figure 1). Percolated water was captured through a bottom drain located near the downstream end of the box. Surface runoff was measured through a weir with the invert located 34 cm above the bottom. A 2-cm layer of gravel was placed at the bottom of the box to facilitate drainage, and a 32-cm layer of sandy-loam topsoil screened to remove larger clods of soil, stones, and roots (74.2% sand, 16.3% silt, and 9.6% clay) was packed in two lifts with a dry bulk density of 1.6 g cm −3 . The surface of the soil was leveled to prevent ponding at the surface. Similar to the setup presented by Kang, Amoozegar, Heitman, and McLaughlin (2014), a rainfall simulator consisting of a single nozzle (FullJet, 1/4GG-12W, Spraying System Company) mounted on a tripod was placed above the box and used to apply uniform and constant rainfall to the soil surface. Three rain gauges placed around the box were used to check the uniformity of rainfall and measure the total rainfall depth at the end of a test. Two ECRN-100 tipping rain gauges connected to EM-50 data loggers (Meter) were used to collect and record surface runoff and bottom drainage. In addition, four rows of soil moisture sensors EC-5 connected to EM-50 data loggers (Meter) were placed in the soil 115, 123, 131, and 139 cm from the upstream end. Four sensors were used for each row and placed 4, 11, 18, and 25 cm below the soil surface. Runoff, bottom drainage discharge, and soil water content were recorded at 60-s intervals. A plain steel round rod with a diameter of 4.76 mm was used to create vertical macropores through the soil. Two macropore configurations with equally spaced vertical macropores going all the way from the surface to the bottom were simulated: Configuration 1 had 39 macropores and Configuration 2 had 59 macropores. The soil was packed at a high soil density, ensuring the relative stability of the macropores. Before each test, macropores were reformed using the rod. It was not possible to control macropore geometry during a test, but no significant blockage was observed. However, a slight increase in the macropore diameter was noticed as the set of experiments progressed. For each macropore configuration, several tests with rainfall rates ranging from 1.9 to 8.5 cm h −1 for Configuration 1 and 2.4 to 12.7 cm h −1 for Configuration 2 were run. Test durations ranged from 1 to 4 h. Data recording started ∼10 min before rainfall application and continued for at least 3 h after the rainfall ended. Runoff was only observed for high rainfall intensities. The higher the rainfall intensity, the sooner bottom drainage started (Table 1). Drainage started to decrease soon after the rainfall was ended. Observations suggested that the flow wave moving in the matrix domain did not reach the bottom of the box and that all the drainage observed was due to preferential flow effects. Sensors located at similar depths exhibited similar temporal dynamics but with large variability in values. Once the rainfall started, sensors showed a rapid increase in water contents that then plateaued, which was interpreted as the preferential flow wave.

Source-responsive model
The more complete formulation of the SR model proposed by Nimmo and Mitchell (2013) was selected for numerical modeling. In this formulation, the contribution of matrix flow to changes in soil water content was neglected. The rate of change of soil water content (θ) in a representative element volume (REV) in the vicinity of a macropore was defined as where M is the macropore internal facial area capable of supporting film flow per unit volume of the bulk medium [L −1 ], D is the effective hydraulic diffusivity of the matrix material [L 2 T −1 ], G is a dimensionless factor representing the effect of the geometrical configuration of macropores, θ e is the equilibrium water content where there are no exchanges between the matrix and macropores, and f is the degree of activation of source responsive transport estimated as where i s is the infiltration rate in the macropore domain at the surface [L T −1 ], i 0 is the maximum rate of SR infiltration [L T −1 ], and t l is the activation lag time [T]. Nimmo and Mitchell (2013) indicated that t l depended on the reciprocal of the film flow velocity and was proportional to −2∕3 s , expressed as where C l is a model parameter. Assuming that i s is constant, the evolution of the water content under constant macropore infiltration conditions could be calculated as where θ 0 is the initial water content. Because this formulation of the SR model neglects the contribution of matrix flow to changes in θ, the surface infiltration in the matrix was not calculated. Therefore, in this study, a simple approach was used to estimate the infiltration rate in the macropore domain at the surface as a fraction of the rainfall rate, I: where ω is a model parameter. According to Equation 4, when ≫ l , the water content θ equals θ e . It should be noted that Equation 4 does not predict changes in soil moisture after the rain has ended. Another important limitation of this formulation is the absence of a mass balance. Physically, exchanges from the macropore to the matrix domain at a specific time and location should be limited by the amount of water available, which is not the case in Equation 4, in which the rate of exchange is controlled by the difference between the equilibrium water content and initial water content in the REV. For this study, a mass balance component was added to the SR model. The downward flux density in the macropore domain q m was estimated from the rate of infiltration in the macropore domain at the surface and abstractions due to exchanges from the macropore to the matrix (Nimmo & Mitchell, 2013): A constraint was placed on the rate of exchange from the macropore to the matrix domain: such that q m remained positive or null at all times.

Dual-permeability model
In the DP model, flow in the matrix domain was simulated using Richards equation: where C(h) = ∂θ mi /∂h [L −1 ] is the water capacity, θ mi [L 3 L −3 ] is the volumetric water content of the matrix domain, h [L] is the soil water pressure head, t is time is the unsaturated hydraulic conductivity, and S [T −1 ] is a source term accounting for water transfer between the macropore and matrix domain and that is positive in the macropore-to-matrix direction. Soil hydraulic properties were described using the van Genuchten-Mualem model (van Genuchten, 1980): where θ r and θ s [L 3 L −3 ] are the residual and saturated water contents, respectively; α [L −1 ] is the inverse of the air-entry pressure head and n is a model parameter. Flow in the macropore domain was simulated using a kinematic wave equation: where θ ma [L 3 L −3 ] is the water content of the macropore domain and Q ma [L T −1 ] is the flow rate in the macropore domain. Similarly to Larsbo et al. (2005), the relationship between Q ma and θ ma is expressed as where S ma = θ ma /β [L 3 L −3 ] is the degree of saturation of the macropore domain, β [L 3 L −3 ] is the saturated water content of the macropore domain, and a and b are model parameters. This equation reduces to Q ma = a when the macropore domain is saturated; therefore, a corresponds to the maximum flux in the macropore domain.
The exchange term S [T −1 ] between the macropore and matrix domains was simulated based on saturation in each domain using a first-order rate equation (Šimůnek et al., 2003). Flow from the macropores into the matrix was only allowed when the saturation of the macropore domain S ma exceeded the saturation of the matrix S mi . Flow from the matrix domain to the macropore domain was only allowed when the matrix was saturated. { = κ ( ma − mi ) if ma > mi or mi = 1 = 0 otherwise (11) where κ [T −1 ] is a first-order rate coefficient.
The representative θ was calculated as θ m = θ mi + θ ma . The infiltration at the surface of the matrix domain depended on the rainfall rate and soil moisture of the surface layer (van Dam & Feddes, 2000). The potential infiltration rate i p in the matrix was estimated using Darcy's law: where K 1/2 [L T −1 ] is the average of the hydraulic conductivity between the two upper nodes, h s [L] represents the height of water ponding at the surface, h 1/2 [L] is the average pressure head between the two upper nodes, and Δz is the distance between the two upper nodes.
The actual rate of infiltration in the soil matrix, i ma , was taken as the minimum between i p and the rainfall on the matrix domain (1 -β)I, where I [L T −1 ] is the rate of rainfall application.
The surface inflow in the macropore domain i mi was calculated as the minimum between the excess rainfall that did not infiltrate in the matrix domain and the infiltration capacity of the macropore domain, corresponding to the parameter a from Equation 10. The bottom boundary condition was treated as free drainage for the macropore domain and as a seepage face for the matrix domain, caused by the interface with the gravel layer.

Simulation of infiltration experiments
In this study, only tests that did not produce surface runoff were simulated. The soil column was represented as a onedimensional vertical grid discretized with a uniform node spacing of 1 cm. Constant time steps of 60 and 5 s were used for the SR model and the DP model, respectively. For the SR model, the simulation duration was taken as the duration of the rainfall. As the DP model was capable of simulating flow recession, simulations were run for a duration equal to the duration of the rainfall event plus 1 h.
Soil properties and corresponding model parameters were assumed to be uniform in depth. The DP model had nine calibration parameters: five for the matrix domain (θ s , θ r , α, n, and K sat [saturated hydraulic conductivity]), three for the macropore domain (β, a, and b) and one for the exchange term (κ). The SR model had seven calibration parameters (ω, M, G, i 0 , C l , D, and θ e ). The number of parameters in the SR model is lower because matrix flow has been neglected. Because the contribution of preferential flow to the bottom drainage was significant, the DP model calibration focused on macropore domain parameters and soil matrix parameters were estimated based on the soil texture and bulk density using pedotransfer functions (Rosetta version 1.1; θ s = 0.38, θ r = 0.07, α = 0.04 cm −1 , n = 1.35). For the SR model, the soil diffusivity, D, and geometric factor for the macropore network, G, were set to reference values from Nimmo and Mitchell (2013): D = 10 −2 m 2 h −1 and G = 0.5.
The bottom drainage discharge, q b , was used as a calibration target because it represented a key model output for application and design-focused simulations. The θ is also a variable of interest for applications such as the design or assessment of the performance of riparian buffers or vegetative filter strips on pollutant removal because it influences the degradation rate of organic compounds such as pesticides (Konopka & Turco, 1991;Munoz-Carpena, Fox, Ritter, Perez-Ovilla, & Rodea-Palomares, 2018). In the SR model, q b was calculated as the downward flux density for the entire soil profile (Equation 6). In the DP model, q b was calculated as the downward flux across the bottom boundary.
The SR model parameters (ω, M, i 0 , C l , and θ e ) and DP model parameters (β, a, b, κ, and K sat ) were inversely calibrated to minimize the weighted average of normalized objective functions (NOF; e.g., Fox, Sabbagh, Chen, & Russell, 2006;Zhou, Wilson, Fox, Rigby, & Dabney, 2016) for q b and θ: where N is the total number of observations of the bottom drainage, q b,obs,i is the ith observation, q b,mod,i is the model estimate for the ith observation and b,obs is the average observation, M is the total number of observations of water content at the measurement depth k, θ k,obs,i is the ith observation at depth k, θ k,mod,i is the model estimate for the ith observation, and θ ,obs, is the average observation. Weights for bottom drainage and water contents b and w θ were set to 0.5 and 0.125, respectively. For the bottom drainage, model performance was also assessed using the Nash-Sutcliffe Efficiency (NSE; Nash & Sutcliffe, 1970): The NSE ranges from −∞ to 1, where 1 indicates a perfect fit. The model performance was considered acceptable if the NSE value was above 0.50, good if the NSE was above 0.65, and very good if the NSE was above 0.75 (Moriasi et al., 2007). For θ, the RMSE of all values was used. To validate the models, the set of parameters that provided the best simulation for the calibration scenario was used to simulate other tests with different rainfall intensities.

Evaluation of model performance
For the calibration of the SR model, NSE was .90 for both configurations, and RMSE was 0.052 and 0.064 for Configurations 1 and 2. The model accurately predicted the total depth of drainage (Table 2). For the DP model, NSE was .78 and .82 and RMSE was 0.070 and 0.079 for Configurations 1 and 2, respectively. Errors on the total depth of drainage were ∼10% (Table 3). For validation scenarios, the calibrated DP and SR models were able to provide relatively accurate predictions of the bottom drainage for different rainfall conditions. The NSE for the SR model ranged from .67 to .79, and errors on the total drainage depth ranged from −12.6 to 9.2% (Table 4). The NSE for the DP model ranged from .75 to .85, and errors on the cumulative drainage depth ranged from 0 to 16.4%. Under a steady rate of infiltration in the macropore domain i s , the SR model predicted that the bottom drainage rate reached a steady state when θ reached θ e (Figure 2). For all the calibration and validation scenar-ios, the bottom drainage rates predicted with the DP model were overestimated during the second half of the test. This was attributed to the decrease of abstractions from the macropore domain to the matrix as the matrix wetted, resulting in a gradual increase in the macropore flux at the bottom. Based on Equation 11, this exchange rate depended on the saturation of both the macropore and matrix domains. However, neglecting the calibration of soil F I G U R E 2 Observed and simulated rates of bottom drainage for validation Scenario 1 (macropore Configuration 2, rainfall intensity = 2.5 cm h −1 ). DP, dual permeability; SR, source responsive.

Vadose Zone Journal
TA B L E 4 Validation performance for the dual-permeability (DP) and source-responsive (SR) models matrix parameters could have led to incorrect estimates of the matrix saturation, and therefore the exchange rate. After the rainfall ended, experimental data showed that bottom drainage continued at a decreasing rate for some time before ending completely. The DP model simulated the flow recession quite accurately (Figure 2). However, the SR model represented by Equation 4 is a simplification for steady-state conditions and therefore cannot simulate water flow recession after the rainfall ended. The original SR model represented by Equation 1 can simulate unsteady conditions, but the step-wise parametrization proposed for the degree of activation of source responsive transport, f, (Equation 4) was not able to simulate unsteady source-responsive fluxes resulting from intermittent rainfall, dynamic infiltration, and exchanges between domains. An improved formulation for f could be obtained from a more mechanistic description of the source responsive flow as film flow along the macropore walls, such as presented by Germann (2020).
In terms of θ, the SR model predicted that θ remained constant until the activation lag time and then reached the equilibrium water content θ e , which was not consistent with observations of final θ varying with depth. Contrary to what was hypothesized, θ e did not seem to be an intrinsic property of the soil. Instead, this parameter was positively correlated to the initial water content θ 0 (R = .94). No significant statistical relationship with the rainfall rate was found. Hincapié and Germann (2009) conducted infiltration experiments in an undisturbed soil column under steady-state surface water input and observed similar water content dynamics. In their experiments, the increase of water content to an approximately constant value called θ max (similar to θ e ) was interpreted as the downward propagation of the wave resulting from the steady-state input. They also found that θ max depended on the initial water content and the amplitude of the water content wave, also dependent on the initial water content. Therefore, θ e represents a local equilibrium in response to the steady-state input and occurs when hydraulic gradients within the area of influence of a macropore became too small to force exchanges between the macropore and the matrix.

F I G U R E 3
Observed and simulated water contents with the dual-permeability (DP) and source-responsive (SR) models for the validation scenario for macropore Configuration 1 (rainfall intensity = 1.5 cm h −1 ): (a) 4-cm depth, (b) 11-cm depth, and (c) 18-cm depth At the 4-cm depth, the DP model simulated the observed increase in θ corresponding to preferential flow and the wetting of the matrix through exchanges with the macropores (Figure 3). The DP model also simulated a second TA B L E 5 Model calibration and validation for the updated source-responsive model (SR2)

Macropore configuration
Rainfall Note. ω, model parameter; M, macropore internal facial area capable of supporting film flow per unit volume of the bulk medium; i 0 , the maximum rate of SR infiltration; C l , model parameter; a t and b t , model parameters; NSE, Nash-Sutcliffe efficiency; q b , bottom drainage discharge; θ, water content.
increase, which was attributed to matrix flow from the surface, but this was not observed by the sensors. This was believed to be caused by a poor parameterization of soil matrix parameters, which were not calibrated but determined from the literature. At 11-and 18-cm depths, the DP model underestimated the increase in θ, which was attributed to the low rate of exchange from the macropores to the matrix simulated by the model (Figure 3). As stated before, this can be caused by a poor parametrization of the soil matrix parameters leading to incorrect exchange rates. In addition, first-order exchange rates are notorious for underestimating highly transient transfer flux at early times (Köhne, Mohanty, Šimůnek, & Gerke, 2004), and the utilization of a constant first-order coefficient κ independent of conditions at the macropore-matrix interface can be questioned.

Refined source-responsive model
The mechanistic prediction of the local equilibrium water content θ e would require matrix flow processes in the vicinity of macropores to be physically described through a two-dimensional radial approach similar to Köhne and Mohanty (2005), which is out of the scope of the SR model. Alternatively, an empirical approach is proposed to estimate θ e in practice. A new formulation for the SR model (SR2) was developed, in which θ e was computed as a linear function of θ 0 : The SR2 model was recalibrated, adjusting parameters ω, M, i 0 , C l , a t , and b t ( Table 5). The SR2 model formulation allowed to accurately simulated the depth of bottom drainage, while also having improved predictions of θ ( Figure 4).

F I G U R E 4
Observed and simulated water contents with the source-responsive (SR) and updated source-responsive (SR2) models for the validation scenario for macropore Configuration 1 (rainfall intensity =1.5 cm h −1 ): (a) 4-cm depth, (b) 11-cm depth, and (c) 18-cm depth Vadose Zone Journal

Interpretation of model parameters
Calibrated parameters for Configuration 1 (fewer macropores) were then compared with those obtained for Configuration 2 (more macropores). For the DP model, the value of the saturated water content in the macropore domain, β, was similar for both configurations (Table 1). It can be noted that these values were one order of magnitude greater than the ratio between the volume of the artificial macropores and the volume of soil, suggesting the existence of other preferential flow pathways in the soil column. The maximum flux in the macropore domain, a, was also similar for both configurations, but the exponent, b, was lower for Configuration 2, suggesting that preferential flow was faster. The first-order coefficient for the exchange rate between the macropores and matrix, κ, was lower for macropore Configuration 1, which was consistent with the lower number of macropores. For the SR2 model, the travel time coefficient C l from the SR2 model was similar between both configurations ( Table 5). The SR fluxes were controlled by the term (M 2 D/Gi 0 )ωI (Equations 1 and 5). The term (M 2 /i 0 )ω could be viewed as a bulk calibration parameter, and it was therefore difficult to interpret specific values obtained for each component i 0 , M, and ω. The term (M 2 /i 0 )ω was greater for Configuration 2 than for Configuration 1. This meant that for the same rainfall intensity, greater SR fluxes would be computed for Configuration 2. Also, values for a t and b t were relatively similar between configurations.

CONCLUSIONS
The performance of the SR model and DP model calibrated for the macropore domain only was evaluated based on experimental data of water flow through a soil with vertical, artificial macropores, for which preferential flow was the dominant process. The DP model provided relatively accurate predictions of bottom drainage rates. However, an inaccurate parametrization of the soil matrix led to poor predictions of the downward matrix flow and of the exchanges between the macropore and matrix domains. The reduced calibration approach appeared valid to simulate subsurface fluxes when matrix flow was negligible; for example, deep fluxes in highly macroporous soils during short storm events. A refined SR model incorporating a mass balance component and an empirical parametrization for the equilibrium water content, θ e , provided accurate predictions of bottom drainage rates and water contents under steadystate rainfall but was not capable of simulating the flow recession once the rainfall ended. The current model for-mulation for the source-responsive flow is derived from steady flow assumptions and cannot simulate intermittent rain and dynamic infiltration scenarios, which limits the model applicability to steady-state scenarios with negligible matrix flow. A mechanistic description of film flow along macropore walls is required in order to improve the SR model applicability to practical scenarios.
In both DP and SR models, the variation in parameter values between macropore configurations was consistent with the variation in macropore properties, but the limited number of configurations tested made it difficult to interpret the specific effect of each parameter on the model response. Additional experiments representing different configurations and contributions of preferential flow would be beneficial to better infer the role of model parameters and draw relationships between parameter values and macropore properties.

A C K N O W L E D G M E N T S
This work was funded by USDA National Institute of Food and Agriculture (NIFA) Project no. 2016-67019-26855. We want to thank Dr. Rich McLaughlin from the Department of Crop and Soil Sciences at North Carolina State University for lending us the rainfall simulator used in this study, and the Research Shop of the Department of Biological and Agricultural Engineering at North Carolina State University for building the column.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.