Top Navigation Bar

Article Request Page ASABE Journal Article

Quantifying the Importance of Preferential Flow in a Riparian Buffer

L. Guertault, G. A. Fox, T. Halihan, R. Muñoz-Carpena


Published in Transactions of the ASABE 64(3): 937-947 (doi: 10.13031/trans.14286). Copyright 2021 American Society of Agricultural and Biological Engineers.


Submitted for review on 30 August 2020 as manuscript number NRES 14286; approved for publication as a Research Article and as part of the Preferential Flow and Piping in Riparian Buffers Collection by the Natural Resources & Environmental Systems Community of ASABE on 26 January 2021.

The authors are Lucie Guertault, Assistant Professor, and Garey A. Fox, Professor and Head, Department of Agricultural and Biological Engineering, North Carolina State University, Raleigh, North Carolina; Todd Halihan, Professor, Boone Pickens School of Geology, Oklahoma State University, Stillwater, Oklahoma; Rafael Muñoz-Carpena, Professor, Department of Agricultural and Biological Engineering, University of Florida, Gainesville, Florida. Corresponding author: Lucie Guertault, D.S. Weaver Labs 178, Raleigh, NC 27695; phone: 919-515-6754; e-mail: lsguerta@ncsu.edu.

Highlights

Abstract. Riparian buffers are uniquely susceptible to preferential flow due to the abundance of root channels, biological activity, and frequent wetting and drying cycles. Previous research has indicated such susceptibility and even measured the connectivity of preferential flow pathways with adjacent streams and rivers. However, limited research has attempted to partition the riparian buffer infiltration between matrix and preferential flow domains. The objectives of this research were to develop an innovative method to quantify soil matrix infiltration at the plot scale, develop a method to partition infiltration into matrix and macropore infiltration at the plot scale, and then use these methods to quantify the significance of macropore infiltration at a riparian buffer site. This research further demonstrated the importance of considering preferential flow processes in design tools and models to evaluate riparian buffer effectiveness. Sprinkler and runon field experiments were conducted at an established riparian buffer site with sandy loam soil. Trenches were installed and instrumented with soil moisture sensors along the width of the riparian buffer (i.e., along the flow path toward the stream) for detecting non-uniform flow patterns due to preferential flow. Riparian buffer parameters, including soil hydraulic parameters, were estimated using HYDRUS-1D for the sprinkler experiments and VFSMOD for the runon experiments. This research partitioned the infiltration into matrix and preferential flow domains by assuming negligible exchange of water between the soil matrix and preferential flow pathways in comparison to the magnitude of soil matrix flow. For these experimental conditions with 0.20 to 0.48 L s-1 of runon and initial soil water contents of 0.29 to 0.32 cm3 cm-3, preferential flow accounted for at least 27% to 32% of the total runon water entering the riparian buffer. This corresponded to approximately 32% to 47% of the total infiltration. While increasing the riparian buffer plot soil hydraulic conductivity in single-porosity models can adequately predict the total infiltration and therefore the surface outflow from the buffer, design tools and models should specifically consider preferential flow processes to improve predictive power regarding the actual infiltration processes and correspondingly the non-equilibrium flow and solute transport mechanisms.

Keywords. Flow partitioning, HYDRUS, Matrix flow, Preferential flow, Riparian buffer, VFSMOD.

Riparian buffers and vegetative filter strips are commonly used mitigation practices to limit the transport of sediment, nutrients, pesticides, bacteria, and other contaminants from upslope fields into receiving surface water bodies (Momm et al., 2019; Fox et al., 2021). Their effectiveness has been demonstrated to depend not on the physical characteristics of the riparian buffer or vegetative filter strip (i.e., width) but rather on the mechanisms of contaminant detention and retainment (Sabbagh et al., 2009; Yu et al., 2011, 2013; Fox and Penn, 2013; Wu et al., 2014a, 2014b; Muñoz-Carpena et al., 2018b; Cole et al., 2020; Walton et al., 2020). More specifically, the mechanisms of infiltration and sedimentation are key drivers (Muñoz-Carpena et al., 1999). Balestrini et al. (2016) monitored nine riparian buffers in northern Italy and found that the best predictors of the nitrate removal efficiency were factors linked to the water residence time, such as the hydraulic conductivity, the soil texture, and the slope of the riparian profile. Global uncertainty and sensitivity analyses also confirmed that hydraulic conductivity is one of, if not the, most important driving parameters regarding overall trapping efficiency (Muñoz-Carpena et al., 2007, 2010; Fox et al., 2010; Lauvernet and Muñoz-Carpena, 2018). However, when analyzing the performance of riparian buffers or vegetative filter strips, researchers commonly report that the observed hydrologic response of the system suggests a much higher infiltration capacity, and correspondingly a much higher hydraulic conductivity, than expected based on soil texture or matrix infiltration experiments alone (Sabbagh et al., 2009, Heeren et al., 2015).

Riparian buffers, as part of riparian floodplain ecosystems, are uniquely susceptible to preferential flow due to characteristics such as geomorphic depositions, abundant roots and other biological activity, and frequent wetting and drying cycles (Mulholland et al., 1990; Carlyle and Hill, 2001; Polyakov et al., 2005; Fuchs et al., 2009; Orozco-López et al., 2018). The impact of preferential flow on hydrology and contaminant transport within riparian buffers and vegetative filter strips remains an important research question. For example, preferential flow acts to increase the infiltration capacity of the soil. This subsequently reduces the runoff potential, thereby allowing greater sediment deposition, and potentially reducing surface exposure to contaminants with high sorption potential and allowing them to degrade and interact with deeper soil layers. Alternatively, if these preferential flow pathways intersect with high hydraulic conductivity lenses, then the preferential flow could provide a short-circuiting mechanism to adjacent streams and rivers (Fox et al., 2011). This bypass mechanism has been demonstrated explicitly in field studies by Heeren et al. (2010) and more recently by Allaire et al. (2015).

Because of the difficulty in identifying the presence of and then being able to quantify the impact of preferential flow, only a few studies in the literature explicitly highlight preferential flow processes in riparian floodplains. For example, Vellidis et al. (2001) identified that preferential ?ow paths allowed nutrients to bypass a riparian buffer. In conducting large-scale infiltration tests paired with tension infiltrometer experiments in an alluvial floodplain, Heeren et al. (2015) observed that macropores accounted for approximately 84% to 99% of the total saturated hydraulic conductivity.

Riparian buffers and vegetated filter strips can be modeled using tools that couple overland flow with single-porosity subsurface flow processes (Muñoz-Carpena et al., 1999; Abu-Zreig, 2001; Dosskey et al., 2002, 2006; Wu et al., 2014a). When analyzing or designing these systems from a mechanistic or quantitative analysis, simply increasing the soil’s saturated hydraulic conductivity may be sufficient for mimicking the observed surface hydrologic response using a single-porosity approach. However, this approach does not accurately represent non-equilibrium soil hydrologic flow and transport processes, which are hypothesized by numerous researchers to significantly influence contaminant transport and degradation rates (Konopka and Turco, 1991; Fuchs et al., 2009; Muñoz-Carpena et al., 2018a) and can lead to an overestimation of the trapping efficiency of the riparian buffer. Linking overland flow with non-equilibrium subsurface flow is critical to better predict the contaminant removal efficiency of riparian buffer systems. As early as the

late 1990s, Gold and Kellogg (1997) specifically called for the development of unique sampling schemes and simulation models for situations where substantial infiltration occurs through preferential flow pathways. As a more recent example, Phogat et al. (2019) attempted to model the impact of various riparian buffer zone widths for controlling the migration of irrigation water and solutes from wine grapes, almonds, and annual horticulture plots using HYDRUS 2D/3D. They specifically noted the need to incorporate the influence of preferential flow for solute dynamics, but this was not part of their modeling study.

Preferential flow in soils is typically simulated using dual-permeability or dual-porosity models that divide the soil into two conceptual domains, corresponding to the soil matrix and preferential flow, that have different hydraulic properties and conduct flow at different rates (Gerke and van Genuchten, 1993; Šimunek et al., 2003; Larsbo et al., 2005). Dual-permeability models have been successfully applied to simulate agrochemical leaching into field drains in a variety of soil types (e.g., Ludwig et al., 1999; Gerke and Köhne, 2004; Marín-Benito et al., 2014). Source-responsive models have been recently proposed as alternatives to dual-permeability and dual-porosity models (Orozco-López et al., 2018). Advances are currently underway in the mechanistic modeling of riparian buffer systems that are considering the most effective and efficient approach to handle preferential flow (Orozco-López et al., 2018; Fox, 2019; Guertault and Fox, 2020; Orozco-López, 2020; Orozco-López and Muñoz-Carpena, 2021) with the goal of improving riparian buffer and vegetative filter strip design tools. Being able to partition between matrix and preferential infiltration will allow us to elucidate the contributions of matrix and preferential flow in riparian buffers.

Previous studies focusing on the in situ characterization of preferential flow showed that infiltration rates measured at the plot scale (areas of 1 m2 or more) were representative of the entire soil system behavior, combining both matrix flow and preferential flow processes (Heeren et al., 2015), whereas point-scale measurements (length scales below 10 cm) were representative of matrix flow only (Di Prima et al., 2018). Promising advances are also underway to characterize fast and complex macropore soil water dynamics in the laboratory (pore-scale and sub-second time resolution), which are necessary to understand the interaction between the soil matrix and macropores. For example, the use of the light-transmission method in thin two-dimensional soil profiles has proven very efficient in measuring point-scale soil moisture, drainage outflow dynamics, and mass balance at the micrometer and smaller scales (Orozco-López et al., 2021).

The objectives of this research were to develop an innovative method to quantify soil matrix infiltration at the plot scale, develop a method to partition infiltration into matrix and macropore infiltration, and use these methods to quantify the significance of macropore infiltration at a riparian buffer site. The results of this partitioning can then provide initial insights on the magnitude of preferential flow and its impacts on infiltration into and outflow from riparian buffers.

Materials and Methods

Field Site

The field site was an established riparian buffer adjacent to a tributary of Walnut Creek, a first-order stream located on the Lonnie Poole golf course in Raleigh, North Carolina (35° 45' 34.7? N, 78° 40' 46.0? W). This riparian buffer is adjacent to and receives runoff from a fairway of the golf course. As such, this buffer does not receive agricultural runoff; therefore, the functionality and quality of the buffer may differ substantially from a similarly established buffer that has been receiving sediment and nutrient-laden agricultural runoff for the duration of its establishment. The soil was a sandy loam composed of 78% sand, 12% silt, and 10% clay. The average dry bulk density of the soil was 1.5 g cm-3. The soil texture and bulk density were used in a pedotransfer function (Rosetta module within RETC; Schaap et al., 2001) to estimate an expected value of the saturated hydraulic conductivity for the riparian buffer soil. This natural riparian buffer had mature vegetation that included trees, bushes, and shorter plants.

A 1.6 m wide (cross-slope direction) and 8.4 m long (running slope direction) plot was delineated between large trees, with the lower end of the plot located 50 cm from the streambank edge (fig. 1). The width was chosen to match the dimension of the runon source. The plot length was the maximum length for which the plot slope was uniform, as the buffer slope became steeper farther from the stream.

Tall vegetation inside the plot was trimmed to a height of approximately 10 cm to enable covering the plot with a tarp during storm events between the experiments. The plot was covered with a tarp on several occasions throughout the experimental period. Short vegetation, such as grass and twigs, was left alone. Leaf litter was cleared from the plot. Plastic edging buried 3 cm into the soil was installed on the edges of the plot to direct surface flow. A flow-collecting funnel was installed at the lower end of the plot, which discharged into a pipe running through the streambank. The plot topography was surveyed, and the average slope of the plot was 4%.

(a)
(b)
Figure 1. (a) Photo of plot and (b) experimental setup.

Just outside of the plot, three trenches with a minimum depth of 60 cm were dug on each side of the plot at stations located 1.8, 4.0, and 7.2 m from the upper end, hereafter referred to as the upper, middle, and lower trenches, respectively (fig. 1). The exact trench positions were influenced by the presence of large bushes and vegetation around the plot. The goal was to dig a trench near the upper end of the plot but not too close to the source, another trench near the center, and a third trench close to the lower end of the plot. The plot sections centered on the 1.8, 4.0, and 7.2 m stations are referred to as the upper, middle, and lower sections, respectively. In the trenches, EC-5 soil moisture sensors (METER Environment, Pullman, Wash.) were inserted below the plot at depths of 10, 20, 30, and 40 cm (fig. 1). EC-5 sensors have a measurement volume of approximately 0.2 L and therefore provide near point-scale estimates of water content. The sensors were connected to EM-50 data loggers (METER Environment), recording soil water contents at 1 min intervals.

Field Experiments

While tension infiltrometers can determine the saturated hydraulic conductivity of a soil matrix, they are limited to the point scale. To characterize soil water flow mechanisms within the plot on a larger scale, sprinkler experiments were first performed over approximately 3 m × 3 m sections of the plot covering the upper, middle, or lower plot areas (table 1). Water was applied using a single nozzle (FullJet, 1/4GG-12W, Spraying System Co., Charlotte, N.C.) mounted on a tripod. Point-scale soil water contents below each plot section were measured to capture the propagation of the wetting front. For these tests, applied rainfall rates ranged from 42 to 51 mm h-1 and initial water contents were on the drier side (0.11 to 0.18 cm3 cm-3). The sprinkler water application rates were selected such that runoff was not generated in the plots and such that the soil did not become saturated to limit activation of macropore flow. Tests that exhibited sequential soil wetting patterns (i.e., those that develop in an orderly sequence) were assumed to represent uniform water flow processes in the soil matrix.

Following the sprinkler experiments, runon experiments were performed to quantify the partitioning of incoming runon into surface and subsurface flow. A line source was used to apply artificial runon at the upper end of the plot. Tests were performed using two constant inflow rates (0.48 and 0.20 L s-1) at different initial soil moisture conditions (table 2). The inflow rates were selected and tested prior to the experiments to verify that they were sufficient to generate surface outflow from the riparian buffer plot, i.e., the inflow exceeded the infiltration by matrix and preferential flow. Experimental durations were determined during the experiments based on the flow conditions reaching a quasi-steady-state condition. Electrical resistivity imagery (ERI) was also collected during tests R2 and R3, which required the runon experiments to be performed for specific and extended periods of time for data collection (Halihan et al., 2021).

Table 2. Characteristics of runon tests conducted on riparian buffer.
TestDate
(in 2019)
and Time
Duration
(min)
Inflow
Rate
(L s-1)
Initial Surface
Soil Water
Content[a]
(cm3 cm-3)
R17 June
9:02 a.m.
610.480.23
R211 June
4:00 p.m.
1500.480.32
R313 June
10:48 a.m.
2050.200.29

    [a] Measured based on the 10 cm soil moisture sensor.

Table 1. Characteristics of sprinkler tests conducted on riparian buffer plot. Tests S1 and S2 had the same application scenario, but water contents were measured in the left and right trenches, respectively.
TestDate
(in 2019)
and Time
Duration
(min)
Location
in Plot
Initial Surface
Soil Water
Content[a]
(cm3 cm-3)
Rainfall
Rates
(mm h-1)
S128 May
8:44 a.m.
226Lower
left
0.1148 (120 min)
0 (46 min)
51 (60 min)
S228 May
8:44 a.m.
226Lower
right
0.1348 (120 min)
0 (46 min)
51 (60 min)
S34 June
11:44 a.m.
345Upper
right
0.1842

    [a] Measured based on the 10 cm soil moisture sensor.

The frequency associated with each runon scenario was determined to characterize the representativity of each experiment. The drainage area contributing flow to the plot was estimated using GIS to be approximately 350 m2, with a maximum length of 220 m and slope of 7%. The time of concentration of the drainage area was estimated to be 20 min. The rational method was used to estimate which rainfall intensities would produce runon rates of 0.48 and 0.20 L s-1. Using a runoff coefficient (C = 0.15) for the golf course, the inflow rates of 0.48 and 0.2 L s-1 corresponded to rainfall intensities of 33.2 and 13.8 mm h-1, respectively. Using the intensity-duration-frequency table for Raleigh, North Carolina, and assuming a storm duration equal to the time of concentration of the drainage area plus the duration of the runoff experiment, test R1 was representative of a 1-year (1.3 h) storm, test R2 corresponded to a 5-year (3.0 h) storm, and test R3 corresponded to a 1-year (3.75 h) storm. If a rainfall duration of only 1 h was considered, the rainfall intensity of 13.8 mm h-1 that resulted in the flow rate of 0.20 L s-1 (R3) would be observed nine times a year on average.

At the lower end of the plot, overland flow was funneled into a pipe and measured manually by collecting the outflow and using a graduated cylinder to quantify the volumetric discharge (q) from the flume. Using a water budget approach, plot-scale infiltration rates integrating both the matrix and preferential flow pathway contributions were determined as the difference between surface inflow and outflow. During test R3, a Rhodamine WT tracer was added to the artificial runon 162 min after the start of the experiment. The surface flow velocity was estimated for this set of conditions from the appearance of the tracer at the end of the buffer to provide additional qualitative information about the extent of preferential flow in the riparian buffer plot. Between the tests, the plot was covered with a tarp during intense rainfall events.

Estimation of Riparian Buffer Properties

Inverse modeling was used to estimate effective in situ soil hydraulic and riparian buffer properties from the measurements obtained during the sprinkler and runon experiments. Sprinkler tests that exhibited sequential wetting patterns (i.e., developed in an orderly sequence) were selected to calibrate local matrix soil water retention curve parameters for the lower left, lower right, and upper right sections of the plot (table 1). Simulations were performed using the single-porosity model of HYDRUS-1D (Šimunek et al., 2013). The single-porosity model simulates one-dimensional water movement in soils by solving the Richards equation coupled with the van Genuchten soil water retention curve for variably saturated water flow:

(1)

(2)

(3)

where ? is the soil moisture content (cm3 cm-3), K(?) is the unsaturated hydraulic conductivity (mm h-1), Ksat,m is the saturated matrix soil hydraulic conductivity (mm h-1), h is the water pressure head (mm), t is time (h), z is the vertical coordinate (mm), Se is the effective saturation (-), ?r is the residual soil moisture content (-), ?s is the saturated soil moisture content (-), and a (1/mm), n, and m = 1-1/n are empirical coefficients affecting the shape of the hydraulic functions (Šimunek et al., 2013). Model inputs including the depth of the soil profile and initial soil water contents along the profile were specified based on the field site characteristics. The surface boundary condition was set to the sprinkler rate with no ponding, and the bottom boundary condition was set to free drainage. For each test, local matrix soil water retention curve parameters ?r, ?s, a, n, and Ksat,m were assumed homogeneous throughout the vertical soil profile and calibrated using inverse modeling to reproduce the point-scale water contents measured during the tests (table 1). Model performance was evaluated using the root mean square error (RMSE) and the average of the Nash-Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970) between observed (using the sensors exhibiting sequential wetting patterns) and simulated ?:

(4a)

(4b)

where ?obs,i is the ith observation of water content, ?sim,i is the model estimate for the ith observation, and is the average observation. NSE greater than 0.65 has been proposed as a stringent threshold for model acceptability (Ritter and Muñoz Carpena, 2013). Calibrated plot-scale Ksat,m values were then averaged to obtain a representative value for the plot.

Inverse modeling of the runon tests (table 2) was performed to determine soil hydraulic properties representative of the entire plot. Simulations were performed using the riparian buffer hydrologic model VFSMOD (Muñoz Carpena et al., 1999), a field-scale, mechanistic model designed to route a storm hydrograph through a vegetative buffer and calculate surface outflow and infiltration. The model couples a finite-element kinematic wave overland flow submodel with a modified Green-Ampt infiltration submodel:

(5)

(6)

where fp is the infiltration rate (mm s-1), Ksat,p is the saturated vertical hydraulic conductivity for the plot (mm s-1), M is the initial soil-water deficit (cm3 cm-3), Sav is the average suction across the wetting front (mm), Fp is the cumulative infiltration after ponding (mm), F is the cumulative infiltration (mm), tp is the time to ponding (s), and to is the shift of the time scale to correct for not having ponded conditions at the start of the event (s). The model has been validated against field observations in many studies (Kuo and Muñoz-Carpena, 2009; Poletika et al., 2009; Sabbath et al., 2009; Fox and Penn, 2013). The model has been accepted as a reference tool for the study and design of vegetative buffers and is used regularly by state and federal agencies and research institutions (Muñoz-Carpena et al., 2010; Fox et al., 2010; Sabbagh et al., 2013).

Model inputs including filter length, width, slope, initial soil water content (?s), and inflow hydrograph were specified based on the field site characteristics and experimental conditions described in table 2. The ?s was estimated as the average of the maximum soil water contents measured with the soil moisture sensors, with a value of 0.44 cm3 cm-3. For each experiment, the plot Manning’s coefficient (np), Ksat,p, and Sav were calibrated to reproduce the flow observed at the lower end of the plot, more specifically the timing of the arrival of the flow wave and the volume of outflow collected during the experiment. Model performance was assessed by comparing observed and simulated total surface outflow volumes as well as NSE(q) for the surface outflow rates (q). The Ksat,p calibrated for the entire plot system represented the combined behavior of the matrix and preferential flow domains and was larger than Ksat,m.

Partitioning between Surface, Matrix, and Preferential Flow

Observed infiltration as a function of time during the runon experiments was calculated as the difference between runon inflow and outflow from the riparian buffer. This observed infiltration therefore included both matrix and macropore infiltration. To calculate matrix infiltration at the plot scale, additional VFSMOD simulations of the runon experiments were performed using the average Ksat,m estimated for the plot. The infiltration rates computed by the model were assumed to be representative of the contribution of solely the soil matrix. The partitioning of infiltration between matrix flow and preferential flow was then determined by comparing the inflow hydrograph with the surface outflow hydrograph simulated using Ksat,m and the observed surface outflow hydrograph. This approach to partition infiltration between matrix and preferential flow assumed that exchanges between the macropores and the matrix in the upper soil layers were negligible. More specifically, it was assumed that exchanges from the macropores to the matrix that contributed to the wetting of the upper layers of the matrix had a smaller magnitude than the overall matrix flow. This assumption is likely more valid for coarser soils, such as the sandy loam soil at this riparian buffer, as opposed to soils with much higher matric suctions. The uncertainty of this approach combines the uncertainty associated with the measurement of surface inflow and outflow rates, the uncertainty related to the modeling of the matrix flow using VFSMOD, including the calibration of Ksat,m, and the uncertainly related to the assumption of negligible exchanges between the macropores and the matrix.

Results and Discussion

Preferential Subsurface Flow Observations

During runon tests, visual evidence of preferential flow was observed, including flow in earthworm and root tunnels discharging to the trenches. Concentrated subsurface discharge was also observed at several locations on the streambank, approximately 40 cm below the surface. The trenches, including those most separated from the runon source, filled with water without any evidence of entering surface flow, but there was flow from identified macropores. Such observations are consistent with previous studies on preferential flow (Graham and Lin, 2011). During test R3, the Rhodamine WT tracer introduced 162 min after the beginning of the test reached the lower end of the plot 6.0 min later, corresponding to an average overland flow velocity of 2.3 cm s-1 for the specific experimental conditions. The dye tracer was observed on the streambank only 7.3 min after being introduced at the upper end of the plot, demonstrating that at least some preferential pathways were highly conductive and connected to the surface.

Soil Matrix Hydraulic Parameters

Plot-scale Ksat,m from the HYDRUS 1-D inverse modeling was 11 mm h-1 for the lower-left location (test S1), 33 mm h-1 for the lower-right location (test S2), and 30 mm h-1 for the upper-right location (test S3) (table 3). The average value was 25 mm h-1. This value generally agreed with the 35 mm h-1 predicted by the pedotransfer function (Rosetta; Schaap et al., 2001) based on the soil texture and bulk density. Optimized values for ?r, ?s, a, and n are reported in table 4. Model results are shown in figure 2.

Table 3. Soil matrix hydraulic parameters calibrated with HYDRUS.
TestRainfall Rates
(mm h-1)
Ksat,m
(mm h-1)
NSE(?)RMSE(?)
S148 (120 min)
0 (46 min)
51 (60 min)
110.830.019
S248 (120 min)
0 (46 min)
51 (60 min)
330.780.034
S342300.940.036

Table 4. Inverse modeling of soil water retention curve parameters: ranges of possible values and calibrated values.
?r?sa
(1/mm)
nKsat,m
(mm h-1)
Range0.04-0.080.20-0.450.001-0.0151.3-2.510-50
Test S10.050.240.011.911
Test S20.040.380.0091.533
Test S30.040.440.0022.330
(a)
(b)
(c)
Figure 2. Observed and simulated HYDRUS-1D plot-scale water contents: (a) test S1, (b) test S2, and (c) test S3. Blue arrows indicate periods of rainfall.

For test S1, except for the beginning of the experiment at the 10 cm depth, the propagation of the wetting front as well as the recession observed between the two rainfall applications were well reproduced by the Richards equation. The difference between observations and model results for the 10 cm and 30 cm depths can be attributed to different soil properties in the vertical profile. The model correctly reproduced the propagation of the wetting front for test S2 but was not able to reproduce the recession between the two rainfall applications. Overall, NSE(?) were approximately 0.78 to 0.95 for the three test locations (table 3), well above the threshold for an acceptable model and corresponding to a model “pedigree” of good to very good (Ritter and Muñoz Carpena, 2013). For test S1 in the lower left trench, the soil moisture never exceeded 0.25 cm3 cm-3 and the calibrated Ksat,m was relatively small, suggesting a limited activation of macropore flow. For test S2, the soil moisture approached the saturated soil moisture content but never saturated, while during test S3 the soil saturated near the end of the test and could have activated some macropores. Overall, the results suggest that surface tension viscous flow was the dominant subsurface process at the scale of the observations.

Plot-Scale Soil Hydraulic Parameters

From the runon tests, Manning’s coefficient (np) and average suction at the wetting front (Sav) were inversely estimated using VFSMOD and were consistent between tests (table 5). The Sav was consistent within the range of values for sandy loam soils reported by Rawls et al. (1983). Similar values of Ksat,p were obtained for tests R1 and R2, but a lower value was calculated for test R3, suggesting a dependence of Ksat,p with the inflow rate. The NSE(q) values for the surface outflow rates were approximately 0.70 across all three tests (table 5) and once again above the threshold for an acceptable model (Ritter and Muñoz Carpena, 2013).

The inconsistency in the calibrated Ksat,p values (table 5) compared to the Ksat,m calibrated from the sprinkler experiments (table 3) shows that using a single-domain, equilibrium infiltration approach to model infiltration in a riparian buffer requires inflated saturated hydraulic conductivity values to account for the matrix and preferential flow. For tests R1 and R2, the model overestimated surface flow rates (underestimated infiltration rates) at early times, and then toward the end of the test, underestimated surface flow rates (overestimated infiltration rates) (illustrated for test R2 in fig. 3a). There was also a disparity in the prediction of the water distribution on the soil surface (10 cm) in the previous matrix flow infiltration experiments (fig. 2), indicating perhaps a higher density of macropores in the top soil layer that affected the beginning of flow generation from the plot. A Ksat,p value of 70 mm h-1 in the model improved the predictions for the first part of the test (illustrated for test R2 in fig. 3b). In contrast, the model was able to reproduce the outflow hydrograph for test R3 (fig. 3c). Overall, the results showed that the activation of infiltration processes depended on the surface flow rate and initial water content for these experimental conditions with inflow rates equivalent to annual or greater storms on dry initial soil moisture conditions.

Table 5. Plot-scale hydraulic parameters calibrated with VFSMOD.
TestnpSav
(m)
Ksat,p
(mm h-1)
Error on Total
Outflow Volume
(%)
NSE(q)
R10.0120.10610.60.68
R20.0120.10612.60.68
R30.0120.10441.60.70

Partitioning Between Matrix and Preferential Infiltration

For test R1 (inflow rate of 0.48 L s-1), the VFSMOD simulations using Ksat,m from the matrix-only experiments computed larger surface outflow volumes (fig. 4a) and smaller infiltration rates compared to the observations (fig. 4b). Observed infiltration rates were estimated based on the difference between surface inflow and outflow rates. Because surface flow rates were not monitored at intermediate locations within the plot, infiltration rates could only be estimated from the time surface flow reached the end of the plot. Before this time, the entire surface flow had infiltrated in areas smaller than the entire plot, resulting in even greater infiltration rates.

In figure 4a, the difference between the inflow and surface outflow simulated using Ksat,m was assumed to correspond to the matrix flow. The difference between the outflow simulated using Ksat,m and the observed outflow was assumed to be primarily attributable to preferential flow. During test R1, which lasted 1 h, the inflow volume was partitioned into 30% surface flow, 38% matrix flow, and 32% preferential flow.

During the first hour of test R2 (inflow rate of 0.48 L s-1), the inflow volume was partitioned into 36% surface flow, 38% matrix flow, and 27% preferential flow (fig. 5a). Test R2 was performed with a similar inflow rate as test R1, but wetter initial conditions, which led to a slightly lower contribution of preferential flow. For the entire duration of test R2, the contribution of surface flow increased to 44%, matrix flow decreased to 30%, and preferential flow remained at 27%.

During the first hour of test R3 (lower inflow rate of 0.20 L s-1), the inflow volume was partitioned into 80% matrix flow and 20% preferential flow, with no surface flow. Test R3 was performed with a smaller inflow rate compared to tests R1 and R2, and the matrix was able to infiltrate a large portion of this inflow, therefore reducing the contribution of preferential flow. For the entire duration of test R3, the contribution of surface flow increased to 5%, matrix flow decreased to 64%, and preferential flow increased to 31%.

(a)(b)
Figure 5. Partitioning between surface outflow, matrix flow, and preferential flow for (a) test R2 and (b) test R3.

Research Implications for Flow and Transport

For these experimental conditions of runon under annual or greater storm events and initially fairly dry soil moisture conditions, preferential flow was at least 27% to 32% of the total runon volume that entered the plot or 32% to 47% of the total infiltration across these three tests, based on the above assumptions. Although other factors might have contributed to the differences estimated, the approach represents a good first-order approximation of the macropore impact on the hydrology of riparian buffers for the conditions represented in this study. It is likely that the preferential flow was even greater, considering the potential for exchanges between the soil matrix and the macropores in the upper soil layers. Although informative, the specific contribution of preferential flow in the riparian buffer observed in this study depended on the hydrologic conditions experienced by the riparian buffer during an event. Such hydrologic-dependent behavior further supports the need for dynamic, mechanistic-based modeling approaches to quantify the performance of riparian buffers, especially for considering non-equilibrium flow and solute transport processes (Fox, 2019).

As shown by the good overland flow simulations obtained using VFSMOD, increasing the Ksat,m to an effective Ksat,p to mimic the observed flow can capture the cumulative infiltration of the plot. For sediment, Fox and Sabbagh (2009) noted an exponential relationship between sediment trapping percentage (?E) and infiltration (?Q) based on a meta-analysis of previous vegetative filter strip and riparian buffer studies:

(7)

When using this approximating equation and considering only matrix flow, the resulting ?E for this riparian buffer was predicted to be 78%, 70%, and 92% for tests R1, R2, and R3, respectively. The corresponding ?E when lumping both matrix and preferential flow was higher, at 94%, 89%, and 98%, respectively. Sedimentation has important implications for highly sorbing solutes such as pesticides or total phosphorus. For nitrate or dissolved reactive phosphorus, the implications of significant preferential flow are even more important, as rapid transport through the soil vadose zone could lead to direct transport to the adjacent stream. However, it is important to recognize that the physical description of the field processes was not consistent in this case, and this has consequences for the surface/subsurface partitioning of contaminants.

Numerous studies have shown preferential subsurface connectivity to streams and rivers that may control the fate and transport of solutes (Fuchs et al., 2009; Heeren et al., 2010; Fox et al., 2011; Allaire et al., 2015). Such observations, when combined with the magnitude of preferential flow demonstrated in this study, support the need for hydrologic and solute transport preferential flow routines within riparian buffer models. Orozco-López et al. (2018) and Guertault and Fox (2020) reviewed and evaluated several of those approaches for future incorporation into models such as VFSMOD. Such advances will require a trade-off between improving model physical description and increasing parameterization efforts, especially regarding the preferential flow models being considered. Correspondingly, this will also require additional field characterization experiments. However, we believe that the research community should commit to this improved representation of these dynamic riparian buffers (Fox et al., 2021), especially considering the contribution of preferential flow reported in this research.

Conclusions

Previous research has noted the role of macroporosity and preferential flow in influencing hydrology and solute transport and connecting riparian buffers to stream systems. However, limited research to date has partitioned the contribution of matrix and preferential flow in riparian buffers. This research developed an innovative method to quantify soil matrix infiltration at the plot scale, developed a method to partition infiltration into matrix and macropore infiltration at the plot scale, and then used these methods to quantify the significance of macropore infiltration at a riparian buffer site. For this established, sandy loam riparian buffer with runon rates that corresponded to annual or greater storm events, preferential flow accounted for approximately 32% to 47% of the total infiltration depending on the specific hydrologic conditions, such as the initial soil water content and the runon rate relative to the total infiltration capacity. Such hydrologic dependency further supports the need for mechanistic modeling in the design of riparian buffers and in the evaluation of their effectiveness, especially when considering solute transport. While preferential flow may be beneficial to promote infiltration and trap highly sorbing solutes, it may also create the opportunity for some solutes to bypass the trapping mechanism of riparian buffers, leading to significantly reduced residence times for denitrification or the processing of dissolved reactive phosphorus. Future research should use the proposed techniques to partition infiltration between matrix and preferential flow and investigate the influence of preferential flow on non-equilibrium solute transport through a greater number of riparian buffers. Imperatively, mechanisms for considering preferential flow and solute transport must be incorporated into design tools and models.

Acknowledgements

This work was funded by USDA National Institute of Food and Agriculture (NIFA) Project No. 2016-67019-26855. The authors acknowledge Brian Green, Director of Golf Course Maintenance at the Lonnie Poole Golf Course at North Carolina State University, for providing access to the field site. We also acknowledge Riley Dowdy-Green for assisting with the field experiments, and Grace Thomas and Michael Ankner for reviewing an earlier version of the manuscript.

References

Abu-Zreig, M. (2001). Factors affecting sediment trapping in vegetated filter strips: Simulation study using VFSMOD. Hydrol. Proc., 15(8), 1477-1488. https://doi.org/10.1002/hyp.220

Allaire, S. E., Sylvain, C., Lange, S. F., Theriault, G., & Lafrance, P. (2015). Potential efficiency of riparian vegetated buffer strips in intercepting soluble compounds in the presence of subsurface preferential flows. PLoS One, 10(7), e0131840. https://doi.org/10.1371/journal.pone.0131840

Balestrini, R., Sacchi, E., Tidili, D., Delconte, C. A., & Buffagni, A. (2016). Factors affecting agricultural nitrogen removal in riparian strips: Examples from groundwater-dependent ecosystems of the Po valley (northern Italy). Agric. Ecosyst. Environ., 221, 132-144. https://doi.org/10.1016/j.agee.2016.01.034

Carlyle, G. C., & Hill, A. R. (2001). Groundwater phosphate dynamics in a river riparian zone: Effects of hydrologic flowpaths, lithology, and redox chemistry. J. Hydrol., 247(3), 151-168. https://doi.org/10.1016/S0022-1694(01)00375-4

Cole, L. J., Stockan, J., & Helliwell, R. (2020). Managing riparian buffer strips to optimise ecosystem services: A review. Agric. Ecosyst. Environ., 296, article 106891. https://doi.org/10.1016/j.agee.2020.106891

Di Prima, S., Marrosu, R., Lassabatere, L., Angulo-Jaramillo, R., & Pirastru, M. (2018). In situ characterization of preferential flow by combining plot- and point-scale infiltration experiments on a hillslope. J. Hydrol., 563, 633-642. https://doi.org/10.1016/j.jhydrol.2018.06.033

Dosskey, M. G., Helmers M., J., & Eisenhauer, D. E. (2006). An approach for using soil surveys to guide the placement of water quality buffers. J. Soil Water Cons., 61(6), 344-354.

Dosskey, M. G., Helmers, M. J., Eisenhauer, D. E., Franti T., G., & Hoagland. K., D. (2002). Assessment of concentrated flow through riparian buffers. J. Soil Water Cons., 57, 336-343.

Fox, G. A. (2019). Process-based design strengthens the analysis of stream and floodplain systems under a changing climate. Trans. ASABE, 62(6), 1735-1742. https://doi.org/10.13031/trans.13594

Fox, G. A., & Penn, C. J. (2013). Empirical model for quantifying total phosphorus reduction by vegetative filter strips. Trans. ASABE, 56(4), 1461-1469. https://doi.org/http://dx.doi.org/10.13031/trans.56.10133

Fox, G. A., & Sabbagh, G. J. (2009). Comment on “Major factors influencing the efficacy of vegetated buffers on sediment trapping: A review and analysis.” J. Environ. Qual., 38(1), 1-3. https://doi.org/10.2134/jeq2009.0001le

Fox, G. A., Heeren, D. M., Miller, R. B., Mittelstet, A. R., & Storm, D. E. (2011). Flow and transport experiments for a streambank seep originating from a preferential flow pathway. J. Hydrol., 403(3), 360-366. https://doi.org/10.1016/j.jhydrol.2011.04.014

Fox, G. A., Muñoz-Carpena, R., & Sabbagh, G. J. (2010). Influence of flow concentration on parameter importance and prediction uncertainty of pesticide trapping by vegetative filter strips. J. Hydrol., 384(1), 164-173. https://doi.org/10.1016/j.jhydrol.2010.01.020

Fox, G. A., Muñoz-Carpena, R., Brooks, B., & Hall, T. (2021). Advancing surface water pesticide exposure assessments for ecosystem protection. Trans. ASABE, 64(2), 377-387. https://doi.org/10.13031/trans.14225

Fuchs, J. W., Fox, G. A., Storm, D. E., Penn, C. J., & Brown, G. O. (2009). Subsurface transport of phosphorus in riparian floodplains: Influence of preferential flow paths. J. Environ. Qual., 38(2), 473-484. https://doi.org/10.2134/jeq2008.0201

Gerke, H. H., & Köhne, J. M. (2004). Dual-permeability modeling of preferential bromide leaching from a tile-drained glacial till agricultural field. J. Hydrol., 289(1), 239-257. https://doi.org/10.1016/j.jhydrol.2003.11.019

Gerke, H. H., & van Genuchten, M. T. (1993). A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media. Water Resour. Res., 29(2), 305-319. https://doi.org/10.1029/92WR02339

Gold, A. J., & Kellogg, D. Q. (1997). Modelling internal processes of riparian buffer zones. In N. E. Haycock, T. P. Burt, K. W. Goulding, & G. Pinay (Eds.), Buffer zones: Their processes and potential in water protection (pp. 192-207). Harpenden, UK: Quest Environmental.

Graham, C. B., & Lin, H. S. (2011). Controls and frequency of preferential flow occurrence: A 175-event analysis. Vadose Zone J., 10(3), 816-831. https://doi.org/10.2136/vzj2010.0119

Guertault, L., & Fox, G. A. (2020). Performance of preferential flow models in predicting infiltration through a remolded soil with artificial macropores. Vadose Zone J., 19(1), e20055. https://doi.org/10.1002/vzj2.20055

Halihan, T., Hager, J. P., Guertault, L., & Fox G., A. (2021). Detecting macropore fingering using temporal electrical resistivity imaging. Applied Eng. Agric., https://doi.org/10.13031/aea.14294, in press.

Heeren, D. M., Fox, G. A., & Storm, D. E. (2015). Heterogeneity of infiltration rates in alluvial floodplains as measured with a berm infiltration technique. Trans. ASABE, 58(3), 733-745. https://doi.org/10.13031/trans.58.11056

Heeren, D. M., Miller, R. B., Fox, G. A., Storm, D. E., Penn, C. J., & Halihan, T. (2009). Preferential flow path effects on subsurface contaminant transport in alluvial floodplains. ASABE Paper No. 095995. St. Joseph, MI: ASABE. https://doi.org/10.13031/2013.27025

Konopka, A., & Turco, R. (1991). Biodegradation of organic compounds in vadose zone and aquifer sediments. Appl. Environ. Microbiol., 57(8), 2260-2268. https://doi.org/10.1128/AEM.57.8.2260-2268.1991

Kuo, Y. M., & Muñoz-Carpena, R. (2009). Simplified modeling of phosphorus removal by vegetative filter strips to control runoff pollution from phosphate mining areas. J. Hydrol., 378(3), 343-354. https://doi.org/10.1016/j.jhydrol.2009.09.039

Larsbo, M., Roulier, S., Stenemo, F., Kasteel, R., & Jarvis, N. (2005). An improved dual-permeability model of water flow and solute transport in the vadose zone. Vadose Zone J., 4(2), 398-406. https://doi.org/10.2136/vzj2004.0137

Lauvernet, C., & Muñoz-Carpena, R. (2018). Shallow water table effects on water, sediment, and pesticide transport in vegetative filter strips: Part 2. Model coupling, application, factor importance, and uncertainty. Hydrol. Earth Syst. Sci., 22(1), 71-87. https://doi.org/10.5194/hess-22-71-2018

Ludwig, R., Gerke, H. H., & Wendroth, O. (1999). Describing water flow in macroporous field soils using the modified macro model. J. Hydrol., 215(1), 135-152. https://doi.org/10.1016/S0022-1694(98)00266-2

Marin-Benito, J. M., Pot, V., Alletto, L., Mamy, L., Bedos, C., Barriuso, E., & Benoit, P. (2014). Comparison of three pesticide fate models with respect to the leaching of two herbicides under field conditions in an irrigated maize cropping system. Sci. Total Environ., 499, 533-545. https://doi.org/10.1016/j.scitotenv.2014.06.143

Momm, H. G., Yasarer, L. M., Bingner, R. L., Wells, R. R., & Kunhle, R. A. (2019). Evaluation of sediment load reduction by natural riparian vegetation in the Goodwin Creek watershed. Trans. ASABE, 62(5), 1325-1342. https://doi.org/10.13031/trans.13492

Mulholland, P. J., Wilson, G. V., & Jardine, P. M. (1990). Hydrogeochemical response of a forested watershed to storms: Effects of preferential flow along shallow and deep pathways. Water Resour. Res., 26(12), 3021-3036. https://doi.org/10.1029/WR026i012p03021

Muñoz-Carpena, R., Fox, G. A., & Sabbagh, G. J. (2010). Parameter importance and uncertainty in predicting runoff pesticide reduction with filter strips. J. Environ. Qual., 39(2), 630-641. https://doi.org/10.2134/jeq2009.0300

Muñoz-Carpena, R., Fox, G. A., Ritter, A., Perez-Ovilla, O., & Rodea-Palomares, I. (2018a). Effect of vegetative filter strip pesticide residue degradation assumptions for environmental exposure assessments. Sci. Total Environ., 619-620, 977-987. https://doi.org/10.1016/j.scitotenv.2017.11.093

Muñoz-Carpena, R., Lauvernet, C., & Carluer, N. (2018b). Shallow water table effects on water, sediment, and pesticide transport in vegetative filter strips: Part 1. Nonuniform infiltration and soil water redistribution. Hydrol. Earth Syst. Sci., 22(1), 53-70. https://doi.org/10.5194/hess-22-53-2018

Muñoz-Carpena, R., Parsons, J. E., & Gilliam, J. W. (1999). Modeling hydrology and sediment transport in vegetative filter strips. J. Hydrol., 214(1), 111-129. https://doi.org/10.1016/S0022-1694(98)00272-8

Muñoz-Carpena, R., Zajac, Z., & Kuo, Y. M. (2007). Global sensitivity and uncertainty analyses of the water quality model VFSMOD-W. Trans. ASABE, 50(5), 1719-1732. https://doi.org/10.13031/2013.23967

Nash, J. E., & Sutcliffe, J. V. (1970). River flow forecasting through conceptual models: Part I. A discussion of principles. J. Hydrol., 10(3), 282-290. https://doi.org/10.1016/0022-1694(70)90255-6

Orozco-López, E., & Muñoz-Carpena, R. (2021). Comparative non-Darcian modeling of subsurface preferential flow: Experimental observations in a riparaian buffer. Trans. ASABE, https://doi.org/10.13031/trans.14559, in review.

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), article 180031. https://doi.org/10.2136/vzj2018.02.0031

Orozco-López, E., Muñoz-Carpena, R., Gao, B., & Fox, G. A. (2021). High-resolution pore-scale water content measurement in a translucent soil profile from light transmission. Trans. ASABE, https://doi.org/10.13031/trans.14292, in press.

Orozco-López., E. (2020). Subsurface preferential flow and transport in riparian buffers. PhD diss. Gainesville, FL: University of Florida.

Phogat, V., Cox, J. W., Kookana, R. S., Šimunek, J., Pitt, T., & Fleming, N. (2019). Optimizing the riparian zone width near a river for controlling lateral migration of irrigation water and solutes. J. Hydrol., 570, 637-646. https://doi.org/10.1016/j.jhydrol.2019.01.026

Poletika, N. N., Coody, P. N., Fox, G. A., Sabbagh, G. J., Dolder, S. C., & White, J. (2009). Chlorpyrifos and atrazine removal from runoff by vegetated filter strips: Experiments and predictive modeling. J. Environ. Qual., 38(3), 1042-1052. https://doi.org/10.2134/jeq2008.0404

Polyakov, V., Fares, A., & Ryder, M. H. (2005). Precision riparian buffers for the control of nonpoint-source pollutant loading into surface water: A review. Environ. Rev., 13(3), 129-144. https://doi.org/10.1139/a05-010

Rawls, W. J., Brakensiek, D. L., & Miller, N. (1983). Green-Ampt infiltration parameters from soil data. J. Hydraul. Eng., 109(1), 62-70. https://doi.org/10.1061/(ASCE)0733-9429(1983)109:1(62)

Ritter, A., & Muñoz-Carpena, R. (2013). Performance evaluation of hydrological models: Statistical significance for reducing subjectivity in goodness-of-fit assessments. J. Hydrol., 480, 33-45. https://doi.org/10.1016/j.jhydrol.2012.12.004

Sabbagh, G. J., Fox, G. A., Kamanzi, A., Roepke, B., & Tang, J.-Z. (2009). Effectiveness of vegetative filter strips in reducing pesticide loading: Quantifying pesticide trapping efficiency. J. Environ. Qual., 38(2), 762-771. https://doi.org/10.2134/jeq2008.0266

Sabbagh, G. J., Muñoz-Carpena, R., & Fox, G. A. (2013). Distinct influence of filter strips on acute and chronic pesticide aquatic environmental exposure assessments across U.S. EPA scenarios. Chemosphere, 90(2), 195-202. https://doi.org/10.1016/j.chemosphere.2012.06.034

Schaap, M. G., Leij, F. J., & van Genuchten, M. T. (2001). Rosetta: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol., 251(3), 163-176. https://doi.org/10.1016/S0022-1694(01)00466-8

Šimunek, J., Jarvis, N. J., van Genuchten, M. T., & Gardenas, A. (2003). Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol., 272(1), 14-35. https://doi.org/10.1016/S0022-1694(02)00252-4

Šimunek, J., Šejna, M., Saito, H., Sakai, M., & van Genuchten, M. T. (2013). The Hydrus-1D software package for simulating the movement of water, heat, and multiple solutes in variably saturated media, Ver. 4.17, HYDRUS software series 3. Riverside, CA: University of California, Department of Environmental Sciences.

Vellidis, G., Lowrance, R., Hubbard, R. K., & Gay, P. (2001). Preferential flow caused by past disturbance in a restored riparian wetland. In D. D. Bosch & K. W. King (Eds.), Proc. 2nd Intl. Symp Preferential Flow, Water Movement, and Chemical Transport in the Environment (pp. 61-64). St. Joseph, MI: ASAE. https://doi.org/10.13031/2013.2141

Walton, C. R., Zak, D., Audet, J., Petersen, R. J., Lange, J., Oehmke, C., ... Hoffmann, C. C. (2020). Wetland buffer zones for nitrogen and phosphorus retention: Impacts of soil type, hydrology, and vegetation. Sci. Total Environ., 727, 138709. https://doi.org/10.1016/j.scitotenv.2020.138709

Wu, L., Gao, B., Tian, Y., & Muñoz-Carpena, R. (2014a). Analytical and experimental analysis of solute transport in heterogeneous porous media. J. Environ. Sci. Health A, 49(3), 338-343. https://doi.org/10.1080/10934529.2014.846686

Wu, L., Muñoz-Carpena, R., Gao, B., Yang, W., & Pachepsky, Y. A. (2014b). Colloid filtration in surface dense vegetation: Experimental results and theoretical predictions. Environ. Sci. Tech., 48(7), 3883-3890. https://doi.org/10.1021/es404603g

Yu, C., Gao, B., Muñoz-Carpena, R., Tian, Y., Wu, L., & Perez-Ovilla, O. (2011). A laboratory study of colloid and solute transport in surface runoff on saturated soil. J. Hydrol., 402(1), 159-164. https://doi.org/10.1016/j.jhydrol.2011.03.011

Yu, C., Muñoz-Carpena, R., Gao, B., & Perez-Ovilla, O. (2013). Effects of ionic strength, particle size, flow rate, and vegetation type on colloid transport through a dense vegetation saturated soil system: Experiments and modeling. J. Hydrol., 499, 316-323. https://doi.org/10.1016/j.jhydrol.2013.07.004