We would like to thank Peter Germann for his comments on our update paper, not least because it gives us an opportunity to further clarify some apparently contentious issues. In particular, Germann writes that we failed to discuss how fast water moves in preferential flow pathways, how much water is involved, and how far it is carried. This is certainly true. We had included such a discussion in an early draft of the paper, but this was deleted to meet the word limit for the format of an update paper. Luckily, we can now include a short discussion of these questions in our response to Germann’s comments.

Germann writes that we ignored Stokes’ equations. This is not strictly true as we did include some discussion on the potential use of the fundamental Navier–Stokes equations as well as the simplified Stokes approach that is limited to laminar flow. However, we seem to have a different perspective than Germann. He advocates the use of analytical solutions that are strictly appropriate only for simplified flow geometries and limited to steady-state applications. We focused our attention on numerical solutions of these equations at the pore scale because they can be applied to the transient and irregular flow geometries that are likely to occur in real (i.e., complex) soil pore networks under field conditions. Indeed, one of our main conclusions was that these approaches should become increasingly important in the future, especially in the context of “hybrid” model applications that combine continuum and pore-scale physics.

*q*(m s

^{−1}), can be expressed with a generalized Kozeny–Carman equation (Scheidegger, 1960) as where

*g*is the acceleration due to gravity (= 9.8 m s

^{−2}),

*G*is a dimensionless flow geometry factor that varies between 2 (cylindrical geometry) and 3 (planar geometry), μ is the kinematic viscosity (= 10

^{−6}m

^{2}s

^{−1}), τ (dimensionless) is the flow path tortuosity, θ is the water content (m

^{3}m

^{−3}), and

*R*

_{h}(m) is an effective “hydraulic radius” defined as θ divided by the wetted specific surface area,

*A*

_{w}(m

^{2}m

^{−3}). Analytical solutions of the Stokes equations that are applicable for simple idealized flow geometries (e.g., the Hagen–Poiseuille equation) can be seen as special cases of this more general model by setting the appropriate values for

*R*

_{h}in Eq. [1]). For example, by definition,

*R*

_{h}takes values of

*d*/4,

*d*/2, and

*d*for saturated cylindrical channels, saturated slits, and free planar films, respectively, where

*d*is the channel diameter, slit width, or film width.

*A*

_{w}(and therefore

*R*

_{h}) will vary with the water content, θ, as the flow rate,

*q*, changes. A general solution to Eq. [1] therefore requires knowledge of the two functional relations

*A*

_{w}(θ) and τ(θ). For the sake of brevity, we have assumed in the following that τ in Eq. [1]) is constant, which may be reasonable for flows through networks of large pores with a limited range of sizes (Jarvis, 2008). A kinematic wave model can be derived if

*A*

_{w}is expressed as a power law function of θ: where γ is a scaling factor and β reflects the spatial configuration of water flowing through the macropore network. Combining Eq. [1] and [2] and setting

*c*=

*g*/(

*G*μτγ

^{2}), we have where

*v*(=

*q/*θ) is the pore water velocity (m s

^{−1}). The pore water velocity can also be expressed in terms of the water flow rate by combining Eq. [3] and [4]:

*c*as the conductance parameter and 3 − 2β as the kinematic exponent. The difficulty of estimating some of the variables included in

*c*(i.e.,

*G*, τ, and γ) means that a modified form of Eq. [3] is more useful in practice (Beven and Germann, 1981). Normalizing the water content by the conducting porosity, ε, and ensuring that

*q*=

*K*

_{s}when θ = ε, where

*K*

_{s}is the saturated hydraulic conductivity of the preferential flow pore network, leads to where the kinematic exponent is now denoted by α (α = 3 − 2β) and

*K*

_{s}=

*c*ε

^{α}. For physically realistic simulations, the value of β should lie between 0 and 1 (i.e., α should vary between 1 and 3) to ensure that both the wetted area and pore water velocity are either constant or increase as the water content and flow rate increase. The kinematic wave equation is therefore not limited to the case when the kinematic exponent equals 3. The value of 3 suggested by Germann in his comment comes from only one of many possible flow geometries (i.e., a constant wetted surface area independent of θ, see Eq. [1]).

In the following, we use the KWE to try to shed some light on one of the important questions raised by Germann in his comment: how fast does water move in preferential flow? Apparently conflicting views have been put forward in recent years concerning this question, which we will here try to resolve. Based on an analysis of a small database of tracer tests performed under field conditions, Nimmo (2007) suggested that preferential flow may occur at a constant maximum pore water velocity in the vadose zone. This is equivalent to assuming that α = β = 1 (i.e., θ and *A*_{w} maintain proportionality during flow), so that the pore water velocity, *v*, is constant at the maximum value (= *c*) regardless of the flow rate and water content (Eq. [4] and [5]). One mechanism that may explain this is the flow of “rivulets” in large macropores (e.g., Bouma and Dekker, 1978; Germann et al., 2007; Cey and Rudolph, 2009). This would give a constant *v* irrespective of *q* if the number of rivulets increased as *q* increased but their thickness remained nearly constant. Hincapié and Germann (2009) found that *v* was proportional to *q*^{2/3} for a soil column irrigated at varying rates in the laboratory, which implies that wetted surface areas changed very little with flow rate (i.e., β ≈ 0, see Eq. [2] and [5]). From this result and those of other studies, Beven and Germann (2013) suggested that β ≈ 0 (and α ≈ 3) might be a useful general working hypothesis. This implies the opposite viewpoint to the concept proposed by Nimmo (2007), since for a given flow rate *q*, *v* is smallest when β = 0. However, neither concept can be considered as a generally valid model of preferential flow, although β = 1 should be a reasonable basis for a simple worst-case characterization, as suggested by Nimmo (2007). Instead, we should expect the velocity of water moving through preferential flow paths in the soil and vadose zone to vary between these two extremes as a function of boundary conditions and the macropore network properties that influence γ, τ, and β. This is also what the experimental data suggest: in the small data set collated by Nimmo (2007), individual values of maximum *v* varied over approximately two orders of magnitude for experiments conducted with continuous water and tracer application. In the data set presented by Germann and Hensel (2006), the measurements of *v* varied over at least one order of magnitude.

In our update paper, we wrote that the sizes of pores supporting preferential flow in structured soil still seems to be a contentious issue some 35 yr after it was first debated. We further suggested that some unnecessary confusion surrounding this question has arisen in recent years from the inappropriate application of the Hagen–Poiseuille equation to profile-scale measurements of wetting front velocities and “mobile” water contents. Germann suggests that we did not state our reasons for objecting to this approach. We agree that our explanation was perhaps rather too concise and therefore incomplete. The Hagen–Poiseuille equation can be applied to reliably calculate flow rates through interconnected networks of saturated tubes, as is done in pore network models (e.g., Köhne et al., 2011), but it is not an appropriate method to estimate conducting pore sizes from field-measured wetting front velocities. This is because it wrongly assumes that the water flows from the soil surface to the measurement depth in vertical, non-intersecting, equal-sized tubes. As we wrote, this will “…overestimate the flow capacity of the complex pore networks found in soil, which is controlled by bottleneck resistances in the system.” This approach also seems to encourage confusion between the size of pores conducting preferential flow and the thicknesses of the water flow pathways within them. This can be illustrated by citing from Germann and Hensel (2006). Introducing their approach, they first wrote: “…that portraying the geometry of flow with idealized structures is inherently different from analyzing the pores that are occupied by the moving water; the latter is not the subject of this manuscript.” However, this vital point was unfortunately forgotten when they summarized their results. In the final paragraph of the conclusions, they wrote “…finally, we note that the radii of most of the equivalent Poiseuille pores do not indicate typical macropores, but more likely belong to the class of mesopores.” Furthermore, in the abstract, they wrote, “Application of Poiseuille’s Law to the velocities resulted in radii and densities of equivalent Poiseuille pores in the ranges of 5 to 30 μm and of 2 × 10^{6} to 2 × 10^{8} m^{−2}, respectively.” For the reasons outlined above, the conclusion that preferential flow occurred in pores of this size cannot be supported by the kinds of measurements they performed.

In a later review, Beven and Germann (2013) suggested that film flow widths of 15 to 30 μm calculated with the same method do not “…correspond with the typical perception of preferential flow along well-visible macropores*.*” In our view, the unsaturated flow of thin water films and rivulets in large (i.e., visible) macropores is entirely compatible with the typical perception of preferential water flow in structured soils. It is interesting that Germann now refers in his comment to the results of these calculations as providing estimates of the “minimal width of preferred-flow paths” rather than to pore sizes. This is probably a more reasonable way of expressing what the calculations might show, although the derived minimum widths could be considerable underestimates of the actual (average) width due to the effects of bottleneck constrictions. Unfortunately, having first referred to minimal flow path widths, Germann adds to the confusion later in his comment by writing that “…the procedure yields densities and diameters of cylindrical *pores* in the ranges from 10^{6} to 10^{8} m^{−2} and from 7 to 30 μm” [our emphasis]. Finally, Germann writes that we failed to refer to any sources when we “…postulated path widths of preferential flow of approximately 300– 500 μm.” First, we did not postulate *path widths* of this size in our paper. Instead, we suggested that *pores* of roughly this thickness (and larger) generally support preferential flow in structured soil (while noting also that the long-range continuity of such pores is also important). However, we agree that it was remiss of us not to cite studies to support this statement. To correct this omission, we would refer to Jarvis (2007), who discussed the various sources of experimental evidence. Most obviously, sprinkler experiments using dye to visualize the flow pathways in structured soils have clearly highlighted the dominant role of large and sparsely distributed macropores in conducting preferential flow. They also do not support the idea that there are typically between 1 million and 100 million microscopically small preferential flow pathways per square meter.

Germann also criticized our treatment of the question of representative elementary volumes (REVs) for preferential flow as being inconclusive. In one sense this is probably true, but this mostly reflects the fact that little work has been done on this topic. Germann writes that “…a REV is required in capillary flow mainly because of the relationship between total water content and capillary potential.” We certainly agree that Richards’ equation as well as geometry-implicit preferential flow models like the KWE require a REV for the soil hydraulic properties (this was noted in our paper). Germann writes that REVs are irrelevant for the analytical models of the Stokes type that he advocates. This is true, but in our view this does not compensate for their limitations. As mentioned above, we think that numerical solutions of these equations for explicitly defined pore networks embedded within a hybrid modeling framework (e.g., Scheibe et al., 2015) will prove more useful in a long-term perspective. However, the mismatch between typical measurement scales and the usual scales of interest for model applications (e.g., fields, hillslopes) means that the scale dependence of water flow rates associated with soil heterogeneities is a relevant and important question even for these pore-scale models (Scheibe et al., 2015; Pachepsky and Hill, 2017). This issue can only be avoided by assuming that the heterogeneity of pore system properties exerts no control on water flow, which seems unreasonable to us.