Preferential flow in the soil and vadose zone refers to the channeling of infiltrating water and solutes through a small fraction of the pore space (Hendrickx and Flury, 2001). Preferential flow significantly affects the hydrological response to precipitation and therefore both groundwater and surface water quality at all scales of relevance for management and decision making (Jarvis, 2007; Clothier et al., 2008; Beven and Germann, 2013). Preferential flow results from significant spatial variations in water velocities due to heterogeneities in soil properties at both pore and Darcy scales. Preferential flow can also occur in highly uniform soils in the form of unstable “fingers,” especially in coarse, single-grained, and water-repellent soils (DiCarlo, 2013; Wallach et al., 2013). Although the underlying causes of preferential flow may differ, one common signature is a convergence of flow toward preferential pathways and nonequilibrium in terms of differences in water pressures and solute concentrations transverse to the flow direction, which means that no single value can describe the state of the system. Under saturated conditions, preferential solute transport occurs without nonequilibrium in water pressures. However, this update focuses on the unsaturated zone, where preferential flow and transport generally occur simultaneously. We therefore use the term preferential flow throughout to imply nonequilibrium conditions for both water pressures and solute concentrations.
Convergence to preferential flow paths can take place at the soil surface under near-ponding conditions or within the soil if the flow encounters less permeable layers (e.g., Kulli et al., 2003; Bogner et al., 2013a; Li et al., 2013) or even above the soil surface due to the redistribution of precipitation by the plant canopy architecture (e.g., Johnson and Lehmann, 2006; Schwärzel et al., 2012; Li et al., 2009, 2013). Not all kinds of heterogeneity promote nonequilibrium (i.e., preferential) flow. The importance of the spatial correlation structure was demonstrated by Schlüter et al. (2012a), who concluded from numerical experiments that “increased isotropic, short-range connectivity reduces non-equilibrium, whereas anisotropic structures that are elongated in the direction of flow enforce it.” This blueprint for preferential flow is certainly matched by large structural pores (termed macropores) that are almost ubiquitous in soil (Jarvis, 2007). Apart from macropore flow, funnel flow or focused recharge at the Darcy scale can also be induced by much larger scale structures of elongated geometry, such as tongue-like soil horizon interfaces (e.g., Kulasekera and Parkin, 2011; Hardie et al., 2011) or inclined layers or lenses of contrasting permeability in the soil or vadose zone (e.g., Kung, 1990; Heilig et al., 2003; Barkle et al., 2014; Sohrt et al., 2014).
Initially stimulated by the pioneering studies of Johan Bouma, Peter Germann, and Keith Beven in the late 1970s and early 1980s, preferential flow has been an important area of research for more than 30 yr, and interest in the subject is certainly not declining (Gerke et al., 2010). In this update, we identify and discuss some of the main themes that have emerged in recent years, review some significant advances, and suggest some promising areas for future work. For convenience, we first discuss empirically grounded research based on experimentation and then the current status and prospects of simulation models, although we recognize that progress in understanding is best served when inductive and deductive approaches go hand-in-hand. Because preferential flow is such a broad subject, this short review is necessarily selective rather than comprehensive.
Experimentation from Pore to Catchment Scales
In their review of experimental methods, Allaire et al. (2009) identified several priority areas for future research. In the following, we focus on three of these topics where we believe significant progress has been made in the intervening years, namely the development of methods to study macropores and preferential flow at the pore scale, the use of mathematical indicators to quantify preferential flow, and the application of novel methods to study preferential flow at the larger field, hillslope, and catchment scales.
Methods for Studying Preferential Flow at the Pore to Core Scales
Noninvasive X-ray imaging is now increasingly used to quantify the properties of soil macropores in three dimensions. Studies using this technique have shown that soil macroporosity generally comprises a large number of highly complex networks (or “clusters”) separated from one another by smaller pores, each network being characterized by multiple pathways of variable thickness, including bottleneck constrictions and “dead ends” (e.g., Perret et al., 1999, 2000; Pierret et al., 2002; Luo et al., 2010a). This work has led to a growing recognition of the need for appropriate methods that can capture the key topological and geometrical properties of these complex multiscale macropore networks, especially their connectivity and continuity. Vogel et al. (2010) showed how the Minkowski functionals calculated from X-ray images can be used to quantify the variation of porosity, surface area, curvature, and the Euler–Poincaré characteristic (a measure of connectivity) with pore size (thickness). The utility of the mathematical theory of networks (e.g., Newman, 2003) to describe the topology of soil pore space imaged by X-ray has also recently been explored (e.g., Santiago et al., 2008; Mooney and Korošak, 2009; Samec et al., 2013). Based on X-ray images of the soil pore space, these researchers developed pore network models characterized by “scale-free” power law distributions of coordination number (the number of pores connected at a node) and pore thickness. With appropriate parameterization, these models can also produce networks with the “small-world” properties (i.e., the presence of “short-circuiting” long-range connections; Watts and Strogatz, 1998; Latora and Marchiori, 2001; Valentini et al., 2007) that are most likely to sustain preferential flow in macroporous soil. It has also often been suggested that concepts from the closely related topic of percolation theory might prove useful to quantify key aspects of the long-range connectivity (i.e., continuity) of macropore networks (e.g., Nieber et al., 2006; Graham and Lin, 2012). This idea has so far hardly been explored, although Jarvis et al. (2016) recently showed that the structural pore networks in cultivated topsoil imaged by X-ray at a resolution of 65 μm displayed the key characteristic features predicted by percolation theory.
Various noninvasive imaging techniques, including X-ray tomography, neutron radiography, magnetic resonance imaging (MRI) and positron emission tomography (PET) and single-photon emission computed tomography (SPECT) scanning have the potential to reveal new insights into preferential flow by directly imaging the spatial arrangements of water, solutes, and air in the pore space during transport experiments (e.g., Luo et al., 2008; Koestel and Larsbo, 2014; Sammartino et al., 2015). Time-lapse X-ray studies have shown that isolated macropores that are disconnected from the rest of the large pore network by small pore necks or throats do not contribute to fast flow (Perret et al., 2000; Sammartino et al., 2015). Using MRI, Sněhota et al. (2010) showed that the thicker parts of the pore network may also be nonconductive due to air entrapment during infiltration. Luo et al. (2008) and Luo and Lin (2009) showed that significant air volumes are entrapped in the large pores even during slow saturation of the soil from below. This explains why generally only a small fraction of the soil macroporosity conducts flow even under high irrigation rates or surface ponding (see Jarvis, 2007, and references therein; Luo et al., 2008; Sammartino et al., 2012, 2015). The results of Sammartino et al. (2015) suggested that a single dominant (percolating) network comprised most of this conducting or “active” fraction of the macroporosity. Under more typical boundary conditions, water flow in large macropores in the unsaturated vadose zone is almost always unsaturated, simply because their flow capacity is likely to far exceed the supply rate (e.g., Jarvis, 2007; Cey and Rudolph, 2009). The spatial configuration of the air and flowing water in these unsaturated macropore networks should play an important role in regulating preferential flow by directly influencing flow velocities and, indirectly, by determining the wetted area for water and solute imbibition into textural pores. Although the technical challenges would not be insignificant, data obtained from direct imaging of flow and transport experiments may help to shed light on these dynamic processes, which are still not fully understood (Jarvis, 2007; Nimmo, 2007, 2010; Beven and Germann, 2013).
Preferential flow is influenced not only by the properties of the pore space but also by the wettability of the solid surfaces in soil. Badorreck et al. (2012) combined X-ray scanning and neutron radiography to investigate water flow through coarse sediments containing beetle burrows. In an initially dry soil sample, water repellency caused a pressure buildup, leading to gravity-driven flow through the burrow, whereas the burrow did not conduct any flow under initially moist soil conditions. Internal macropore surfaces may often be hydrophobic, which can reduce the rate of imbibition of water into smaller pores and further enhance nonequilibrium flow (see Jarvis, 2007, and references therein; Lipiec et al., 2015; Leue et al., 2015). Recent work has demonstrated that the water repellency of macropore surfaces is largely determined by the chemical composition of organic matter coatings and linings, especially the content of alkyl functional groups (e.g., Leue et al., 2015).
Mathematical Indicators to Quantify Preferential Flow at Core to Pedon Scales
The strength or degree of preferential flow can be quantified by summarizing data on fluxes or contents of water or solute with various mathematical indices and indicators. The dilution index proposed by Kitanidis (1994) is a fundamental entropy-based measure of the degree of preferential flow calculated from the spatial distribution of solute concentrations (or water contents) within a control volume. This index has only rarely been applied to real experiments because these kinds of detailed measurements have not been available. However, as noted above, methods for real-time noninvasive imaging of column flow and transport experiments are developing rapidly and could, in principle, supply the necessary data. For example, Koestel and Larsbo (2014) calculated the temporal evolution of the dilution index from time-lapsed three-dimensional X-ray images of tracer transport through a saturated soil sample containing two large macropores.
Most widely used measurement techniques capture only either spatial (e.g., dye staining) or temporal variations in water and solute transport (e.g., tracer breakthrough curves, BTCs). In the past, BTC data were usually fitted with transport models such as the mobile–immobile version of the convection–dispersion equation (CDE), but simpler model-independent measures that describe the shape of the BTC (i.e., the distribution of solute travel times) are now more often preferred because they make no a priori assumptions about flow mechanisms. In recent years, various BTC shape measures, including temporal moments, the so-called “holdback factor,” and the normalized 5% arrival time (Koestel et al., 2011) have been used to assess how preferential flow is affected by land use and soil management practices (e.g., Mossadeghi-Björklund et al., 2016) or soil properties such as texture (Ghafoor et al., 2013), bulk density (Koestel et al., 2013; Soares et al., 2015; Katuwal et al., 2015a), and organic C content (Larsbo et al., 2016).
Several recent studies have related indicators of the strength of preferential flow derived from solute breakthrough experiments to the pore structures quantified by X-ray imaging on the same samples. Luo et al. (2010b) fitted transport models to tracer breakthrough data obtained from experiments performed under saturated flow conditions and found that the CDE fitted the BTCs for all columns except one that contained a single large continuous (i.e., percolating) macropore. Luo et al. (2010b) also found smaller dispersivities (i.e., more uniform flow) in columns with larger imaged macroporosity. Other recent studies performed under near-saturated conditions also found that preferential flow at a given water flow rate was weaker in columns with larger macroporosity (e.g., Larsbo et al., 2014; Katuwal et al., 2015b; Paradelo et al., 2013, 2016). Larsbo et al. (2014) explained this finding qualitatively in terms of a multiscale percolating pore network: the columns of smaller imaged porosity had smaller near-saturated hydraulic conductivities, which presumably resulted in increased water pressures that were sufficient to generate flow in larger (but still percolating) pores. This reasoning was also used by Mossadeghi-Björklund et al. (2016) to explain why preferential flow may be exacerbated by reductions in macroporosity induced by compaction (see Etana et al., 2013), at least until the percolation threshold for pores large enough to sustain nonequilibrium flow is reached, at which point preferential flow should be significantly curtailed. Column studies performed under near-saturated conditions have also emphasized the role of well-connected dense networks of smaller pores in preventing pronounced preferential flow from being “triggered” in larger macropores. For example, Larsbo et al. (2016) found that initial tracer breakthrough was faster in columns with smaller volumes of pores between 0.2 and 0.6 mm in diameter but was not correlated with the properties of macropores larger than 1.2 mm in diameter. Other studies found that the strength of preferential flow was significantly and positively correlated with the average gray-scale value of the region segmented as matrix, which should reflect the volume of pores smaller than the X-ray resolution, but not with the imaged macroporosity (Soares et al., 2015; Katuwal et al., 2015a; Paradelo et al., 2016).
Shape measures derived from tracer BTCs can serve as indicators of the strength of preferential flow, but they give no information on spatial patterns and thus the dominant flow mechanisms (i.e., fingering or macropore flow). Multi-compartment samplers provide additional data on the spatial distributions of water and solute fluxes for a two-dimensional surface orthogonal to the main flow direction as well as BTCs for each individual compartment (Bloem et al., 2010). Thus, spatiotemporal “leaching surfaces” and statistical indices calculated from these data can be used to both quantify the strength of preferential flow as well as give some indication of the prevailing preferential flow process (e.g., Schotanus et al., 2012). In contrast, dye staining experiments only give two-dimensional “snapshots” of spatial variations in water flow patterns but at the larger pedon scale (Flury et al., 1994). Bogner et al. (2013b) recently summarized several statistical indices that can be calculated from this information to characterize preferential flow, such as the stained area, the widths of stained features, fractal dimensions, Shannon information entropy, and measures of fragmentation. They used statistical clustering techniques to delineate regions in the soil profile with similar flow patterns, which can then be related to variations in soil properties and soil profile morphology. Dye staining experiments are time consuming and destructive, which means that replication is often quite limited. Nevertheless, treatment effects can sometimes be distinguished despite spatial variability in soil properties. For example, in one recent study, Gimbel et al. (2016) used dye staining to investigate the probable effects of climate change on flow processes in forest soils. They found that the strength of preferential flow was enhanced by the development of significant water repellency induced by an imposed 2-yr drought. Aside from these quantitative aspects, simple qualitative observations of the preferential flow processes operating at a site are also of considerable pedagogic value. In this respect, several recent dye staining studies have emphasized the diversity of preferential flow processes that potentially can occur in any one soil. For example, both Kramers et al. (2009) and Etana et al. (2013) observed preferential flow through both macropores and coarse sand lenses in loamy glacial till soils. Another good example is the texture-contrast soils that are common in Australia, which were found to exhibit fingering due to water repellency in the coarse-textured topsoil, macropore flow in the strongly aggregated clayey subsoil, as well as funnel flow at the irregular horizon interface (Hardie et al., 2011, 2013).
Methods for Studying Preferential Flow at Field, Catchment, and Regional Scales
One major challenge is to bridge the gap between the scales at which preferential flow takes place (i.e., single pores to the Darcy scale) and the larger scales relevant for management (i.e., hillslopes, fields, catchments). In some extreme cases, individual macropores can themselves extend across large distances. For example, Wilson et al. (2016) performed BTC experiments with fluorescein dye to demonstrate the connectivity of hillslope pipe networks across distances of more than 190 m from the upper ridges of the catchment to the streams. These “self-organized” networks of large pipes form from the decay or combustion of large tree roots (Leslie and Heinse, 2013) or the internal erosion of preexisting macropores (Wilson et al., 2013) and are commonly found in near-saturated and saturated zones overlying impermeable soil horizons or rocks on steep slopes under natural vegetation (e.g., Sidle et al., 2001; Weiler and McDonnell, 2007; Anderson et al., 2009).
In recent years, geophysical techniques such as ground-penetrating radar (GPR) and electrical resistivity tomography (ERT) have been increasingly applied to investigate flow and transport processes in heterogeneous soils (Vereecken et al., 2014; Binley et al., 2015). These techniques are best suited to studying Darcy-scale flow processes because the spatial resolution is usually not sufficient to image pore-scale preferential flow (e.g., Oberdörster et al., 2010), although this may sometimes be possible for soils characterized by large and sparsely distributed macropores. For example, Greve et al. (2011) used ERT to illustrate the contrasting wetting patterns that were found following the irrigation of dry or moist cracking clay soils, while both ERT (Leslie and Heinse, 2013) and GPR (Holden, 2004; Guo et al., 2014) have also been used to identify functioning soil pipes. The occurrence of preferential flow at large scales can also be assessed from measurements of spatial and temporal variations in soil moisture or solute concentrations obtained with more traditional methods such as time-domain reflectometry (TDR), moisture capacitance probes, tensiometers, suction cups, and soil coring. For example, Yang et al. (2013, 2014) investigated the controls of land use, topography, soil properties, and rainfall characteristics on preferential flow using a novel scale-dependent distribution of treatments within a 144-m2 field plot. Wessolek et al. (2009) combined soil moisture sampling (TDR and gravimetric) with measurements of water drop penetration times and a double tracer experiment to show how preferential finger flow patterns in a 150- by 25-m plot changed seasonally with the degree of water repellency induced by soil wetting and drying cycles. Rye and Smettem (2015) also used intensive soil moisture monitoring to identify seasonal variations in the cross-sectional area of flow through a water-repellent soil at the plot scale.
Extensive arrays of automated (wireless) soil moisture sensors have also been used at the much larger catchment scale (Vereecken et al., 2014). For example, from multiyear data obtained from 144 moisture sensors distributed across a small catchment, Liu and Lin (2015) showed that the occurrence of preferential flow was controlled by topography and soil type, being more prevalent at drier hilltop sites and wetter valley bottoms. In contrast, Wiekenkamp et al. (2016) could not identify any systematic effects of soil type or topography on preferential flow from a similar network of sensors distributed across a 38-ha forested catchment in Germany. At that site, preferential flow was predominantly controlled by precipitation amounts and antecedent soil moisture, such that the degree of preferential flow was enhanced by large storms falling on dry soil, presumably due to water repellency.
Ghafoor et al. (2013) measured tracer BTCs on 45 topsoil columns sampled from conventionally tilled fields in an agricultural catchment 16 km2 in size. Under these conditions of similar management and climate, preferential flow was controlled by variations in basic soil properties, mainly soil texture but also organic C content. At the larger landscape scale, all the factors of soil formation (climate, parent material, topography, vegetation and fauna, land use, and management) can potentially influence preferential flow via processes of soil structural development and degradation (Lin, 2003; Jarvis et al., 2012). Some attempts have been made in recent years to disentangle the complexities of the factors that control preferential flow at the landscape scale. For example, Jarvis et al. (2009) developed a simple decision tree to classify susceptibility to macropore flow from basic soil properties and site attributes, which they then tested against a small dataset of BTCs taken from the literature. Later, Koestel and Jorda (2014) took a more inductive approach and tested whether the application of machine learning techniques to the global meta-database of literature BTC studies collated by Koestel et al. (2012) could lead to new insights into how site attributes, soil properties, and experimental conditions might influence preferential flow. This large dataset (591 data entries) still suffered from some significant correlations between predictor variables as well as data gaps. Nevertheless, analyses with random forests did identify some key factors controlling macropore flow, especially the clay content above a threshold value of approximately 8 to 10%, which was thought to correspond to the value required to initiate aggregation (Koestel and Jorda, 2014).
At regional scales, spatial variation in the climate, especially in precipitation regimes, should strongly influence the occurrence of preferential flow and leaching risks, especially for adsorbing and degrading contaminants such as pesticides (e.g., Nolan et al., 2008; McGrath et al., 2010). These relationships have been explored using models supported by empirical pedotransfer functions to estimate the parameters. For example, Iversen et al. (2011) mapped the vulnerability to macropore flow in Denmark using pedotransfer functions for near-saturated hydraulic conductivity developed from soil survey data, coupled with simulation modeling to account for the effects of climate. Steffens et al. (2015) used a simulation tool based on a macropore flow model parameterized by pedotransfer functions to assess the likely effects of direct and indirect effects of climate change on herbicide leaching risks to groundwater at the regional scale. McGrath et al. (2009) developed a simpler preferential flow index that could be used at large spatial scales to predict site vulnerability to pesticide leaching depending on soil properties and rainfall characteristics.
Model Concepts and Simulation
Heterogeneous Flow and Fingering at the Darcy Scale
The effects of Darcy-scale heterogeneity can be simulated by appropriate parameterizations of Richards’ equation for water flow and the CDE for solute transport in two- or three-dimensional simulations (e.g., Kasteel et al., 2000; Javaux et al., 2006; Schlüter et al., 2012b). Preferential flow then becomes an emergent property of the smaller scale heterogeneities (Vogel and Roth, 2003). One challenge with these physics-based approaches is model parameterization, which may be one reason why one-dimensional empirical models (e.g., Ritsema et al., 2005; Liu et al., 2005; Sheng et al., 2009) are still often used. However, geophysical techniques such as GPR and ERT could be utilized to support the parameterization of spatially explicit two- and three-dimensional models at plot to field scales (Vereecken et al., 2014). For even larger scales (e.g., hillslopes and catchments), a lack of computational power may lead to a spatial resolution that is too coarse for accurate solutions (e.g., Vogel and Ippisch, 2008; Or et al., 2015). In the future, these limitations may be alleviated by advances in numerical techniques and computational power (Orgogozo et al., 2014), but model parameterization will surely remain a major challenge.
Deurer and Bachmann (2007) showed that heterogeneous flow patterns could be simulated for water-repellent soils using Richards’ equation in two dimensions by assuming a dynamic contact angle that varies spatially and temporally as a function of soil water status. However, this approach cannot simulate unstable finger flow, which also commonly occurs in water-repellent soil, because Richards’ equation cannot predict the “saturation overshoots” (i.e., the reversal of saturation and pressure potential gradients at the wetting front) that actually cause the flow instabilities (DiCarlo, 2010, 2013). Several different empirical modifications to Richards’ equation have been suggested to enable simulation of dynamic saturation overshoots and unstable finger flow, but they do not seem to match experimental observations of pressure and saturation during fingering (DiCarlo, 2010). Instead, recent advances in the understanding of the physics of flow instability suggest that although unstable “fingering” manifests itself at the Darcy scale, the process can only be properly understood and modeled as a pore-scale wetting process (DiCarlo, 2013; Baver et al., 2014).
Models that can account for nonequilibrium flow at the pore scale can be classified according to whether macropores are described implicitly or explicitly. Geometry-implicit dual-permeability (DP) models deal with the effects of pore-scale heterogeneity on preferential flow by treating the soil macropores and matrix as two separate interacting flow domains (Jarvis et al., 1991; Gerke and van Genuchten, 1993; Šimůnek et al., 2003; Larsbo et al., 2005; Šimůnek and van Genuchten, 2008). In some DP models, Darcy’s law is used to calculate water flow in both domains. One potential problem with this is that the underlying assumptions of Darcy’s law (e.g., that the flow is laminar and that the resistance to flow is dominated by the friction force between the water and solids) may not be appropriate for the flow of water in macropores (Jarvis, 2007). A parsimonious alternative is to simulate gravity-driven variably saturated flow in the macropore domain with a kinematic wave equation that can be derived from a generalization of pore-scale flow equations (Beven and Germann, 2013). Dual-permeability models implicitly assume that a representative elementary volume (REV) can be defined for macropores at the relevant scale of interest, for example the thickness of the unsaturated zone, which may not always hold true (Vogel and Roth, 2003). Nevertheless, DP models can often match measurements reasonably well (Köhne et al., 2009a, 2009b); they are also fast to run and numerically stable, and their structure is simple, which means that many chemical and biological processes can easily be incorporated into the models. For these reasons, DP models are now sometimes used as management and decision-making tools, for example in pesticide risk assessment (Jarvis et al., 2012; Jarvis and Larsbo, 2012). Two- and three-dimensional versions of DP models have also been developed that have been used to explore the impacts of pore-scale processes on soil hydrology and water quality in both drained and undrained soils at hillslope and field scales (e.g., Gärdenäs et al., 2006; Warsta et al., 2013; Gerke et al., 2013; Bishop et al., 2015; Dusek and Vogel, 2016). One disadvantage of DP models is that the properties of the macropore system are only implicitly reflected in model parameters that cannot easily be directly measured. This has stimulated research in recent years on methods to estimate these parameters by inverse modeling (e.g., Larsbo and Jarvis, 2005; Kodešová et al., 2010; Arora et al., 2012; Jarvis and Larsbo, 2012) and empirical pedotransfer functions (e.g., Moeys et al., 2012; Jarvis et al., 2012). Recent developments in noninvasive imaging techniques might also allow more direct estimation of some important model parameters. For example, Gerke (2012) suggested that macropore surface areas quantified by X-ray might help to constrain the values of the empirical mass transfer coefficient in a DP model.
The potential advantage of geometry-explicit models is that preferential flow at the Darcy scale is simulated as an emergent phenomenon dependent on the nature of the pore-scale features. An explicit representation of soil macropores strictly also requires pore-scale equations to simulate flow (e.g., the fundamental Navier–Stokes equations, the Lattice–Boltzmann approximations of them, or the Stokes equations that assume laminar flow). These equations have been applied directly to the pore space quantified by three-dimensional X-ray imaging techniques, both for porous rocks (e.g., Blunt et al., 2013; Bultreys et al., 2016) and occasionally also for soils (Hyväluoma et al., 2012; Khan et al., 2012). However, this approach has so far been limited to simulation domains that are much smaller than would constitute an REV for flow in macroporous soils because solving these pore-scale equations is computationally extremely demanding. One possible exception is the recent study of Scheibe et al. (2015). They used a hybrid multiscale approach on a supercomputer, combining X-ray imaging with pore-scale (Stokes) and Darcy flow equations to simulate water flow and tracer transport through a soil column. This kind of approach may prove to be a significant advance that should become increasingly viable at even larger scales in the future, as the necessary computer power becomes more widely available.
The practical limitations inherent in pore-scale computations have been circumvented by applying equations that are strictly only appropriate to collections of pores at the Darcy scale to describe flow and transport in individually defined macropores. In this way, two- and three-dimensional geometry-explicit modeling approaches have been applied to much larger scales (pedons and even hillslopes) than can be handled in practice by pore-scale flow equations (e.g., Nieber et al., 2006; Vogel et al., 2006; Nieber and Sidle, 2010; Klaus and Zehe, 2011; Wienhöfer and Zehe, 2014). Conceptually, this is a fundamental mismatch in scale, which can cause practical problems. As discussed by Nieber and Sidle (2010) and Alberti and Cey (2011), great care should be exercised in estimating the hydraulic properties of the macropores to avoid unrealistic parameterizations that would produce spurious simulation results. For example, water retention properties appropriate for coarse sand have been used to represent macropores that are supposedly several millimeters in diameter. It would not be surprising to find that such “virtual macropores” do not show the same hydrological behavior as real macropores of this size. Another challenge in parameterizing these models is how to define realistic macropore networks. One approach is to use “genetic” algorithms that aim to mimic the relevant structure-forming process. For example, Vogel et al. (2006) generated virtual earthworm burrows from the shortest path on the “backbone” of random percolation networks, while Klaus and Zehe (2011) generated stochastic distributions of earthworm channels using a correlated random walk algorithm parameterized from field observations. In numerical experiments performed to illustrate the effects of macropore continuity on preferential flow, individual macropores have been subjectively placed in the simulation domain (e.g., Allaire et al., 2002; Nieber et al., 2006; Akay et al., 2008; Nieber and Sidle, 2010; Alberti and Cey, 2011). Thus, in most applications of these models to date, rather simple systems have been defined that bear little resemblance to the complex multiscale networks found in natural soils. Furthermore, by necessity, the significant effects of micrometer-scale features such as macropore coatings and linings on flow and transport are also neglected. More realistic parameterizations of geometry-explicit models at the larger scales of interest for management must await substantial improvements in computational power.
Emerging Trends, Lessons Learned, and Future Prospects
The study of macropore flow at the pore scale has emerged as an important research theme in recent years, stimulated by advances in imaging technologies (and their increasing availability) and motivated by the desire to better understand the process at the most relevant scale (i.e., the scale at which it occurs). In this respect, we have not yet fully explored the possibilities offered by direct, noninvasive imaging of flow and transport experiments to improve our understanding of the flow mechanisms. It should be possible, for example, to directly test macropore flow models against pore-scale measurements of water contents and solute concentrations made on the same soil sample in flow and transport experiments under contrasting initial and boundary conditions. However, it is often problematic to repeat experiments on the same soil sample because it may cause temporal changes in pore space properties. In such cases, it may prove advantageous to make use of more durable three-dimensional printed copies of soil macropore networks imaged by X-ray (e.g., Bacher et al., 2015). Further improvements in process understanding can also be expected from studies that use two or more complementary noninvasive imaging techniques. For example, Badorreck et al. (2012) combined X-ray imaging, to visualize and quantify soil macroporosity as well as zones of compacted soil surrounding macrofaunal burrows, with neutron radiography, to directly image water infiltration at a high temporal resolution. Another possibility would be to combine X-ray with PET scanning to enable direct imaging of the transport of specific solutes at the small concentrations that usually occur in natural soils.
Direct imaging experiments might also help to resolve the question of the size of pores that enable nonequilibrium flow, something that remains a contentious issue even after 35 yr of research (Beven and Germann, 2013). Until now, the sizes of pores conducting preferential flow have been estimated indirectly using models to interpret flow and transport measurements. Hunt et al. (2013) noted that errors and misunderstandings can arise from such indirect approaches if they are based on oversimplified or incorrect conceptual models. This can be illustrated by studies that have made use of the Hagen–Poiseuille equation to calculate the sizes of conducting pores in soil from flow velocities inferred from TDR measurements of soil water content (Germann and Hensel, 2006; Beven and Germann, 2013). The impression could easily be gained from these studies that fast preferential flow through structured soils can occur in pores as small as 30 μm in diameter, but such an interpretation would be misleading. This is because the Hagen–Poiseuille equation overestimates the flow capacity of the complex pore networks found in soil (Bodhinayake and Si, 2004; Jarvis, 2008; Alberti and Cey, 2011), which is controlled by “bottleneck” resistances in the system (Hunt et al., 2013).
Another important unresolved question is whether current tomographic techniques (e.g., industrial and medical X-ray scanners) can combine the image resolution needed to quantify macropore networks (~300–500 μm) with a sufficiently large sample size that ensures an REV for the important network properties that control preferential flow (e.g., connectivity and continuity). To properly represent the topology of pore networks, the ratio between the field of view and the image resolution should be much larger than the ratio between the length and thickness of the pores (Vogel and Roth, 2003). This should not be an issue in topsoil (e.g., Jarvis et al., 2016), but these two ratios may be of a similar order of magnitude in many subsoil horizons of strongly anisotropic structure, where preferential flow may take place through earthworm burrows and shrinkage cracks that are continuous across large distances. Quite surprisingly, this question has hardly been investigated, even though it is critical for reliable parameter estimation and model prediction (though less so for testing and process understanding).
The connectivity of fast flow pathways has emerged as an important topic of research in recent years. The fact that preferential flow in soil is such a widespread phenomenon suggests that a continuous network of fast flow pathways exists in most soils, at least at the observation scale of these experiments, which is usually within the depth of the vadose zone that is influenced by structure-forming processes (i.e., the depth influenced by soil wetting–drying and biological activity, the plant rooting depth; Jarvis et al., 2012). Thus, what remains to be quantified is the size and properties of the connected or active fast-flow portion of the soil pore network. Smaller macropores are likely to play a critical role in this respect (Sidle et al., 2001) because they are likely to be activated more often than larger macropores under natural boundary conditions in the field and they are also more densely distributed and better connected. To date, the effects of connectivity have only been investigated in a superficial way in geometry-explicit models and have not been considered at all in dual-permeability models. It seems likely that some key concepts from percolation theory may prove useful to quantify the effects of connectivity in improved preferential flow models in an analogous way to models of hillslope runoff (e.g., Lehmann et al., 2007; Janzen and McDonnell, 2015). In addition to the network topology, the highly dynamic nature of soil macropore networks is also not well represented in the current generation of preferential flow models. Many different processes influence the dynamics of soil macroporosity and preferential flow, including soil freezing and thawing in cold climates (e.g., Frey et al., 2012; Hayashi, 2013; Lundberg et al., 2016), swell–shrink in clay soils (Warsta et al., 2013; Coppola et al., 2015), earthworm activity (van Schaik et al., 2014), and soil tillage and traffic (Schwen et al., 2011; Sandin et al., 2017). Furthermore, no existing simulation model accounts for all preferential flow mechanisms at both pore and Darcy scales, even though in many soils more than one may occur simultaneously (e.g., Jarvis et al., 2012; Hardie et al., 2013). We are also not aware of any model that considers the effects of water repellency on the generation of macropore flow, although many recent studies have demonstrated its importance, especially for uncultivated soils under natural vegetation (e.g., Jarvis et al., 2008; Nyman et al., 2010; Carrick et al., 2011; Badorreck et al., 2012; Hardie et al., 2012; Nimmo, 2012).
Thus, despite some important advances in recent years, for various reasons, current-generation simulation models still do not fully reflect the present state of empirical knowledge of preferential flow. However, we agree with the cautious optimism expressed a few years ago by Clothier et al. (2008). Even if some important questions still remain unresolved, our understanding of preferential flow processes in the soil and vadose zone is continuing to steadily improve, facilitated by the application of novel methods of measurement and data analysis both at the pore scale and at the larger scales that are relevant for management. Together with the significant advances that we expect to see in computational techniques, computer hardware, and measurement technologies, this improved process understanding should eventually lead to more reliable predictive modeling tools.