|
Assessment of Macropore Component of RZWQM2 in Simulating Hourly Subsurface Drainage and Peaks
Ziwei Li1, Changchi Xian1, Zhiming Qi1,*, Liwang Ma2, Matthew W. Sima3, Matthew Helmers4, Tiequan Zhang5, Rob Malone6, Quanxiao Fang7
Published in Journal of the ASABE 66(5): 1303-1315 (doi: 10.13031/ja.15530). 2023 American Society of Agricultural and Biological Engineers.
1Department of Bioresource Engineering, McGill University, Sanite Anne de Bellevue, Quebec, Canada.
2Rangeland Resources and Systems Research Unit, USDA ARS, Fort Collins, Colorado, USA.
3Department of Civil and Environmental Engineering, Princeton University, Princeton, New Jersey, USA.
4Agricultural and Biosystems Engineering, Iowa State University, Ames, Iowa, USA.
5Harrow Research and Development Center, Agriculture and Agri Food Canada, Harrow, Ontario, Canada.
6National Laboratory for Agriculture and the Environment, USDA ARS, Ames, Iowa, USA.
7Qingdao Agricultural University, Qingdao, Shandong, China.
*Correspondence: zhiming.qi@mcgill.ca
Submitted for review on 11 January 2023 as manuscript number NRES 15530; approved for publication as a Research Article and as part of the “Advances in Drainage: Selected Works from the 11th International Drainage Symposium” Collection by Community Editor Dr. Mohamed Youssef of the Natural Resources & Environmental Systems Community of ASABE on 25 July 2023.
Mention of company or trade names is for description only and does not imply endorsement by the USDA. The USDA is an equal opportunity provider and employer.
Highlights
Abstract. Understanding preferential flow through soil macropores is critical to effectively managing subsurface drainage water quantity and quality. This study aims to assess the macropore component of the Root Zone Water Quality Model (RZWQM2) in simulating subsurface tile flow with a high time resolution. Observed hourly tile flow rates from two experimental sites in Ontario, Canada (2008-2011) and Iowa, USA (2007-2008) were used to evaluate the importance of including a macropore flow component in subsurface drainage simulation. Activating the macropore component in the model improved the simulation of hourly drainage peaks, especially peak amplitude. Still, it did not improve the simulation of the total drainage amount for each rainfall event. Simulation of the drainage peak recession varied from peak to peak, suggesting that further studies are warranted for drainage flow in the model. In general, the macropore component in the RZWQM2 model improved subsurface peak subsurface simulation at the hourly resolution. However, further investigation and model modifications are needed to improve the drainage simulation’s timing and quality for RZWQM2’s hydrologic simulation of macropore flow and subsurface drainage.
Preferential flow is a common phenomenon in structured soils through macropores or due to flow instability (e.g., finger flow) (Weiler, 2017) and is one of the key factors dominating nutrient translocation in soils (Julich et al., 2022). In agricultural soils, macropores are formed under various circumstances, including burrowing animals, plant roots, soil drying cracks, etc. (Beven and Germann, 2013). In the presence of soil macropores, the capillary tube analogy in soil hydraulics and the steady equilibrium assumptions of Darcy and Richards’s equations may not be applicable (Beven and Germann, 2013; Weiler, 2017). However, the Richards equation is still the most commonly used algorithm in system models due to its basic idea of balancing all potentials in an elementary volume. Weiler (2017) called for a new theory on macropore flow and new views to see preferential flow as the modeling norm, not the exception.
The dual-permeability or dual-porosity model is generally used and adopted in the literature (Orozco-López et al., 2018; Weiler, 2017; Jarvis et al., 2016). Each domain (mobile or immobile pore) is simulated separately in most models, but with interactions among the pores during infiltration (Malone et al., 2001; Orozco-López et al., 2018). Practically, simulation of macropore flow is difficult because: (1) macropores are dynamic and not easy to quantify (Reck et al., 2018; Liu et al., 2014); (2) connectivity and continuity of macropores are in general unknown (Jarvis et al., 2016); (3) effectiveness of macropores in transporting water and chemicals changes with antecedent soil water content, rainfall intensity, and management (Hardie et al., 2011; Malone et al., 2001, 2004a,b; Nimmo, 2016); (4) initialization of macropore flow is less quantified (Nimmo, 2016; Klaus et al., 2013); and (5) interactions between macropores and soil matrix are not fully understood (Malone et al., 2004a,b; Nimmo, 2016; Arora et al., 2012).
Macropore flow directly contributes to subsurface drainage and chemical losses in tile flow (Malone et al., 2001; Klaus et al., 2013; Nazari et al., 2020) and is responsible for high peak concentration and the early arrival of chemicals in tile drainage (Malone et al., 2004a,b; Tiktak et al., 2012). Although Klaus et al. (2013) found that up to 20% of tile drainage was from macropore flow while the rest was from soil matrix flow, chemical and bacterial losses in tile drainage were mainly from macropore flow (Fox et al., 2004, 2007, 2012). A lysimeter experiment showed that the solute exchange between mobile and immobile domains was impeded under shallow water table conditions, thus increasing the risk of groundwater pollution (Mencaroni et al., 2021). Some macropores may connect to drainage tiles directly, creating an express path for water and chemicals (Malone et al., 2001; Fox et al., 2004). In general, accurate simulations of chemical transport in tile flow are dependent on accurate subsurface drainage (Malone et al., 2004a,b), and the latter is affected by rainfall intensity, soil matrix properties, and macroporosity (Fox et al., 2007; Nimmo, 2016; Que et al., 2018; Gao et al., 2018).
Macropore flow velocity may be simulated by Poiseulle’s law (Ahuja et al., 2000), kinematic wave (Que et al., 2018; Orozco-López et al., 2018), Stokes’s law (Orozco-López et al., 2018), and power law (Jarvis et al., 1994). However, Gao et al. (2018) found that Poiseulle’s and the kinematic wave (Manning equation) overestimated macropore flow velocity. Furthermore, models vary considerably in macropore dynamics, macropore flow mechanism, macropore flow initialization, macropore-soil matrix interaction, and macropore connectivity. Macropore flow is considered in many hydrological models and applied to tile drainage simulation, such as APEX (Ford et al., 2017), HYDRUS (Šimunek et al., 2012), MACRO (Jarvis and Larsbo, 2012), MIKE SHE (Zhou et al., 2013), PEARL (Tiktak et al., 2012), and RZWQM (Malone et al., 2001; Fox et al., 2004).
Among the above-mentioned process-based models, RZWQM2 is a comprehensive one-dimensional agricultural systems model that covers a detailed description of a wide array of processes from the upper boundary (atmosphere) to the bottom layer (soil). RZWQM2 could also be used to study the interaction of physical, chemical, and biological processes within the soil profile under various management practices (Ahuja et al., 2000), including controlled drainage, different types of irrigation, inorganic and manure fertilization, tillage, pesticide use, cover crop, crop rotation, and residual management (Ma et al., 2012). Accurate simulation of subsurface drainage with RZWQM2 could facilitate further evaluation of the effect of various agricultural best management practices on agricultural subsurface drainage (Heilman et al., 2006; Jeong and Bahttarai, 2018; Jiang et al., 2018; Sadhukhan et al., 2019). RZWQM2 uses a simple algorithm for macropore flow simulation to initialize macropore flow by assuming rainfall intensity above saturated soil hydraulic conductivity to be macropore flow. In addition, users can specify dead macropores for each soil layer to take care of connectivity issues and express fractions connecting directly to drainage tiles for fast chemical transport. The Radial Green-Ampt equation is used for macropore-soil matrix water exchange (Ahuja et al., 2000; Fox et al., 2007).
The macropore component in RZWQM2 has been proven to be satisfactory in simulating tile drainage and associated pesticide and bacteria losses (Malone et al., 2004a,b; Fox et al., 2004, 2007, 2012), but it also suggested that hourly rainfall should be used to improve macropore flow simulation, especially for the initiation of macropore flow (Malone et al., 2004b; Fox et al., 2007). However, few studies have compared simulated subsurface drainage under hourly rainfall with and without macropores at multiple locations. Therefore, the overall objective of this study was to test the macropore component of RZWQM2 at an hourly time step, focusing on drainage peak periods where the macropore component would be activated, using observed data from Iowa, USA, and Ontario, Canada, and to understand the role of macropore flow in hourly tile drainage flow.
Materials and Methods
Ontario and Iowa Drainage Sites
Two sets of data were used to evaluate the macropore component of the RZWQM2 model. The first experiment was conducted by Tan et al. (2009) at an experiment field located on the north shore of Lake Erie near Harrow, Ontario. The second field study was conducted at the Agricultural Drainage Water Quality – Research and Demonstration Site (ADWQ-RDS) near Gilmore City in Pocahontas County, north-central Iowa. To assess the effects of considering macropores on the accuracy of high-resolution time scale drainage peak simulation, we selected three periods (two from Ontario and one from Iowa) during which macropore flow occurred. Since high infiltration and preferential flows were not frequent in the scenarios, to better assess and demonstrate the macropore component, only the periods having drain flow peaks equal to or exceeding 4.0 mm hr-1 were selected and evaluated (table 1).
Initiated in the spring of 2008, the experimental site in Ontario consisted of 16 drained plots (15.2 m × 67.1 m each, Tan et al., 2009). Buffer zones and impermeable barriers were set to prevent interaction among plots. In this study, two plots with regular drainage were used to test the macropore component of RZWQM2. The predominant soil in these fields was Brookston clay loam (fine-loamy, mixed, superactive, mesic Typic Argiaquolls), with 28% sand, 37% silt, and 35% clay and a mean bulk density of 1.34 Mg m-3. The measured soil porosity was 52.4% and 
ranged from 17 to 119 mm d-1, averaging 50 cm d-1 (Liu et al., 2011). Hourly observed drainage outflow from the two experimental plots recorded during 2008 and 2012 was averaged across the plots. Hourly weather data, including precipitation, air temperature, solar radiation, wind speed, and relative humidity, were recorded at the Whelan weather station, located less than 0.5 km from the experiment site. However, precipitation measurements at the Whelan station were unavailable from 1 October 2008 to 30 April 2009, and therefore precipitation measured at the Harrow weather station located 16.6 km from the experiment site was used.
Initiated in the fall of 2004 and continuing for five years, the Iowa field experiment was arranged in a completely randomized block design with 78 individually drained plots (15.2 m × 38.0 m each, Lawlor et al., 2008) replicated in each of the four blocks selected for their contrasting long-term drainage performance (e.g., high, medium-high, medium-low, and low drainage). In each block, plots were randomly assigned a specific land use/land cover treatment, resulting in each treatment, in this case, corn (Zea mays L., odd years)-soybean (Glycine max L. Merrill, even years) rotation with winter rye (Secale cereale L.) cover crop growing between the harvest and planting of summer crops with four replicates. Hourly drainage from each block was monitored, and the mean values across the four blocks served as the observed hourly drainage data used in this study.
| Table 1. Peak drainage flow periods and individual events employed in the hourly analysis of RZWQM2 model accuracy (with and without the macropore component) for peak flow scenarios at Iowa, USA, and Ontario, Canada sites. | |||||||
| Site | Peak Designation | Period | Peak | ||||
| No. | Date of Period | No. | Time | Observed Peak Drain Flow Rate (mm hr-1) | |||
| Iowa |  | 1 | 19 – 23 August 2007 | 1 | 21 August 2007 20:00 | 6.0 | |
| Ontario |  | 1 | 27 – 29 June 2008 | 1 | 28 June 2008 1:00 | 4.6 | |
 | 1 | 27 – 29 June 2008 | 2 | 28 June 2008 15:00 | 5.2 | ||
 | 2 | 25 – 27 May 2011 | 1 | 26 May 2011 10:00 | 4.0 | ||
In the Iowa scenario, the hydraulic properties, soil bulk density (?), particle size distribution, saturated hydraulic conductivity (ksat), and soil water characteristic curves of the predominantly fine loam soils were determined from undisturbed soil cores. Residual water content and soil water content (?) at different matric potentials (10, 33, and 1500 kPa) were interpolated or extrapolated from the soil water characteristic curve. Bubbling pressure and pore size distribution were fitted using the Brooks-Corey equation (Brooks and Corey, 1964). An onsite weather station monitored hourly meteorological data during the five-year experiment period, including precipitation, air temperature, solar radiation, relative humidity, and wind speed.
Macropore Component in RZWQM2
Developed in 1993, the macropore component in RZWQM2 was included to allow the simulation of preferential flow in agricultural lands (Ahuja et al., 2000). The macropore soil parameters include sorptivity factor, macropore to tile drain expression fraction, and effective lateral infiltration wetting thickness (Malone et al., 2004a,b). In addition, for each soil horizon, macropore-related parameters include macroporosity, pore radius, width and length of the crack, length of aggregate, and the fraction of dead-end pores. Macropore flow that was simulated by the RZWQM2 was assumed to be rapid and unaffected by the pore tortuosity. The topsoil horizons are assumed to have cylindrical macropore channels. The bottom soil horizons are assumed to have planar cracks (Ahuja et al., 2000). The macroporosity and pore radius in each layer describe the macropores quantitatively and determine the total number of pores in the field. Depending on these two parameters, RZWQM2 computes the density of macropores per unit area as (Ahuja et al., 2000):
(1)
where
nmacro = number of macropores per unit area, which are assumed to be evenly distributed in the soil layers
r = radius (cm) in each soil horizon within which macropores are found
Pmac = macroporosity in each soil horizon (cm3 cm-3).
According to the model, rainfall or irrigation first infiltrates the soil matrix, according to the Green-Ampt equation. The excess water is then allowed to enter macropores, and the reminder goes off as surface runoff. Based on Poiseulle’s law, the maximum capacity of macropore flow (cm h-1) for the cylindrical holes 
and planar cracks 
were computed assuming gravity flow as (Ahuja et al., 2000):
(2)
(3)
where
r = radius of cylindrical pores (cm)
d = width of planar cracks (cm)
g = gravitational constant (cm h-2)
Pmac = continuous macroporosity (cm3 cm-3) as a fraction of the soil volume
= dynamic viscosity of water (g cm-1 h-1)
? = density of water (g cm-3).
Continuous macropores are assumed to extend vertically in soil, and a portion of them may be blocked at each soil horizon. The flow is sequentially routed downward at each timestep through the continuous macropores (1 cm increments). Macropore flow is allowed to flow into the dead-end macropores. The soil matrix could absorb macropore flow that streams into both the dead-end and continuous macropores by radial or lateral infiltration. Lateral flow by the soil matrix adsorption is time-dependent and only occurs in the drier soil matrix below the transient wetting front. The infiltration rates (cm h-1) for the cylindrical macropores (Vc) at the first-time step and the rest of the timesteps, and planar cracks (Vp) in the soil layers were computed as (Ahuja et al., 2000):
(4)
(5)
(6)
where
ksat = saturated hydraulic conductivity (cm h-1)
?t1 = first-time step in the model calculation (h)
t = cumulative time for lateral flow (h)
Hc = capillary drive term for the soil matrix (cm)
rwf = wetted radius at any given time (cm)
rp = macropore radius
(?sat – ?i) = soil moisture deficit in the particular soil depth range, namely the soil moisture content at saturation (?sat) minus the soil moisture under present conditions (?i) (cm3 cm-3).
The water reaching a dead end can be absorbed or stored, while the remaining water will be routed to the water table. When the macropore flow first reaches the top of the water table, a small fraction of water, containing the initial solute concentration from the soil surface, is assumed to go directly to the tile drains (Fox et al., 2004; Malone et al., 2001).
Multiple field studies demonstrated that some macropores could have direct hydraulic connections towards the subsurface drains (Shipitalo et al., 2004; Akay and Fox, 2007; Fox et al., 2012). Express fraction (EF) is a parameter designed to represent the percentage of macropores that are directly connected to the subsurface drainage and was included in the RZWQM2 by Fox et al. (2004). The parameter has been demonstrated to be a crucial parameter in subsurface drainage in RZWQM2 simulation, especially for simulating the chemical transports in the subsurface drainage (Kumar et al., 1998; Fox et al., 2004; Fox et al., 2007).
Parameterization for the Ontario Scenario
Lu (2015) calibrated RZWQM2 based on measured subsurface drainage and surface runoff data at the Ontario site, in which the experimental period (from 1 June to 22 December 2011) was divided into 17 time periods. Soil parameters were calibrated to match the accumulative observations in those 17 time periods, which might be inadequate to simulate hourly observation. Therefore, we recalibrated the hydraulic parameters to match hourly observations (table 2).
| Table 2. Measured and calibrated soil hydraulic properties at the Ontario site.[a] | |||||||||||
| Measured | Calibrated | ||||||||||
Soil Moisture Content, ? (mm3 mm-3)Soil Matric Potential,  (kPa) | |||||||||||
| Depth (m) | ? (Mg m-3) | Sand (g g-1) | Silt (g g-1) |  (mm h-1) |  (mm h-1) | Tsat 0 kPa | T10 -10kPa | Tfc -33 kPa | Tp -1500kPa | Tr -8 kPa | |
| 0-0.25 | 1.326 | 0.299 | 0.363 | 9.2 | 17 | 0.5 | 0.383 | 0.325 | 0.198 | 0.040 | |
| 0.24-0.45 | 1.391 | 0.238 | 0.349 | 38 | 68 | 0.475 | 0.378 | 0.3363 | 0.240 | 0.090 | |
| 0.45-0.80 | 1.391 | 0.257 | 0.33 | 30 | 60 | 0.475 | 0.371 | 0.3299 | 0.236 | 0.090 | |
| 0.80-1.20 | 1.391 | 0.243 | 0.359 | 20 | 40 | 0.475 | 0.390 | 0.347 | 0.246 | 0.090 | |
| 1.20-3.00 | 1.391 | 0.243 | 0.359 | 5 | 20 | 0.475 | 0.390 | 0.347 | 0.246 | 0.090 | |
| 3.00-3.09 | 1.391 | 0.243 | 0.359 | 0.1 | 20 | 0.475 | 0.390 | 0.347 | 0.246 | 0.090 | |
[a] ? = soil bulk density; sand (50–2000 µm) and silt (2–50 µm) fractions; | |||||||||||
and 
, vertical and lateral saturated hydraulic conductivity, respectively. ?sat, ?fc, and ?pwp are soil moisture at saturation, field capacity, and permanent wilting point, respectively. ?r is the residual soil water content.The average lengths of macropores were assumed to decrease, and the dead-end pore fraction was assumed to increase with increasing soil depth. The sorptivity factor parameter was set to 0.5, based on the results of a previous study (Kumar et al., 1998). EF was calibrated against seven values: 0 (no macropore had direct connection towards drainage), 0.01, 0.02 (previously calibrated value by Fox et al. [2004] and Kumar et al. [1998]), 0.03, 0.04, 0.05, and 0.1. Hourly drainage simulation results under different EF values were similar. The change in the mean drainage was ± 0.01 mm for periods 1 and 2. NSE performances for hourly drainage simulation were also similar. The difference could only be observed from the third digit, with EF of 0.1 having the worst NSE (period 1: 0.456, period 2: 0.248) and EF of 0.02 giving the best NSE (period 1: 0.464, period 2: 0.252). For the long-term simulation (compared with 17 observation periods), the NSE performance remained the same at 0.46. EF of 0.02 gave a slightly better result during the hourly drainage simulation. As a result, similar to Kumar et al. (1998) and Fox et al. (2004), an EF value of 0.02 was used for the Ontario site. All other parameters (macroporosity, pore radius, crack lengths, crack widths, and proportions of dead–end pores) were calibrated to match hourly observed drainage flow peaks over the full experiment period (table 3) based on the PBIAS, NSE, and IoA statistics.
Parameterization for the Iowa Scenario
Soil hydraulic parameters for RZWQM2 were measured or calibrated by Qi et al. (2011) on a daily scale without consideration of macropore flow (table 4). Soil bulk density, particle size distribution and ksat (i.e., vertical hydraulic conductivity 
) were derived from field measurements and lateral saturated hydraulic conductivities 
in each layer were adjusted to 
, in order to match the peaks of daily drain flow. Soil water retention curve parameters were measured in the laboratory using soil columns but were subsequently calibrated to match the measured soil water content better. Parameter recalibration was not needed for the Iowa site since the current calibrated parameters were satisfactory for hourly simulation as determined by PBIAS, NSE, and IoA.
| Table 3. Calibrated parameters for the RZWQM2 macropore component at the Ontario site. | |||||||||||||||||
| Soil Depth (m) | Macroporosity (mm3 mm-3) | Radius of Cylindrical Pores (mm) | Width of Cracks (mm) | Length of Cracks (mm) | Average Length of Aggregate (mm) | Fraction of Dead-End Pores | |||||||||||
| 0-0.25 | 0.0003 | 1 | 0 | 0 | 100 | 0.01 | |||||||||||
| 0.25-0.45 | 0.0003 | 0 | 1 | 100 | 100 | 0.3 | |||||||||||
| 0.45-0.80 | 0.0003 | 0 | 1 | 50 | 50 | 0.5 | |||||||||||
| 0.80-1.20 | 0.0003 | 0 | 1 | 50 | 50 | 0.8 | |||||||||||
Table 4. Measured and calibrated soil hydraulic properties at the Iowa site.[a] | |||||||||||||||||
| Measured | Calibrated | ||||||||||||||||
Soil moisture content, ? (mm3 mm-3)Soil matric potential,  (kPa) | |||||||||||||||||
| Depth(m) | ? (Mg m-3) | Sand (g g-1) | Silt (g g-1) |  (mm h-1) |  (mm h-1) | Tsat 0 kPa | T10 -10kPa | Tfc -33 kPa | Tpwp -1500kPa | Tr -8 kPa | |||||||
| 0-0.10 | 1.37 | 0.32 | 0.36 | 48 | 97 | 0.482 | 0.383 | 0.376 | 0.189 | 0.071 | |||||||
| 0.10-0.20 | 1.38 | 0.32 | 0.36 | 33 | 66 | 0.476 | 0.384 | 0.376 | 0.230 | 0.072 | |||||||
| 0.20-0.30 | 1.39 | 0.33 | 0.53 | 51 | 101 | 0.473 | 0.384 | 0.376 | 0.200 | 0.079 | |||||||
| 0.30-0.40 | 1.39 | 0.40 | 0.30 | 40.8 | 82 | 0.474 | 0.384 | 0.399 | 0.212 | 0.072 | |||||||
| 0.40-0.60 | 1.39 | 0.46 | 0.30 | 40.8 | 82 | 0.474 | 0.408 | 0.368 | 0.218 | 0.065 | |||||||
| 0.60-0.90 | 1.45 | 0.44 | 0.34 | 26.4 | 53 | 0.450 | 0.380 | 0.368 | 0.204 | 0.034 | |||||||
| 0.90-1.20 | 1.46 | 0.44 | 0.34 | 26.4 | 53 | 0.450 | 0.312 | 0.299 | 0.184 | 0.033 | |||||||
| 1.20-2.00 | 1.46 | 0.44 | 0.34 | 26.4 | 53 | 0.450 | 0.310 | 0.299 | 0.168 | 0.033 | |||||||
| 2.00-3.00 | 1.50 | 0.44 | 0.34 | 26.4 | 53 | 0.450 | 0.310 | 0.299 | 0.168 | 0.033 | |||||||
| 3.00-3.90 | 1.50 | 0.44 | 0.34 | 0.001 | 50 | 0.450 | 0.310 | 0.299 | 0.168 | 0.033 | |||||||
[a] ? = soil bulk density; sand (50–2000 µm) and silt (2–50 µm) fractions; | |||||||||||||||||
and 
, vertical and lateral saturated hydraulic conductivity, respectively. ?sat, ?fc, ?pwp, are soil moisture at saturation, field capacity, and permanent wilting point, respectively. ?r is the residual soil water content.Similar to the Ontario site, the sorptivity factor was set to 0.5 based on the calibrated values from Kumar et al. (1998) and Fox et al. (2004). Seven EF values were tested, the same as the Ontario site’s EF calibration. Similar to the Ontario site, changes in subsurface drainage for the Iowa site under different EF values were minimal. For the hourly subsurface peak drainage simulation from 19 August to 23 August, 2007, the change in the mean and sum of the simulated hourly subsurface drainage during the period was from -0.014 mm and -1.9 mm (EF of 0) to 0.012 mm and 1.7 mm (EF of 0.1) compared to using an EF of 0.02. With the increase of EF, the drainage also increased but was only noticeable (increase greater than 0.1 mm) where drainage peak occurred (Aug/20th and Aug/21st). The NSE performance for the cumulative drainage for the period was similar for EF of 0, EF of 0.01, and EF of 0.02 (0.37, 0.36, 0.36), which were better than the larger EF values applied (worst NSE value by EF of 0.1 of 0.34). At daily resolution, NSE remained the same (0.73) for long-term subsurface drainage across different EF settings. Considering that EF of 0 provided a slightly better NSE result, it was considered that EF was not needed for the Iowa site. A similar situation was reported by Fox et al. (2007), where no EF was needed for the Baton Rouge site for subsurface drainage simulation using RZWQM2. The length of cracks and aggregates was assumed to decrease while the dead-end pore fractions increased. Other parameters were calibrated better to characterize the hourly drainage peak (table 5).
Statistical Methods for Model Accuracy Evaluation
To evaluate the performance of the model, PBIAS, NSE, and Index of Agreement IoA were calculated:
(6)
(7)
(8)
where
n = number of measured (or simulated) values
Oi = ith observed value
= mean observed value
Si = ith simulated value.
In this study, PBIAS within ±15% (eq. 6), NSE > 0.70 (eq. 7), and IoA close to 1 (eq. 8) indicate that the model performance is good for a daily time step in RZWQM2 (Moriasiet al., 2015).
Results and Discussion
Comparison Between Simulations with and without Macropores in Ontario
Simulated drain flows during the 17 observation periods were not significantly different with the macropore component (MP) and without (NoMP) activated (fig. 1). Compared with observed data, total drainage volume was underestimated with NoMP and MP. Activating MP, which provided channels for more water to penetrate, only slightly improved predicted total drain flow, with a 6% increase in total predicted drainage. The NSE and IoA statistics were slightly decreased after activating the macropore component, with NSE dropping from 0.48 to 0.46 and IoA from 0.74 to 0.72 (table 6). However, the changes in these two statistics were insignificant (within 5% of the original values) because only a few events over the entire simulation had macropore flow.
 |
MP significantly improved peak drainage simulations when using hourly drainage data to evaluate model performance during drainage peak periods (fig. 2). For the two peaks 
in Period 1 (27-29 June 2008; table 1) at the Ontario site, the saturated surface soil layer and limited infiltration rate led to a severe underestimation of the two peaks for the NoMP option (PBIAS = -39.97% and -51.59%, respectively). With MP, a portion of excess water can bypass the soil profile and enter the drainage tile directly. As a result, the MP simulation showed BIAS values of -11.28% and -20.85% for the two peaks (table 7). Similarly, the PBIAS of the single predicted peak 
in Period 2 (25-27 May 2011) decreased from -50.69% to -31.35% under MP (fig. 3, table 7).
| Table 6. Model accuracy statistics for no macropore (NoMP) and with macropore (MP) simulation of long-term drainage with RZWQM2, Ontario, and Iowa sites. | |||||
| Model Accuracy Statistics[a] | Ontario | Iowa | |||
| NoMP | MP | NoMP | MP | ||
| PBIAS | -18.20% | -13.19% | 6.47% | 10.43% | |
| NSE | 0.48 | 0.46 | 0.71 | 0.73 | |
| IoA | 0.74 | 0.72 | 0.76 | 0.76 | |
[a] PBIAS, percent bias; NSE, Nash-Sutcliffe Efficiency; IoA, Index of Agreement | |||||
However, the overall quantity and timing of drain flow simulated with the MP option were not as good as those simulated with the NoMP option. The PBIAS increased from 56.52% to 89.49% in Period 1 (peaks 
), and from -0.76% to 5.46% in Period 2 (peak 
) (table 7). There was a larger overestimation in the drainage peak in MP because of the preferential flow in macropores. Both measured and simulated drainage showed a peak followed by a recession (figs. 2 and 3). However, the simulated peak height and recession duration vary from event to event and location to location. For example, the tailing of peaks was overpredicted for Period 1 (fig. 2), whereas it was adequately predicted for Period 2 at the Ontario site (fig. 3).
Comparison Between MP and NoMP in Iowa
First, the effects of macropore flow on drainage simulations were analyzed daily (fig. 4). The predicted drainage values from the NoMP and MP options were almost identical except for the one high drainage day (21 August 2007). Macropore flow was only generated when the model simulated “runoff,” the only case on 21 August 2007. Thus, the simulated 36 mm of macropore flow resulted in 32 mm greater tile drainage. As such, the simulation with the MP component activated in RZWQM2 over predicted cumulative discharge over the two years by 3.7% compared to when MP was not activated (table 6). However, large drainage peaks were better simulated by the MP than the NoMP option, resulting in a slightly better performance of the MP option in predicting drainage peaks, raising the NSE from 0.71 to 0.73 and the IoA from 0.757 to 0.762 (table 6).
Based on the hourly results from the Iowa site (table 7), MP simulated a higher peak value than NoMP, which improved the IoA. Due to the preferential flow, the MP model predicted a sharp peak of about 10 mm on 21 August after an intensive rainfall, while the NoMP simulation underestimated this peak at only 3 mm. During Period 1 in Iowa (fig. 5), the second peak showed the most significant difference between MP and NoMP simulations during a continuous rainfall event. The first peak (early morning on 20 August 2007) was well simulated, and there was no difference between MP and NoMP because there was no runoff. However, when the second intensive rainfall event occurred (about 44 hours later), the excess water beyond infiltration capacity was directed into macropore flow under MP but as surface runoff under NoMP. The former increased drainage peak occurred on 20 August 2007.
 |
Figure 2. Hourly rainfall, observed drainage, and drainage simulated by RZWQM2 without macropore (NoMP) and with macropore (MP) component for the Ontario site, Period 1, peaks  and  . |
| Table 7. Model accuracy statistics for predicted drainage peak, amount, and timing by RZWQM2 with no macropore component (NoMP) and with macropore component (MP) Ontario and Iowa sites. | |||||
| Peak(s) | Model Accuracy Statistic[a] | Simulated | |||
| NoMP | MP | ||||
 | PBIAS | -39.97% | -11.28% | ||
 | PBIAS | -51.59% | -20.85% | ||
 | PBIAS | -50.69% | -31.35% | ||
 | PBIAS | -51.92% | 67.63% | ||
| Cumulative value (mm over Period) | |||||
 | PBIAS | 0.76% | 5.46% | ||
| NSE | 0.48 | 0.25 | |||
| IoA | 0.58 | 0.52 | |||
 [b] | PBIAS | -0.76% | 5.46% | ||
| NSE | 0.47 | 0.46 | |||
| IoA | 0.77 | 0.75 | |||
 | PBIAS | -25.07% | -7.40% | ||
| NSE | 0.54 | 0.36 | |||
| IoA | 0.68 | 0.71 | |||
[a] PBIAS, percent bias; NSE, Nash-Sutcliffe Efficiency; IoA, Index of Agreement. [b] Includes | |||||
and peak at 18:00 on 25 May 2011.Although PBIAS for simulated peak values was not substantially improved by using the MP vs.NoMP model (PBIAS = 67.63% vs. -51.92%, respectively, table 7), the cumulative drainage volume for this period was better predicted with MP than withNoMP (PBIAS = -7.40% vs. -25.07%, IoA = 0.71 vs. 0.68, NSE = 0.36 vs. 0.54, respectively). Because the NSE value is dependent on the hourly comparisons of observed and simulated values, the results indicated that NoMP better simulated the timing of the peaks than MP. The simulated drainage peak with MP (10.5 mm) occurred one hour earlier than the observed peak (when the observed drainage was low, 2.3 mm) and has a much lower NSE. The one-hour advance in the simulated peak obtained with MP could be due to the macropore flow velocity calculated by Poiseulle’s law being faster than it was in the soil due to tortuosity and noncontinuous macropores.
Sensitivity Analysis
Sensitivity analyses were conducted to investigate the effects of macroporosity (macro) and pore radius (r) on total macropore flow. The Ontario and Iowa site simulation results were similar (figs. 6 and 7). The predicted macropore flow was relatively insensitive to the two macropore parameters. Thus, the model’s macropore component was not sensitive enough to quantitatively mimic the actual macropore structure. According to the equations in RZWQM2, increasing macroporosity can provide more open channels for surface water to flow through, thereby increasing macropore flow. Conversely, as shown in equation 1, a larger pore radius can result in fewer macropores, which can accept excess water and lead to lesser macropore flow. At the Ontario site, reducing macroporosity by 10%-30% resulted in a 0.5% increase in total macropore flow. Increasing the macropore radius by 10%-30% would, in theory, reduce macropore flow; instead, the flow increased by 0.5%. Similar fluctuations occurred for the Iowa simulations. With a ±30% change in either parameter, there were only slight changes (<0.5%) to macropore flow.
 |
Figure 3. Hourly rainfall, observed drainage, and drainage simulated by RZWQM without (NoMP) and with (MP) macropore component for the Ontario site, Period 2, peak  . Note that  is the peak occurring on 26/5/2011 at 16:00. |
 |
| Figure 4. Observed daily drainage flow and daily drainage flow simulated by RZWQM2 without (NoMP) and with (MP) macropore component for the Iowa site. |
Possible reasons for this unexpected effect include the following:
- The balance between lateral absorption and vertical macropore flow may be altered. When macroporosity increases or the radius decreases, the total number of macropores increases. Thus, there is more lateral absorption of water by the surrounding soil. Conversely, if the change in absorption is more significant than the change in water flow, the total macropore flow could be reduced.
- Changing the parameters might not only result in a different amount of water in the macropore but might also affect the velocity of the macropore flow. At each time step, water in the macropore channels is routed downward to deeper soil layers, and then the lateral absorption rate will depend on the soil moisture deficit in that specific layer. According to equations 2 and 3, macroporosity and pore radius affect the maximum macropore flow capacity, resulting in different depths at which lateral absorption occurs. Therefore, the total macropore flow would be influenced.
- Since the system components in the model are interactive, changing parameters in the macropore component of the model will also affect soil matrix water flow. As a result, plant root uptake, soil surface evaporation, and plant transpiration could be altered, and these factors may also impact macropore flow.
 |
Figure 5. Hourly rainfall, observed drainage, and drainage simulated by RZWQM2 without (NoMP) and with (MP) macropore component for the Iowa site, Period 1 (19-23 August 2007), peak  . |
Although macropore flow fluctuations were caused by small changes (within 30%) in parameters, the main effects from macroporosity and radius changes can still be observed. Those fluctuations were always within 0.5% when changes in parameters were <30%; however, when the adjustments in parameters increased to 50%, the effects of these two parameters became obvious. For example, the macropore flow changes increased at 50%, ranging from 1% to 4% in Ontario. The total amount of flow tended to increase with higher macroporosity and decrease with a larger pore radius. Nonetheless, total macropore flow was generally not sensitive to macroporosity and pore radius, with a 50% change in these parameters resulting in a < 4% change in total macropore flow.
 |
| Figure 6. Macropore flow sensitivity to changes in RZWQM2 macropore component parameters at the Ontario site. |
 |
| Figure 7. Macropore flow sensitivity to changes in RZWQM2 macropore component parameters at the Iowa site. |
Sensitivity analyses were also conducted for drainage peak duration on an hourly scale. Results suggest that the main effects of input parameters (macroporosity and radius) on predicted drainage are in the rising limb (front) of the drainage flow hydrograph. For the two periods monitored in Ontario and one in Iowa when macropore flow occurred, parameters ranging from 50% to 150% of the original calibrated values were factorially combined as input for the model, and the simulated mean values and the standard errors were shown in figure 8. An error bar was displayed on each hourly simulation (±3× standard errors). In most
 |
| Figure 8. Hourly subsurface drainage sensitivity to changes in RZWQM2 macropore component parameters at the Ontario and Iowa site. (a) for the Ontario site’s drainage peak period 1, (b) for the Ontario site’s drainage peak 2, and (c) for the Iowa site’s drainage peak period. |
simulations, the ratio of 3× standard error to mean (R3s.e./mean) was lower than 1%, which indicated that the variability of simulated drainage was small even when input parameters varied across a relatively larger range (50% to 150%). For Ontario, Period 1 (fig. 8a), the simulated values of the two peaks (a, b) showed the greatest variability, with R3s.e./mean of 1.3% and 0.8%, respectively, while others were < 0.5%. For Ontario Period 2 (fig. 8b), the largest R3s.e./mean also occurred in the peak or rising hours (c, d, e, f), with R3s.e./mean values of 3.9%, 6.1%, 0.8%, and 4.1% for the four peaks, respectively. The results for the Iowa scenario were similar, where the three simulated peaks most affected (peaks g, h, i) are shown in figure 8c; their R3s.e./mean values were 19.36%, 1.9%, and 1.3%, respectively.
The reason for the effectiveness of macropores during the rising period arises from the fact that preferential flow is faster than soil matrix water movement and reaches the water table first, which contributes to the early arrival of drainage flow. However, this effect is not sensitive to the two macropore parameters, with variations in predicted drainage peaks being within 5% when parameters are varied by 50% to 150%. In addition, after peak time, drainage is dominated by equilibrium flow, so macropores have little impact on the recession period of the drainage hydrograph. Therefore, there were nearly no responses in simulated drainage during the recession to varying macropore parameters.
Conclusions
This study calibrated the macropore component of RZWQM2 against observed hourly tile drainage data from two drainage sites in Ontario and Iowa. At each site, by comparing the performance of model options with and without activating the macropore flow component, we evaluated the overall effects of this component over long-term periods and the model’s ability to simulate drainage peaks on an hourly time scale. Over the simulation periods of high rainfall events, the effects of implementing the macropore component on drainage were insignificant. The drainage distribution remained unchanged except for several drainage peak periods during the long-term simulation. Considering macropores, there was a 6% increase in predicted drainage volume over the four-year period in Ontario and a 3.7% increase over the two-year period in Iowa.
When focusing on drainage peak periods on an hourly scale, the macropore component did improve the model performance in simulating hourly drainage peaks. These peaks, usually underestimated with the NoMP option, generally occurred during an intense rainfall event or under conditions with a saturated soil surface. In such cases, the predicted peak values were closer to the observed values with MP than NoMP. This implied that activating the macropore component would represent actual field conditions more closely. However, in most cases, the total drainage amount and drainage timing predicted by the MP option were worse than those simulated by the NoMP option. A more precise simulation of peak flow would result in a better drainage design, as peak flow is usually used to determine pipe sizes. The preferential flow was dominant in the rising drainage period at the beginning of the rain, causing an initially sharp ascent, similar to observations. In contrast, the curve was flatter during the recession period due to the subsequent equilibrium flow, which did not match the observed hydrograph.
Some previous studies (Jarvis et al., 2008; Larsbo et al., 2014) indicated that preferential flow could also be initiated in unsaturated soils and subsurface soils, which was not simulated in RZWQM2. Although RZWQM2 has the capability to simulate macropore flow initiated on the soil surface, more studies are needed to investigate the hydrologic process of preferential flow in macropore channels and appropriately modify the macropore component of RZWQM2 to provide more precise drainage simulations. Future investigations might include: (1) developing better soil water distributions between soil matrix flow and macropore flow and improving macropore channel flow to shorten drainage peak recession time; and (2) enhancing macropore flow simulation by calibrating the macropore parameters with sub-hourly time resolution drainage data.
References
Ahuja, L.R. (2022). Water and Chemical Transport in the Soil Matrix and Macropores. In Ahuja,L.R., Kersebaum, K.C., & Wendroth, O (Eds), Modeling Processes and Their Interactions in Cropping Systems. Madison, WI: American Society of Agronomy.
Akay, O. & Fox, G.A. (2007). Experimental Investigation of Direct Connectivity between Macropores and Subsurface Drains during Infiltration. Soil Sci. Soc. Am. J., 71: 1600-1606. https://doi.org/10.2136/sssaj2006.0359
Arora, B., Mohanty, B. P., & McGuire, J. T. (2012). Uncertainty in dual permeability model parameters for structured soils. Water Resour. Res., 48(1). https://doi.org/10.1029/2011WR010500
Beven, K., & Germann, P. (2013). Macropores and water flow in soils revisited. Water Resour. Res., 49(6), 3071-3092. https://doi.org/10.1002/wrcr.20156
Brooks, R. H., & Corey, A. T. (1964). Hydraulic properties of porous media and their relation to drainage design. Trans. ASAE, 7(1), 26-28. https://doi.org/10.13031/2013.40684
Ford, W. I., King, K. W., Williams, M. R., & Confesor Jr., R. B. (2017). Modified APEX model for simulating macropore phosphorus contributions to tile drains. J. Environ. Qual., 46(6), 1413-1423. https://doi.org/10.2134/jeq2016.06.0218
Fox, G. A., Malone, R., Sabbagh, G. J., & Rojas, K. (2004). Interrelationship of macropores and subsurface drainage for conservative tracer and pesticide transport. J. Environ. Qual., 33(6), 2281-2289. https://doi.org/10.2134/jeq2004.2281
Fox, G. A., Marvin, M. M., Guzman, J. A., Hoang, C. K., Malone, R. W., Kanwar, R. S., & Shipitalo, M. J. (2012). E. coli transport through surface-connected biopores identified from smoke injection tests. Trans. ASABE, 55(6), 2185-2194. https://doi.org/10.13031/2013.42511
Fox, G. A., Pulijala, S. H., & Sabbagh, G. J. (2007). Influence of rainfall distribution on simulations of atrazine, metolachlor, and isoxaflutole/metabolite transport in subsurface drained fields. J. Agric. Food. Chem., 55(14), 5399-5407. https://doi.org/10.1021/jf063753z
Gao, M., Li, H.-Y., Liu, D., Tang, J., Chen, X., Chen, X.,... Leung, L. R. (2018). Identifying the dominant controls on macropore flow velocity in soils: A meta-analysis. J. Hydrol., 567, 590-604. https://doi.org/10.1016/j.jhydrol.2018.10.044
Hardie, M. A., Cotching, W. E., Doyle, R. B., Holz, G., Lisson, S., & Mattern, K. (2011). Effect of antecedent soil moisture on preferential flow in a texture-contrast soil. J. Hydrol., 398(3), 191-201. https://doi.org/10.1016/j.jhydrol.2010.12.008
Heilman, P., Malone, R. W., Ma, L., Hatfield, J. L., Ahuja, L. R., Ayen, J.,... Kanwar, R. (2006). Decision support for nitrogen management in tile-drained agriculture. Proc. 3rd Int. Congress on Environmental Modelling and Software. USDA-ARS.
Jarvis, N. J., Stähli, M., Bergström, L., & Johnsson, H. (1994). Simulation of dichlorprop and bentazon leaching in soils of contrasting texture using the MACRO model. J. Environ. Sci. Health. Part A: Environ. Sci. Eng. Toxicol., 29(6), 1255-1277. https://doi.org/10.1080/10934529409376106
Jarvis, N., & Larsbo, M. (2012). MACRO (v5.2): Model use, calibration, and validation. Trans. ASABE, 55(4), 1413-1423. https://doi.org/10.13031/2013.42251
Jarvis, N., Etana, A., & Stagnitti, F. (2008). Water repellency, near-saturated infiltration and preferential solute transport in a macroporous clay soil. Geoderma, 143(3), 223-230. https://doi.org/10.1016/j.geoderma.2007.11.015
Jarvis, N., Koestel, J., & Larsbo, M. (2016). Understanding preferential flow in the vadose zone: Recent advances and future prospects. Vadose Zone J., 15(12). https://doi.org/10.2136/vzj2016.09.0075
Jeong, H., & Bhattarai, R. (2018). Exploring the effects of nitrogen fertilization management alternatives on nitrate loss and crop yields in tile-drained fields in Illinois. J. Environ. Manag., 213, 341-352. https://doi.org/10.1016/j.jenvman.2018.02.062
Jiang, Q., Qi, Z., Madramootoo, C. A., & Singh, A. K. (2018). Simulating hydrologic cycle and crop production in a subsurface drained and sub-irrigated field in Southern Quebec using RZWQM2. Comput. Electron. Agric., 146, 31-42. https://doi.org/10.1016/j.compag.2018.01.021
Julich, D., Makowski, V., Feger, K.-H., & Julich, S. (2022). Phosphorus fluxes in two contrasting forest soils along preferential pathways after experimental N and P additions. Biogeochemistry, 157(3), 399-417. https://doi.org/10.1007/s10533-021-00881-w
Klaus, J., Zehe, E., Elsner, M., Külls, C., & McDonnell, J. J. (2013). Macropore flow of old water revisited: Experimental insights from a tile-drained hillslope. Hydrol. Earth Syst. Sci., 17(1), 103-118. https://doi.org/10.5194/hess-17-103-2013
Kumar, A., Kanwar, R. S., & Ahuja, L. R. (1998). Evaluation of preferential flow component of rzwqm in simulatingwater and atrazine transport to subsurface drains. Trans. ASAE, 41(3), 627-637. https://doi.org/10.13031/2013.17231
Larsbo, M., Koestel, J., & Jarvis, N. (2014). Relations between macropore network characteristics and the degree of preferential solute transport. Hydrol. Earth Syst. Sci., 18(12), 5255-5269. https://doi.org/10.5194/hess-18-5255-2014
Lawlor, P. A., Helmers, M. J., Baker, J. L., Melvin, S. W., & Lemke, D. W. (2008). Nitrogen application rate effect on nitrate-nitrogen concentration and loss in subsurface drainage for a corn-soybean rotation. Trans. ASABE, 51(1), 83-94. https://doi.org/10.13031/2013.24229
Liu, H. L., Yang, J. Y., Tan, C. S., Drury, C. F., Reynolds, W. D., Zhang, T. Q.,... Hoogenboom, G. (2011). Simulating water content, crop yield and nitrate-N loss under free and controlled tile drainage with subsurface irrigation using the DSSAT model. Agric. Water Manag., 98(6), 1105-1111. https://doi.org/10.1016/j.agwat.2011.01.017
Liu, M. X., Cui, W. H., Wu, D., Liao, L. J., & Du, W. Z. (2014). Soil Macropore structures and their effect on preferential flow. Appl. Mech. Mater., 522-524, 990-994. https://doi.org/10.4028/www.scientific.net/AMM.522-524.990
Lu, C. (2015). Modeling surface runoff and subsurface tile drainage under drainage and controlled drainage with sub-irrigation in southern Ontario. MS thesis. Montreal, Quebec: McGill University, Department of Bioresource Engineering.
Ma, L., Ahuja, L. R., Nolan, B. T., Malone, R. W., Trout, T. J., & Qi, Z. (2012). Root Zone Water Quality Model (RZWQM2): Model use, calibration, and validation. Trans. ASABE, 55(4), 1425-1446. https://doi.org/10.13031/2013.42252
Malone, R. W., Ma, L., Wauchope, R. D., Ahuja, L. R., Rojas, K. W., Ma, Q.,... Byers, M. (2004a). Modeling hydrology, metribuzin degradation and metribuzin transport in macroporous tilled and no-till silt loam soil using RZWQM. Pest Manag. Sci., 60(3), 253-266. https://doi.org/10.1002/ps.738
Malone, R. W., Shipitalo, M. J., Ma, L., Ahuja, L. R., & Rojas, K. W. (2001). Macropore component assessment of the Root Zone Water Quality Model (RZWQM) using no-till soil blocks. Trans. ASAE, 44(4), 843-852. https://doi.org/10.13031/2013.6249
Malone, R. W., Weatherington-Rice, J., Shipitalo, M. J., Fausey, N., Ma, L., Ahuja, L. R.,... Ma, Q. (2004b). Herbicide leaching as affected by macropore flow and within-storm rainfall intensity variation: A RZWQM simulation. Pest Manag. Sci., 60(3), 277-285. https://doi.org/10.1002/ps.791
Mencaroni, M., Dal Ferro, N., Radcliffe, D. E., & Morari, F. (2021). Preferential solute transport under variably saturated conditions in a silty loam soil: Is the shallow water table a driving factor? J. Hydrol., 602, 126733. https://doi.org/10.1016/j.jhydrol.2021.126733
Moriasi, D. N., Gitau, M. W., Pai, N., & Daggupati, P. (2015). Hydrologic and water quality models: Performance measures and evaluation criteria. Trans. ASABE, 58(6), 1763-1785. https://doi.org/10.13031/trans.58.10715
Nazari, S., Ford, W. I., & King, K. W. (2020). Impacts of preferential flow and agroecosystem management on subsurface particulate phosphorus loadings in tile-drained landscapes. J. Environ. Qual., 49(5), 1370-1383. https://doi.org/10.1002/jeq2.20116
Nimmo, J. R. (2016). Quantitative framework for preferential flow initiation and partitioning. Vadose Zone J., 15(2). https://doi.org/10.2136/vzj2015.05.0079
Orozco-López, E., Muñoz-Carpena, R., Gao, B., & Fox, G. A. (2018). Riparian vadose zone preferential flow: Review of concepts, limitations, and perspectives. Vadose Zone J., 17(1), 180031. https://doi.org/10.2136/vzj2018.02.0031
Qi, Z., Helmers, M. J., Malone, R. W., & Thorp, K. R. (2011). Simulating long-term impacts of winter rye cover crop on hydrologic cycling and nitrogen dynamics for a corn-soybean crop system. Trans. ASABE, 54(5), 1575-1588. https://doi.org/10.13031/2013.39836
Que, Y., Lin, P., & Lin, D. (2018). Integrative analysis of surface runoff and macropore flow for slopes under rainfall conditions. Math. Probl. Eng., pp.1-13. https://doi.org/10.1155/2018/9458410.
Reck, A., Jackisch, C., Hohenbrink, T. L., Schröder, B., Zangerlé, A., & van Schaik, L. (2018). Impact of temporal macropore dynamics on infiltration: Field experiments and model simulations. Vadose Zone J., 17(1), 170147. https://doi.org/10.2136/vzj2017.08.0147
Sadhukhan, D., Qi, Z., Zhang, T., Tan, C. S., Ma, L., & Andales, A. A. (2019). Development and evaluation of a phosphorus (P) module in RZWQM2 for phosphorus management in agricultural fields. Environ. Model. Softw., 113, 48-58. https://doi.org/10.1016/j.envsoft.2018.12.007
Shipitalo, M. J., Nuutinen, V., & Butt, K. R. (2004). Interaction of earthworm burrows and cracks in a clayey, subsurface-drained, soil. Appl. Soil Ecol., 26(3), 209-217. https://doi.org/10.1016/j.apsoil.2004.01.004
Šimunek, J., van Genuchten, M. T., & Šejna, M. (2012). HYDRUS: Model use, calibration, and validation. Trans. ASABE, 55(4), 1263-1274. https://doi.org/10.13031/2013.42239
Tan, C. S., Zhang, T. Q., Drury, C. F., Welacky, T. W., & Reynolds, W. D. (2009). Impacts of liquid and solid manures under controlled drainage with sub-irrigation recycling systems on water quality and crop production. Paper No. 095903. Proc. 2009 ASABE Annual Int. Meeting. St. Joseph, MI: ASABE. https://doi.org/10.13031/2013.26991
Tiktak, A., Hendriks, R. F., & Boesten, J. J. (2012). Simulation of movement of pesticides towards drains with a preferential flow version of PEARL. Pest Manag. Sci., 68(2), 290-302. https://doi.org/10.1002/ps.2262
Weiler, M. (2017). Macropores and preferential flow — a love-hate relationship. Hydrol. Process., 31(1), 15-19. https://doi.org/10.1002/hyp.11074
Zhou, X., Helmers, M., & Qi, Z. (2013). Modeling of subsurface tile drainage using MIKE SHE. Appl. Eng. Agric., 29(6), 865-873. https://doi.org/10.13031/aea.29.9568