**Analysis of the Interconnections between Classic Vortex Models of Coherent Structures Based on DNS Data**

#### **Hao Wang \*, Guoping Peng, Ming Chen and Jieling Fan**

College of Civil Engineering, Fuzhou University, Fuzhou 350116, China; pengguoping\_1314@163.com (G.P.); korry\_chen@163.com (M.C.); fanjieling318@163.com (J.F.)

**\*** Correspondence: wanghao\_0512@163.com

Received: 7 August 2019; Accepted: 23 September 2019; Published: 26 September 2019

**Abstract:** Low- and high-speed streaks (ejection, Q2, and sweep, Q4, events in quadrant analysis) are significant features of coherent structures in turbulent flow. Streak formation is closely related to turbulent structures in several vortex models, such as attached eddy models, streamwise vortex analysis models, and hairpin vortex models, which are all standard models. Vortex models are complex, whereby the relationships among the different vortex models are unclear; thus, further studies are still needed to complete our understanding of vortices. In this study, 30 sets of direct numerical simulation (DNS) data were obtained to analyze the mechanics of the formation of coherent structures. Image processing techniques and statistical analysis were used to identify and quantify streak characteristics. We used a method of vortex recognition to extract spanwise vortices in the *x–z* plane. Analysis of the interactions among coherent structures showed that the three standard vortex models all gave reasonably close results. The attached eddy vortex model provides a good explanation of the linear dimensions of streaky structures with respect to the water depth and Q2 and Q4 events, whereby it can be augmented to form the quasi-streamwise vortex model. The legs of a hairpin vortex envelop low-speed streaky structures and so move in the streamwise direction; lower-velocity vortex legs also gradually accumulate into a streamwise vortex. Statistical analysis allowed us to combine our present results with some previous research results to propose a mechanism for the formation of streaky structures. This study provides a deeper understanding of the interrelationships among different vortex models.

**Keywords:** image processing; streaky structures; hairpin vortex; attached-eddy vortex; streamwise vortex

#### **1. Introduction**

Turbulence is generally not altogether chaotic, whereas there are many regular coherent structures in a turbulent flow. The coherent structures include streaky structures formed by low- and high-speed streaks, the bursting phenomenon that includes ejection and sweep events (in quadrants Q2 and Q4), vortex structure models (streamwise vortex model, attached eddy vortex model, hairpin vortex model and hairpin vortex groups), as well as superscale structures.

Low- and high-speed streaks are important in turbulence dynamics because of their large scale [1]. Experimental research into low- and high-speed streaks using hydrogen bubbles was first conducted by Kline et al. [2]. The characteristic scales of streaky structures were also identified by many researchers as the average nondimensional width *W* = 20–40*y*\* and spanwise distance *D* = 100*y*\* in the boundary layer region [3–5]. Note that *y*\* = *v*/*u*\* defines the inner scale, where *v* is kinematic viscosity and *u*\* is friction velocity, which represents the shear stress velocity. Lin et al. [6] used particle image velocimetry

(PIV) to capture the flow fields. Their results show that the spatial distribution of high-speed streaks is similar to that of low-speed streaks.

Zhong et al. [1] identified elongated streamwise low- and high-speed streaks near the free surface in open-channel flows by spanwise correlation analysis. The presence of large-scale streaks across the whole flow depth has been confirmed by many researchers. Previous evidence indicates that the distance between neighboring low-speed streaks is the water depth scale (*H*–2*H*) [7–10]. Sukhodolov et al. [11] found that streamwise streak length could exceed 3*H* while Zhong et al. [1] found the length to be greater than 10*H*. The existence of streaky structures throughout the whole turbulent layer is now commonly accepted [7,12–14].

Various hypotheses and models of vortices have been created to explain the formation of low- and high-speed streaky structures. Many studies proposed super-streamwise vortex models of Q2/Q4 events, which included alternating low- and high-speed streaks in the spanwise direction [9,15,16]. The attached eddy hypothesis developed by Townsend [17] explained Q2/Q4 events and the development of streaky structures, which scale linearly with their water depth from the inner region to the outer region. Adrian and Marusic [18] advocated a model using hairpin vortices and packets: hairpins and packets cause the ejection of low-speed streaks between the two legs of the hairpin vortex when the super-streamwise vortices feed themselves by sweeping low-momentum hairpins and packets into the low-speed regions. Secondary flow cells have also been modelled as vortices which originate in the vicinity of the side walls [19,20].

The existing research indicates that vortex models have limited use. Researchers accept the existence of super-streamwise vortices theoretically, but the literature reviewed above shows that there is no consensus among researchers concerning the formation of vortices. For example, the streamwise vortex model cannot explain how streak length varies linearly from inner region to outer region. The attached eddy vortex is a single structure, which does not explain the distribution and organization of the many funnel vortices in turbulent flow. Hairpin vortex models are usually developed for a single flow field and vortex structure in the *x*–*y* or *x–z* plane. However, current understanding of the characteristics of hairpin vortices is insufficient to generate a robust interpretational theory. There are relatively few studies of vortex models, and thus there is a lack of systematic quantitative vortex model analysis; vortex models can still be improved.

We used models to investigate vortices as coherent structures in turbulent flow, using direct numerical simulation (DNS) data. We identified the positions of low- and high-speed streaks using image processing and calculated the characteristic dimensions of streaky structures in both the inner and outer layers using a statistical method. We identified streamwise vortices, attached eddy vortices, and hairpin vortices by analyzing the variation in streak dimensions with respect to water depth and analyzed the spatial relationships between streaky structures and spanwise vortex position to explain the relationship between the three vortex models. Finally, we propose a new hypothesis.

The remainder of this paper is organized as follows. Section 2 describes the methods used to analyze the DNS data and to identify and calculate the characteristic dimensions of streaky structures. Section 3 presents an analysis of the regular variation in streaky structures and the mechanics of the three vortex models. Section 4 offers a summary and a brief discussion of our major findings and the conclusions we draw from them.

#### **2. Materials and Methods**

#### *2.1. Closed Channel Flow: DNS*

Particle image velocimetry (PIV) is the principal experimental method of measuring the flow field. The area captured by the camera is relatively small, due to the limited intensity of laser light, as the physical width (*z* direction) of the image. Thus the number of low- and high-speed streaks sampled is relatively small, whereby the characteristic scale of streaks is not particularly accurate. Therefore, we used the numerical data of Del Alamo et al. [21] to investigate coherent structures in turbulent flow,

and thereby obtained complete flow field information for closed channel flow. Data series L950, which we used extensively, contains data for almost all recognized large-scale coherent structures scaled by water depth [5,21,22]. Figure 1 shows that the dimensionless length and width of the DNS flow field are both large, whereby the characteristic scales of the streaks are relatively accurate. We give a brief introductory summary here. Detailed information can be found in Del Alamo et al. [21].

**Figure 1.** The computation region of direct numerical simulation (DNS) in closed channel flow.

The friction Reynolds number for the flow was 934, which indicates that the range of temporal and spatial fluid scales involved in turbulence was considered to be relatively large. The simulation covers a spatial domain (*x*, *y*, *z*) of 16π*h*/3 × 1*h* × 2π*h*, where *h* is the half-channel height and the domain is discretized into an array (*x* × *y* × *z*) of 2048 × 193 × 1536 points. Each grid point contains three velocity components corresponding to nine velocity gradient data points. The streamwise (*x*), vertical (*y*), and spanwise (*z*) dimensions are shown in Figure 1, which summarizes of the DNS data; *u*, *v*, and *w* represent the instantaneous velocities in the *x*, *y*, and *z* directions, respectively. Major parameters of the DNS data are summarized in Table 1.

**Table 1.** Parameters of the DNS (data from Del Alamo et al. [21]).


In Table 1, *H* is the water depth; *Lx*, *Ly*, and *Lz* are the spatial domains along the *x*, *y* and *z* directions, respectively; Δ*x* and Δ*z* are the grid resolutions in the *x* and *z* directions, respectively; *Nx* and *Nz* correspond to the grid numbers; Δ*yc* is wall-normal grid spacing at the channel center; *Ny* represents the grid numbers along the *y* direction; the superscript <sup>+</sup> denotes normalization by the inner scale (*u*\* and *v*); and *u*\* is the friction velocity and represents the shear stress velocity, for example, Δ*x*<sup>+</sup> = Δ*xu*\*/*v.*

For each of the three-dimensional instantaneous velocity fields, totals of 153 × 30 *x–y* planes, 204 × 30 *y–z* planes, and 193 × 30 *x–z* planes were extracted for analysis. There were 2048 × 124 (*x–y* plane), 245 × 1536 (*y–z* plane), and 2048 × 1536 (*x–z* plane) grid points. We extracted 193 × 30 *x–z* planes for analysis, and there are 2048 × 1536 grid points in each *x–z* plane.

#### *2.2. Detection of Streaky Structures*

The formation of low- and high-speed streaks are related to instantaneous turbulence fluctuations. Three steps were followed to study the characteristics scale of streaky structures: (1) the detection function was used to identify the high- and low-speed streaks; (2) image processing, including binarization and morphological operations, was used to extract the image structure of both low- and high-speed streaks [6,10]; and (3) statistical analysis was used to calculate the characteristic scales of streaks.

#### 2.2.1. Detection Function

The method, after modification, used the following two functions,

$$F\_d(m, n, y^+, t) = \frac{u'(m, n, y^+, t)}{u\_{std}(y^+)}\tag{1}$$

$$\mathbb{C}t(\boldsymbol{y}^{+}) = \mathbb{C} \times \max[\boldsymbol{u}\_{\text{std}}] / \boldsymbol{u}\_{\text{std}}(\boldsymbol{y}^{+}) \tag{2}$$

where (*m*, *n*) is the grid position in the *x–z* plane; *u'* is the streamwise velocity fluctuation; *ustd*(*y*+) is the standard deviation of the streamwise velocity at *y*+; *Ct*(*y*+) is the water depth threshold at *y*+; *Fd* is the dimensionless value of detection function; *C* is a constant, equal to 0.6, as recommended by Lin et al. [6]; and max[*ustd*] is the maximum value of *ustd* in the flow domain. *Fd* > *Ct* (high-speed) and *Fd* < −*Ct* (low-speed) identify the streaks. Justification for the two equations and specific details are provided in Wang et al. [10].

Figure 2a shows the contours of *Fd* for low-speed streaks at *y*<sup>+</sup> = 21.05. The positive and negative values of *Fd* indicate the existence of instantaneous streamline fluctuations, forming the low- and high-speed streaks. Low-speed regions (brown), high-speed regions (blue), and other flow regions (green) can be recognized by applying a threshold value of *Ct*(*y*+) to the contour map, as shown in Figure 2b.

(**b**) After applying threshold to *Fd.*

**Figure 2.** Visualization of streaks represented by the dimensionless value of detection function *Fd*: (**a**) original *Fd* with the range of the color bar set from −2.5 to 2.5; (**b**) after applying the threshold value to *Fd*.

#### 2.2.2. Image Processing

To better quantitatively analyze the low- and high-speed streaks, a binary procedure was used to extract the streaks: values less than −*Ct* were assigned the value 1, whereas values greater than –*Ct* were assigned a value 0. Figure 3 shows the image processing procedure for extracting low-speed streaks. The procedure for extracting high-speed streaks is similar, but uses a different *Ct* threshold value.

**Figure 3.** *Cont.*

**Figure 3.** Image processing: (**a**) binary image, (**b**) opening operator, (**c**) closing operator, and (**d**) clean image.

The original *Fd* image was binarized (Figure 3a), and a basic morphological transformation was used to filter out noise in the binary images. This transformation was done in two steps [23]. First, the opening operator (Figure 3b) and closing operator (Figure 3c) were used to delete some isolated regions and fill some holes. The opening operator is derived from the fundamental morphological operations of dilation as well as erosion and was used to break the adhesion between objects and remove small particle noise; the closing operator combines the operations of erosion and dilation and can be used to connect neighboring regions and fill in small holes. The area of the streak graph does not change significantly during calculation when using the opening and closing operators.

Second, some isolated objects were deleted from the binary image with the *bwareaopen* function. After these two steps were performed, the streaky structures were clearly visible (Figure 3d). The selection of specific parameters and values is described by Lin et al. and Wang et al. [6,10].

#### 2.2.3. Model of Streaky Structures

We used streak width (*w*) and the distance between adjacent streaks (*d*) to scale and characterize streaky structures. As the structures vary spatiotemporally, the image was parsed line-by-line to quantify both *w* and *d*. We assumed that the number of streaky structures in a line of the image was *ns*. The streak widths are denoted by *w*1, ... , *wi*−1, *wi*, *wi*+1, ... , *wns* when the streak distances are denoted by *d*1, ... , *di*−1, *di*, *di*+1, ... , *dns*−1.

For the random row in the image (*r*th line), the mean streak width of the *r*th line, *w*(*r*), and the mean streak width of the whole velocity field at each *y*+, *W* (*y*<sup>+</sup>), were obtained by Equations (3) and (4):

$$w\_{\binom{r}{r}} = \frac{\sum\_{i=1}^{ns} w\_i}{ns}, (r = 1, 2, \dots, m) \qquad \mathcal{W}\_{y^+} = \frac{\sum\_{r=1}^{m} w\_{\binom{r}{r}}}{m} \tag{3}$$

where *m* is the total number of rows of the flow field. Similarly, the mean spanwise distance of each row, *d*(*r*), and the mean spanwise distance of the whole velocity field at each *y*+, *D*(*<sup>y</sup>* <sup>+</sup>), were calculated by

$$d\_{\boldsymbol{\alpha}\_{(r)}} = \frac{\sum\_{i=1}^{ns-1} (d\_{i+1} - d\_i)}{ns - 1}, (r = 1, 2, \dots, m) \qquad D\_{(y^+)} = \frac{\sum\_{r=1}^m d\_{\boldsymbol{\alpha}\_{(r)}}}{m} \tag{4}$$

The dimensionless characteristic scales of the streaks can be obtained by *D*<sup>+</sup> = *Du*\*/*v* and *W*<sup>+</sup> = *Wu*\*/*v*. The 30 instantaneous *x–z* velocity fields (DNS data) were captured at each *y* position in all cases.

Comparisons between characteristic scales and previous data are given in Figure 8 of Wang et al. [10]. The variations in the mean spanwise distance relative to calculated wall distance are completely feasible, and the above method can be used to analyze the characteristic scales of streaky structures.

#### **3. Results**

The relationship between the spanwise distance between streaks and the vortex model used is significant for analysis of the entire phenomenon. The spanwise interstreak distances for both low- and high-speed streaks *D*/*H* and water depth *y*/*H* were plotted. Figure 4 shows that the trend of high-speed streaks is similar to that of low-speed streaks. As water depth *y*/*H* increases, *D*/*H* reaches a turning point close to half the water depth of the closed-channel flow. The increased amplitude decreases significantly near the half water depth due to a weak boundary layer.

**Figure 4.** Spanwise distances of streaks; *D*/*H* varies with water depth *y*/*H*.

Most research on streaky structures has been concerned with the inner region (*y*+) and outer region (*H*) scales. The relationship between two streaky structures of different scales is unclear. We calculated the spanwise distances over the entire water depth continuously for both the inner and outer regions. The results show that the development of streaky structures along the entire flow depth is a continuous process.

*D*/*H* increases linearly with *y*/*H* in the outer layer (i.e., when 0.1 < *y*/*H* < 0.8), and the slope is approximately 2. Streaky structures are closely linked to the vortex model used. Our results are interpreted in the context of the streamwise and attached eddy vortex models as follows.

#### *3.1. Streamwise Vortex Model*

The spanwise distances are approximately twice the water depth, and the formation of streaky structures in the outer layer is related only to water depth. This is consistent with streamwise vortex structures being generated automatically from the self-organization of wall-bound turbulence [7]. In this case, the streamwise vortex also shows that the strong pumping action of low-speed fluid creates an ejection event in the associated second quadrant, Q2 (*u* < 0, *v* > 0). Fluid moving at high speed from the water surface toward the bed creates a sweep event in Q4 (*v* < 0, *u* > 0). The low- and high-speed streaks are located at the downwelling and upwelling sides of the streamwise vortices respectively.

The conceptual streamwise model in Figure 5 was built to explain the formation of streaky structures and Q2/Q4 events; the distance between adjacent streaks near the water surface is twice the water depth [1,14,18].

**Figure 5.** Representation of the streamwise vortex model.

#### *3.2. Attached Eddy Vortex Model*

Figure 6a shows that the increases in the characteristic scales of the streaks are linear, which is a core assumption of Townsend's attached eddy vortex hypothesis [17,23]. Our results indicate the suitability of the attached eddy vortex model.

Figure 6a shows that, according to the model, the attached eddy vortex develops from the near-wall region into a conical vortex in the streamflow direction. The formation mechanics of low- and high-speed streaks and the Q2/Q4 events are similar to those of the streamwise vortex. In particular, the spanwise distances between adjacent low-speed (or high-speed) streaks are closely related to the size of the streamwise vortices and thus linearly proportional to *y*. Three cross-stream plane sections (*slices*) of the attached eddies at different water depths *y* are shown in Figure 6; the blue and red backgrounds indicate the low-speed and high-speed streaks, respectively. Figure 6a shows the Q2\Q4 events when Figure 6b shows that the dimensions of the vortex increase linearly as water depth increases and the relationship between *D* and *h* remains basically constant with *D* ≈ 2*h*. Overall, the attached eddy vortex model explains the distance between adjacent streaks from the inner region to the outer region.

(**a**) Diagram of the streaky structures and a pair of attached-eddy hypothesis.

#### **Figure 6.** *Cont.*

(**b**) High- and low-speed streaky structures for each slice.

**Figure 6.** Streaky structures described in terms of the attached eddy hypothesis: (**a**) A pair of vortices and the (**b**) corresponding high- and low-speed streaky structures for each slice.

#### *3.3. Hairpin Vortex Model*

The hairpin vortex model developed by Adrian is a relatively new vortex model [24]. We investigated the positional relationships between spanwise vortices and streaky structures in the *x–z* plane to determine the suitability of the hairpin vortex model.

#### 3.3.1. Vortex Extraction in the X–Z Plane

We used the two-dimensional swirling-strength λ*ci*-criterion [25]. Streaky structures form in the in *x–z* plane, so a brief introduction to extracting a vortex in the *x–z* plane is given.

The swirling strength λ*ci* is given by

$$
\lambda\_{ci} = \begin{cases}
\sqrt{R - P^2/4}, & 4R - P^2 > 0 \; ; \\
0 & 4R - P^2 \le 0
\end{cases}
\tag{5}
$$

where

$$P = -\frac{\partial \mathbf{u}}{\partial \mathbf{x}} - \frac{\partial w}{\partial z}, R = \frac{\partial \mathbf{u}}{\partial \mathbf{x}} \frac{\partial w}{\partial z} - \frac{\partial \mathbf{u}}{\partial y} \frac{\partial w}{\partial z} \tag{6}$$

Following Wu and Christensen [26], we defined the normalized swirling strength Λ*ci* as Λ*ci* = λ*ci*ω*z*/|ω*z*|, where ω*<sup>z</sup>* is the fluctuating spanwise vorticity and λ*ci* and Λ*ci* are swirling strength discriminators. Λ*rms ci* (*y*) is the local root mean square of Λ*ci* at the wall-normal position *y*, and we defined the normalized swirling strength Ω*ci* by

$$
\Omega\_{ci}(x, y) = \frac{\Lambda\_{ci}(x, y)}{\Lambda\_{ci}^{rms}(y)} \tag{7}
$$

In an ideal fluid, there is a clear boundary between rotating and irrotational fluid. Zero (0) can be used as a threshold to easily extract the vortex. However, in a nonideal (actual) fluid, viscosity causes dissipation of the vortex, which greatly complicates vortex identification. We used a non-zero threshold of 1.5 to identify a vortex, following the recommendation of Wu and Christensen [26], so that

$$|\Omega\_{\rm ci}| \ge 1.5\tag{8}$$

Negative or positive values of Ω*ci* in Equation (8) represent a clockwise or counterclockwise vortex, respectively.

#### 3.3.2. Vortex Density

When the vortex structure has been determined by the preceding methods, the vortex population density Π<sup>+</sup> is calculated by

$$\Pi^{+} = \frac{N\_{\text{vortex}}(y^{+})}{Nx(y^{+}) \cdot Nz(y^{+}) \cdot \Delta x^{+} \cdot \Delta z^{+}} \tag{9}$$

where *Nvortex* is the spanwise number of vortices at position *y*+.

The prograde and retrograde vortices in the *x–z* plane were separated, and the population densities of the vortices were computed for each position of *y*+. Figure 7 shows that the population density of 2D vortices varies with water depth, reaching a maximum in the near-wall region at *y*<sup>+</sup> = 40.22. This result may be partly due to the number of streamwise vortices in the deeper water. In the outer region, vortex density decreases gradually as *y*<sup>+</sup> increases. As the shear stress in the *x–z* plane is approximately zero, the population density of prograde vortices is equal to that of retrograde vortices at each *y*<sup>+</sup> position. These results agree with those obtained by Chen et al. [27], which confirms the logic of our vortex extraction method.

**Figure 7.** Vortex population density in the *x–z* plane.

#### 3.3.3. Location of Vortices and Streaks

Figure 8 shows the cores of spanwise vortices surrounded by nine velocity vectors (red), high-speed streaks (yellow), low-speed streaks (blue), and the in-between region (green) that were obtained by the preceding methods. We extracted and analyzed 400×800 grid points in the *x–z* plane. The dimensionless area is 3040 <sup>×</sup> 3040, and *x*<sup>+</sup> and *z*<sup>+</sup> represent the dimensionless length along the streamwise and spanwise directions, respectively.

**Figure 8.** Positional distributions of streaks and spanwise vortices.

To further investigate the relationship between streaky structures and spanwise vortices, the vortex core (*ui*,*j*) velocity and the eight surrounding streamwise velocities (Figure 9) were averaged using Equation (10). Equation (11) was then used to calculate *Fd* at the core of each vortex. We use *Ct*(*y*) to identify the region of the distribution of vortex cores.

$$
\overline{u} = \frac{1}{9} \sum\_{i=1}^{9} u\_i \tag{10}
$$

$$\overline{F\_d} = \frac{\overline{u}'}{\mu\_{std}(y^+)} \tag{11}$$

where *u¯* is the average velocity of spanwise vortices in the *x–z* plane, *u¯* is the fluctuating velocity, and *Fd* is the detection function of average vortex velocity. We can now calculate the statistical measures of the vortices in different streaky structures.

**Figure 9.** Grid for streamwise velocity *u* showing the core of the vortex.

In previous research, streaks have generally been divided into low- and high-speed. However, to obtain a more detailed analysis of the relationships between spanwise vortices and streaky structures, we divided the *x–z* plane into three types using the threshold *Ct*(*y*+): low-speed streaks, high-speed streaks, and in-between regions.

Figure 10 shows that the numbers of vortices differ greatly within the different streaks. The spanwise vortices in the in-between region occur in the greatest numbers, followed by low-speed streaks, with the least numbers in the high-speed streaks. When *y*/*H* > 0.1 the number of vortices in high-speed streaks is approximately equal to 0.

**Figure 10.** Vortex numbers in low- and high-speed streaks and in-between regions.

It should be noted that the three areas occupied by the three types are also different. By dividing the number of vortices by the area containing them, we eliminated the influence of area from our analysis to better understand the distributions of spanwise vortices located in streaky structures.

#### 3.3.4. Vortex Density in Different Streaks

We first calculated the areas of the low- and high-speed streaks. The values of both low- and high-speed streak widths are influenced by the threshold value (*Ct*) and are difficult to recognize relative to the spanwise distance. If *Ct* is too large, the streak width will not include the whole with of the streak; nevertheless, if *Ct* is too small, the streak width will contain parts of the in-between region. However, even with different threshold values, the change tendencies are basically consistent. Here, we used the threshold value suggested by Lin et al. [6] to identify the width of low- and high-speed streaks. Figure 11 shows that as *y*<sup>+</sup> increases, streak width first increases and then decreases. As the water depth increases from inner region to outer region, the streak scale also increases. When water depth is close to the surface (about 0.7*H*), the weak boundary layer restrains the streak scale, and the streak width will decrease. This trend is stable and clearly demonstrated.

**Figure 11.** Characteristic dimensions (width) of low- and high-speed streaks varying with *y*+.

The area percentages of low- and high-speed streaks in the *x–z* plane at each position *y*+ are also important characteristic dimensions and can be regarded as the normalized areas of streaks. The area percentage, defined as *Ps*, can be obtained by

$$A\_{\delta} = \sum\_{i=1}^{m} n s\_i \cdot w\_i \tag{12}$$

$$P\_s = \frac{A\_s}{A\_t} \times 100\% \tag{13}$$

where *As* is the total area of low- or high-speed streaks and *At* is the total area of the flow field (2048 × 1536). Figure 12 shows that the percentage area of both low- and high-speed streaks decreases as *y*<sup>+</sup> increases. In the near-wall region (*y*/*h* < 0.1), the gradients are steep; in the outer region (*y*/*h* > 0.3), the decreasing trends become less steep. This result shows that the streaks occur mainly in the near-wall region where there is high shear stress. Mean shear stress decreases as *y*/*h* increases, and so does its effect on the streaks. Streaks in the outer region decrease in number and so the percentages of both low- and high-speed streaks in the *x–z* plane also decrease.

**Figure 12.** Percentage area of low- and high-speed streaks along the *x–z* plane.

Figures 11 and 12 both show that the formation of streaks between the near-wall region and the outer layer is a continuous process, as also shown in Figure 4.

#### 3.3.5. Calculation of Vortex Density

Equation (10) was used to calculate the density of spanwise vortices in different streaky structures, as shown in Figure 13. The density of spanwise vortices is highest in low-speed streaks, intermediate in the in-between region and least in the high-speed streaks. Thus, there are big differences between the number of vortices and vortex density located in different streaky structures.

**Figure 13.** Vortex densities located in different streaky areas.

The hairpin vortex model developed by Adrian [24] has gained widespread acceptance, due to experimental visualization using particle image velocimetry and direct numerical simulation. Figure 14 shows the standard coherent structure model of a hairpin vortex developed by Adrian [24]. The alignment of coherent vortices induces a low-speed fluid region inside the hairpin packets. Due to the closed-loop feedback cycle between hairpin vortex cells and streamwise vortices [1,28], the streamwise vortices are stable.

**Figure 14.** Representation of hairpins and hairpin packets by Adrian [24].

Research into hairpin vortex behavior has become an important direction of research. However, the sample sizes used in the research are fairly small, so the regularity of the relationship between streaky structures and spanwise vortices in the *x–z* plane must be further researched by analyzing large samples.

We obtained statistics from large DNS data samples, and we found that spanwise vortex density in low-speed streaks is greater than in high-speed streaks. This result indicates that hairpin vortex legs are closer to the low-speed streaks and further from the high-speed streaks. Thus, the results we obtained exhibit an important feature of hairpin vortex legs when they envelop low-speed streaks to move along the quasi-streamwise direction, as shown in Figure 14. The legs of the hairpin vortex are spanwise vortices in the *x–z* plane, as shown in Figure 15. Spanwise vortices are mainly distributed in the region of low-speed streaks, consistent with the structure of hairpin vortices. Our results support the logic of the hairpin vortex model and reveal mechanisms of hairpin vortex behavior more explicitly.

**Figure 15.** Representation of hairpin vortex and spanwise vortex.

#### **4. Discussion and Conclusions**

#### *4.1. Discussion*

A large-scale vortex is a conceptual model, or representation, of a natural phenomenon intended to be used in the provision of logical explanations of all kinds of coherent structures in turbulent fluids. Classical and prevailing views of vortices have led to many vortex models being developed. We identified both low- and high-speed streaks from the wall to the surface using image processing technology; the meandering large scale motions are impossible to ignore. The low- and high-speed streaks are formed by an ejection event (Q2, *u* < 0, *v* > 0) and a sweep event (Q4, *v* < 0, *u* > 0) [14,18]. We used the super-streamwise vortex (Figure 5) as the interpretative model to explain the preceding results and the spanwise distance of nearby streaks (2*H*). We found that the scale of the streaks increased in proportion to their distance from the wall. The result is consistent with the classical model, which combines length growth with growth in eddies, developed by Townsend [17]. Our results also explain the logarithmic growth in open channel flow. The distributions of spanwise vortex density in low- and high-speed streaky structures suggest further research into hairpin vortices. We statistically sampled large datasets to compare and analyze three vortex models. Our analysis of the results shows the benefit of explaining coherent structures from the three different model perspectives.

The literature contains little record of large dataset statistical sampling, but it is urgently needed to demonstrate the suitability of different vortex models and to clarify the relationships between them. As stated in the introduction, the large-scale streamwise vortex model provides a good explanation of the coherent structures of Q2/Q4 events and the spanwise distances between adjacent streaky structures near the water surface (which is about 2*H*). However, the large-scale streamwise vortex model is relatively coarse and represents a large structure (Figure 5), and it cannot accurately explain the continuous development of streaks from the inner region to the outer region. The attached eddy vortex model cannot provide a precise organized structure for the large vortices that accumulate in turbulence. The hairpin vortex model requires more usage and analysis to show its suitability.

Vortex models are limited. However, our research into the characteristic dimensions of streaky structures across the entire water depth, described in this study, leads us to conclude: the streamwise vortex model, the attached eddy vortex model, and the hairpin vortex model are all suitable models in certain circumstances.

We used quantitative analysis to develop a theoretical model in which packets of attached eddy vortices self-organize and accumulate along the flow direction, thereby forming a cumulative vortex structure, the streamwise vortex. Figure 16 shows that many attached eddy vortices are connected along the flow direction to form the large-scale structure of a streamwise vortex. This behavior provides more details about the formation of large-scale streamwise vortexes. The model explains the characteristics of streamwise vortices (Q2/Q4 events and low- and high-speed streaky structures) and linear streaks based on water depth. Our investigation into the spatial relationships between spanwise vortex density and streaky structures shows that the legs of the hairpin vortex model envelop low-speed streaks. These low-speed hairpin vortex legs can be organized and accumulated into larger-scale quasi-streamwise vortices (Figure 16).

**Figure 16.** Theoretical models of a streamwise vortex, an attached eddy vortex, and a hairpin vortex in the *x–z* plane.

Further analysis of the details of the vortex models led us to propose a simple hypothesis: the three coherent structures, modeled individually as a streamwise vortex, while an attached eddy vortex and a hairpin vortex both exist separately in turbulent flow. It is likely that they are all manifestations of the same turbulent structure under different paradigmatic approaches, as shown in Figure 16.

#### *4.2. Conclusions*

We based this study on 30 sets of closed-channel flow DNS data. Image processing was employed to identify low- and high-speed streaks, using a detection function and a threshold value, *Ct*. Statistical methods were used to calculate the characteristic dimensions of both low- and high-speed streaks. We investigated three models of coherent structures (streamwise vortex, attached eddy vortex, and hairpin vortices) and demonstrated their application. Analysis of the characteristic dimensions of streaky structures and vortices and further analysis of the relationships among the three vortex models led us to suggest a straightforward hypothesis. The results we obtained are summarized as follows.


**Author Contributions:** H.W. performed the data collection, data processing, and manuscript preparation; G.P. assisted with the data analysis. All authors participated in the revision of the manuscript and graphics.

**Funding:** The study is financially supported by the National Natural Science Foundation of China (Grant No. 51709047), the CRSRI Open Research Program (Grant No. CKWV2018461/KY), and the Open Research Fund of State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin (China institute of Water Resources and Hydropower Research, Grant No. IWHR-SKL-KF201815).

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Abbreviations**



#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Viscosity Controls Rapid Infiltration and Drainage, Not the Macropores**

#### **Peter Germann**

Faculty of Science, Institute of Geography, University of Bern, 3012 Bern, Switzerland; pf.germann@bluewin.ch Received: 12 November 2019; Accepted: 20 January 2020; Published: 24 January 2020

**Abstract:** The paper argues that universal approaches to infiltration and drainage in permeable media pivoting around capillarity and leading to dual porosity, non-equilibrium, or preferential flow need to be replaced by a dual process approach. One process has to account for relatively fast infiltration and drainage based on Newton's viscous shear flow, while the other one draws from capillarity and is responsible for storage and relatively slow redistribution of soil water. Already in the second half of the 19th Century were two separate processes postulated, however, Buckingham's and Richards' apparent universal capillarity-based approaches to the flow and storage of water in soils dominated. The paper introduces the basics of Newton's shear flow in permeable media. It then presents experimental applications, and explores the relationships of Newton's shear flow with Darcy's law, Forchheimer's and Richards' equations, and finally extends to the transport of solutes and particles.

**Keywords:** wetting shock fronts; shear flow; viscosity; capillarity; kinematic waves

#### **1. Introduction**

Infiltration is the transgression of liquid water from above the surface of the permeable lithosphere to its interior, while drainage refers to liquid water leaving some of its bulk. For example, within about three weeks after liquid manure applications, bad odors appeared in drinking water wells at depths between 10 to 50 m in English Chalk [1]. In the same geological formations, the annual surplus of the water balance moved down with about 1 m year−<sup>1</sup> as profiles of isotope ratios of 18O/ 16O revealed [2]. The contrasting observations lead to ratios between 150 and 700 of the velocities of the 'odor' front vs. the isotope front that illustrate well the difference between rapid and slow infiltration.

Fast flow and transport are usually attributed to preferential flow in the macropore domain, while slower movements are supposedly due to flow in the capillary domain as expressed with Richards' [3], allowing for exchanges between the two domains as, for instance [4] and [5] have recently compiled. This contribution favors a dual-process approach to infiltration and drainage over dual-porosity approaches, thus avoiding the a-priori delineation between 'macropores' and the remaining porosity.

Beginning in the middle of nineteenth century, the next section summarizes the evolution of permeable media flow concepts for both, saturated and unsaturated flows. The third section provides the base for Newton's shear flow, stressing vertical viscous flow during rapid infiltration, while capillarity is assigned to the much slower redistribution of soil moisture but in all directions. The fourth section provides experimental evidence, and is followed by the conclusions.

#### **2. Evolution of Infiltration and Drainage Concepts**

This section provides some milestones on the way to a dual-process approach to infiltration and drainage, where Newton's shear flow covers fast and gravity driven flows, while Richards' capillary flow is considered to deal mainly with slower capillary rises and redistributions. By no means is this section intended to cover the history of infiltration and drainage.

#### *2.1. Steady Saturated Flow*

In the mid-19th century, there was an increasing interest in flows in saturated soils and similarly permeable media. Hagen, a German hydraulic engineer, and Poiseuille (1846) [6] a French physiologist, independently analyzed laminar flow in thin capillary tubes. Darcy (1856) [7], in the quest of designing a technical filtration system for the city of Dijon, empirically developed the concept of hydraulic conductivity as proportionality factor of flow's linear dependence on the pressure gradient. Dupuit (1863) [8] expanded Darcy's law to two dimensions as perpendicular and radial flow between two parallel drainage ditches and towards a groundwater well, respectively. Forchheimer (1901) [9] added an expression to Darcy's law that accounts for high pressures, steep gradients, and high flow velocities in the liquid.

#### *2.2. Early Studies on Flow in Unsaturated Soils*

Schumacher (1864) [10], a German agronomist, was probably the first who considered capillarity as the cause for simultaneous flows of water and air in partially water-saturated soils. He qualitatively compared the rise of wetting fronts in soil columns with the rise of water in capillary-sized glass tubes, and concluded that the wetting fronts rise higher but slower in finer textured soils compared with coarser materials. He also infiltrated water in columns of undisturbed soil and found that infiltration fronts progressed much faster than the rising wetting fronts. He suggested two separate processes for the two flow types: (i) slower capillary rise and (ii) faster infiltration, however, without further dwelling on the latter. Lawes et al. (1882) [11] concluded from the chemical composition of the drain from large lysimeters at the Rothamsted Research Station that "The drainage water of a soil may thus be of two kinds (1) of rainwater that has passed with but little change in composition down the open channels of the soil; or (2) of the water discharged from the pores of a saturated soil." [11] prioritized two separate flow paths to explain the observations.

#### *2.3. Universal Capillary Flow in Unsaturated Soils*

During the second half of the 19th century, modern irrigation agriculture spread in semi-arid areas, thus increasing the demand for better understanding the soil-water regime. Buckingham (1907) [12], working on a universal approach to the simultaneous storage and flow of water and air in soils, postulated the relationship between the capillary potential ψ Pa and the volumetric water content θ m<sup>3</sup> m<sup>−</sup>3, also known as the water retention function, retention curve, or water release curve. The capillary potential follows from the Young–Laplace (1805) [13] relationship, stating that the pressure difference between a liquid and the adjacent gas phase increases inversely proportional to the radius of the interface. Capillary potential emerges as energy per unit volume of water in the permeable medium due to the water's surface tension. Canceling energy and volume with one length leads to force per area. Because the menisci's surfaces are bent into the liquid, ψ < 0, where ψ = 0 corresponds to the air pressure as reference. In addition to the specific weight of the soil water, [12] introduced the spatial gradient of ψ as the other major driving force. Besides infiltration, this stroke of a genius accounts for the redistribution of soil water in all directions, evaporation across the soil surface, transpiration via roots, and capillary rise from perched water including groundwater tables. In analogy to Fourier's (1822) [14] and Ohm's (1825) [15] laws for heat flow and electrical current, and Darcy's law for water flow in saturated porous media, [12] postulated the hydraulic conductivity for flow in unsaturated porous media as function of either *K*(θ) or *K*(ψ)ms<sup>−</sup>1. According to Or (2018) [16], the British meteorologist Richardson (1922) [17] was most likely the first who introduced a diffusion type of *K*-ψ-θ-relationship in the quest of quantifying water exchange between the atmosphere and the soil as lower boundary of the meteorological system. A second-order partial differential expression became necessary because ψ depends on θ, and both their temporal variations on flow, while flow itself is driven by the gradient of ψ.

The race was on to the experimental determination of the *K*-ψ-θ-relationships. For instance, Gardner et al. (1922) [18] used plates and blocks of fired clay with water-saturated pores fine enough to hydraulically connect the capillary-bound water within soil samples with systems outside them. [3] applied the technique to the construction of tensiometers that directly measure ψ within an approximate range of 0 > ψ >≈ −80 kPa. With the pressure plate apparatus he measured ψ-θ-relationships and determined hydraulic conductivity *K*(ψ or θ). Similar to [17], he presented a diffusion-type approach to the transient water flow in unsaturated soils. Numerous analytical procedures evolved for solving the Richards equation, among them a prominent series of papers by J.R. Philip [19]. Moreover, van Genuchten (1980) [20] developed mathematically closed forms of *K-*ψ*-*θ-relationships that provide the base for the many hues of HYDRUS i.e., various numerical simulation packages dealing with flow and storage of water and solutes in unsaturated soils (Simunek et al., 1999, [21]).

#### *2.4. Early Alternatives to Universal Capillary Flow*

In his quest of demonstrating the benefit of forests and reforestations on controlling floods and debris flows from steep catchments in the Swiss Alps and Pre-Alps, Burger (1922) [22] measured in situ the time lapses Δ*t*<sup>100</sup> for the infiltration of 100 mm of water into soil columns of the same length. In the laboratory, he determined the air capacity *AC* m<sup>3</sup> m−<sup>3</sup> of undisturbed samples collected near the infiltration measurements, where *AC* is the difference of the specific water volume after standardized drainage on a gravel bed and complete saturation.

Veihmeyer (1927) [23] investigated water storage in soils with the aim of scheduling irrigation. He proposed the water contents at field capacity (FC) and at the permanent wilting point (PWP) as upper and lower thresholds of plant-available soil water. Free drainage establishes FC a couple of days after a soil was saturated and evapo-transpiration was prevented. The water loss before achieving FC is referred to as 'drainable' or 'gravitational' soil water. PWP is considered θ at ψ = −150 kPa (= −15 bar).

#### *2.5. Dual Porosity Approaches*

It became unavoidable that concepts based on Buckingham's (1907) [12] fundamental and seminal work contradicted with practical and field-oriented research. Veihmeyer (1954) [24], for instance, stated "Since the distinction between capillary and other 'kinds' of water in soils cannot be made with exactness, obviously a term such as non-capillary porosity cannot be defined precisely since by definition it is determined by the amount of 'capillary' water in the soils". Additionally, progress in field instrumentation as well as in computing techniques allowed for producing and processing large data sets including the numerical solution of Richards equation. In the late 1970s, the development increasingly unveiled substantial discrepancies between measurements and the numerous approaches to water movement in unsaturated soils based on [3] capillarity-dominated theory. Particularly disturbing were observations on wetting fronts advancing much faster than expected from the Richards approach.

The mid-1970s initiated the thread of dual-porosity approaches to preferential flows. Bouma et al. (1977) [25] were among the first to introduce the term macropores in view of non-equilibrium in ψ*-*θ-relationships, while the compilation by Beven and Germann (1982) [26] on the subject is still referred to today. The discussion has gradually moved from macropore flow to preferential flow that summarizes all the flows in unsaturated porous media not obeying the Richards equation [4]. Numerous reviews on macropore flows, preferential flows, non-equilibrium flows, and non-uniform flows in permeable media appear periodically, and only the latest are here referenced [5,27].

Beven (2018) [28] argued that, for about a century, the hardly questioned preference given to capillarity denied recognition of concepts considering flow along macropores, pipes, and cracks. Indeed, there is an increasing number of contributions focusing on the dimensions and shapes of flow paths, their 3-D imaging, and trials to derive flows from them [4]. However, there is hardly an approach capable of applying the wealth of information about the paths to the quantification of flow. Ignoring Veihmeyer's (1954) [23] warning, the attraction of research on flow paths is so dominant that, for instance, Jarvis et al. (2016) [4] flatly denied the applicability of Hagen–Poiseuille concepts to flow in soils. Moreover, advanced techniques of infiltration with non-Newtonian fluids led so far just to the description of path structures rather than more directly to the flow process [29]. Wide-spread research in the types, dimensions, and shapes of 'macropores' and their apparent relationships to flow and transport mostly pivot around Richards equation that is numerically applied to either macropore-/micropore-domains or by modelling flow and solute transport in the macropore domain with separate rules yet still maintaining a Richards-type approach to micropore-flow. Both types of approaches allow for exchanges of flow and solutes between the two domains. Imaging procedures visualize flow in 2-D and 3-D in voids as narrow as some 10 μm, rising hope that the wealth of information gained so far at the hydro-dynamic scale will eventually lead to macroscopic models at the soil profile scale of meters [4]. Thus, Beven's (2018) [28] denial of progress in infiltration research is here carried a step further. The obsession with pores, channels, flow paths, and their connectivity, tortuosity, and necks actually retarded research progress towards more general infiltration that has to be based on hydro-mechanical principles.

#### *2.6. Early Search for Alternatives*

The alternative approach should be based on the same principles as Hagen–Poiseuille's and Darcy's laws, whereas Newton's (1729) [30] shear flow (Nsf) appears as a solid and suitable foundation. The approach then should close the gap of one to two orders of magnitude of hydraulic conductivity between saturated flow and flow close to saturation [31].

Unearthing Burger's (1922) [22] data, [32] found an encouraging coefficient of determination of *r*<sup>2</sup> = 0.77 when correlating via a Hagen–Poiseuille approach 76 pairs of Δ*t*100- and *AC*-values. Beven and Germann (1981) [33] modelled laminar flow in tubes and planar cracks, and proposed kinematic wave theory according to Lighthill and Whitham (1955) [34] as analytical tool for handling Nsf. Germann (1985) [35] applied the theory successfully to data from an infiltration-drainage experiment carried out on a block of polyester consolidated coarse sand that experimentally support the flow-(*q*-) version of Nsf, where *q* m s−<sup>1</sup> is volume flux density.

Dye experiments confirmed expectations, not the least by purposefully setting the initial and boundary conditions. [36], for instance, demonstrated preferential flow along earthworm burrows down to the 0.4 m depth as well as radially away from the channels by inundating with rhodamine dye and a bromide solution the tops of columns of undisturbed soil. Following [11], who reported fast drainage, Germann (1986) [37] assessed the arrival times of precipitation fronts in the Coshocton lysimeters. Accordingly, rains of 10 mm d−<sup>1</sup> sufficed to initiate or increase drainage flow within 24 h at the 2.4 m depth when θ in the upper 1.0 m of the soil was at or above 0.3 m3 m−3. The observations result in wetting front velocities greater than 3 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m s<sup>−</sup>1.

#### *2.7. Development of Newton's Shear Flow Approach to Infiltration and Drainage*

The advent of TDR- and data recording techniques fostered the water-content-(*w*-) version of Nsf (where *w* m3 m−<sup>3</sup> refers to the mobile water content) applicable to in-situ infiltration experiments using controlled sprinkler irrigation and producing time series of θ at various depths in soil profiles. Discussions with L. Dipietro and V. P. Singh linked experiments with basic concepts [38]. Investigations of acoustic velocities across a column of an undisturbed soil during infiltrations, [39] proved independently from the capillarity narrative that rapidly infiltrating water is under atmospheric pressure. Originally, macropore flow was the prominent research topic. Shortfalls of the capillarity-dominated concepts emerged when juxtaposed to matured Nsf-versions [40]. Hence, this paper is the première of consequently carrying through the dual-process approach to infiltration and drainage.

There were projects leading off the track. [41], for instance, proposed a kinematic wave approach to infiltration that included a sink function accounting for water abstraction from the moving to the sessile parts. Later rigorous testing of the approach against real data unveiled the misconception. Additionally, [42] superimposed a bundle of rivulets to model data. That again turned out as impracticable. [43] were looking for FC in soil profiles after the cessation of infiltration. However, repeated experiments showed quite variations in expected constant final FCs.

#### *2.8. State of Nsf*

Gradually, some of the subordinate concepts lined out above get referenced. However, as of yet, outside the author's immediate environment of research, there is hardly an elaborated application of Nsf to infiltration and drainage. Unavoidably, this leads to the impression of self-indulgence, particularly in view of the paper's lopsided list of references. The Guest Editor's invitation to its submission may ease the charge. Moreover, although [26] mainly outlined the need for research on 'macropore flow' and did not offer much of research guidance, that paper is still frequently referenced, indicating that researchers have hardly agreed on a basic approach to the problem as they did and still do, for instance, on the Richards equation.

#### **3. Newton's Shear Flow in Permeable Media**

The section summarizes the basics of Newton's shear flow (Nsf) in permeable media. Detailed derivations of the relationships have previously been published [44–46], thus only the major expressions are here presented. The approach is laid out at the hydro-mechanical scale of spatio-temporal process integration, allowing for its easy handling with analytical expressions, yet under strict observance of the balances of energy, momentum, and mass (i.e., the continuity requirements).

The interior of a permeable solid medium contains connected flow paths that are wide enough to let liquids pass across its considered volume. The definition purposefully avoids further specification of the flow paths' shapes and dimensions. Water supply to the surface is thought of a pulse *P*(*qS*,*TB*,*TE*), where *qS* m s−<sup>1</sup> is constant volume flux density from the pulse's beginning at *TB* to its ending at *TE*, both s. (The subscript *S* refers to the surface of the permeable medium). The pulse initiates a water content wave (*WCW*) of mobile water that is conceptualized as a film gliding down the paths of a permeable medium according to the rules of Newton's shear flow. The *WCW* is the basic unit of shear flow whose hydro-mechanical properties are going to be presented below.

Referring to Figure 1, the parameters film thickness *F* m and specific contact length *L* m m−<sup>2</sup> per unit cross-sectional area *A* m2 of the medium specify a *WCW*. Regardless of the thickness of *F*, atmospheric pressure prevails within the film. A *WCW* supposedly runs along the flow paths while forming a discontinuous and sharp wetting shock front at *zW*(*t*). The *WCW* partially fills the voids of the upper part of the medium within 0 <sup>≤</sup> *<sup>z</sup>* <sup>≤</sup> *zW*(*t*) with the mobile water content *<sup>w</sup>*(*z*,*t*) m<sup>3</sup> <sup>m</sup><sup>−</sup>3, where *w* <sup>&</sup>lt; <sup>ε</sup> <sup>−</sup> <sup>θ</sup>ante, where <sup>ε</sup> and <sup>θ</sup>ante are porosity and antecedent <sup>θ</sup>, respectively, both m<sup>3</sup> m<sup>−</sup>3. The lower part *z* > *zW*(*t*) remains at θante. The coordinate *z* m originates at the surface and points positively down. The water film of the *WCW* consists of an assembly of parallel laminae, each d*f* m thick and moving with the celerity of *c*(*f*)ms−1. Celerity refers to the wave velocity, for instance of a lamina. Newton (1729) [30] defined viscosity as "The resistance, arising from the want of lubricity in the parts of a fluid, is, caeteris paribus, proportional to the velocity with which the parts of the fluid are separated from each other." In our case, the definition translates to the shear stress ϕ(*f*) Pa in the unit area of *L* × *A* <sup>×</sup> *zW*(*t*) m<sup>2</sup> per unit volume *A* <sup>×</sup> *zw*(*t*) m<sup>3</sup> of the medium at distance 0 <sup>≤</sup> *f* <sup>≤</sup> *F* m from the soil-water interface (*SWI*), acting in the direction opposite to gravity. ϕ(*f*) balances the weight of the film from *f* to *F*. Integration of the shear stress balance from 0 to *f* under consideration of the non-slip condition of *c*(0) = 0 leads to the parabolic celerity profile of the laminae within the film. The differential flow of a lamina is d*q*(*f*) = *L* × *c*(*f*) × d(*f*). Its integration from 0 to *F* i.e., from the *SWI* to the air-water interface (*AWI*), leads to the volume flux density of the entire film as:

$$q(F, L) = \frac{\mathcal{G}}{\mathfrak{d} \cdot \eta} \cdot L \cdot F^3 \tag{1}$$

m s<sup>−</sup>1, where *<sup>g</sup>* (=9.81 m s<sup>−</sup>2) is acceleration due to gravity and <sup>η</sup> (≈10−<sup>6</sup> <sup>m</sup><sup>2</sup> <sup>s</sup><sup>−</sup>1) is temperature-dependent kinematic viscosity. The volume of mobile water per unit volume of the permeable medium from the surface to *zW*(*t*) amounts to:

$$w(F, L) = F \cdot L \tag{2}$$

m<sup>3</sup> m<sup>−</sup>3. The constant velocity of the wetting shock front follows from the volume balance as:

$$w\_W(F) = \frac{q(F, L)}{w(F, L)} = \frac{\mathcal{g}}{\mathfrak{d} \cdot \eta} \cdot F^2 \tag{3}$$

**Figure 1.** Film flow along a vertical plane, where *F* m is the film thickness, *f* m is the thickness variable, d*f* m is the thickness of a lamina, *zW*(*t*) m is the vertical position of the wetting shock front as function of time *t* s, and *L* m m−<sup>2</sup> is the specific contact length between the moving water film and the sessile parts of the permeable medium per its unit horizontal cross-sectional area *A* m2. *AWI* and *SWI* are the air-water and the solid-water interfaces, respectively ([44]; with permission from the publisher).

Thus, the position of the wetting shock front as function of time becomes:

$$
\omega\_W(t) = \upsilon\_W \cdot (t - T\_B) = (t - T\_B) \cdot \frac{\mathcal{G}}{\mathfrak{d} \cdot \eta} \cdot F^2 \tag{4}
$$

m s−1. Accordingly, *L* m m−<sup>2</sup> is not only the specific contact length of the film per *A* but, under consideration of *zW*(*t*), *L* m<sup>2</sup> m−<sup>3</sup> evolves as the specific vertical contact area of the *WCW* with the sessile parts of the permeable medium per unit volume. It thus turns into the locus of exchange of momentum, heat, capillary potential, water, solutes, and particles between the *WCW* and the sessile parts of the medium.

*Water* **2020**, *12*, 337

Equations (1)–(4) hold during infiltration i.e., *TB* ≤ *t* ≤ *TE*. Input ends abruptly at *TE* and at *z* = 0 i.e., *qS*→0, when and where the *WCW* collapses from *f* = *F* to *f* = 0. All the rear ends of the laminae are released at once at *z* = 0. The outermost lamina moves the fastest with the celerity of the draining front as:

$$\mathbf{c}\_{D} = \frac{\mathbf{d}q(F)}{\mathbf{d}w} = \frac{\mathbf{d}q(F)}{L \cdot \mathbf{d}f} = \frac{\mathbf{g}}{\eta} \cdot F^2 = \mathbf{3} \cdot \mathbf{v}\_{W} \tag{5}$$

m s<sup>−</sup>1. While *vW* follows from the volume balance Equation (3), the three times faster *cD* follows from the parabolic function *c*(*f*) (see, for instance, Germann and Karlen, 2016). The outer most lamina moves with *cD* also during 0 ≤ *t* ≤ *TE*. Therefore, the slower moving wetting shock front continuously intercepts the faster moving laminae. This results in a sharp wetting shock front at *zW*(*t*) with θante + *w* between the surface and the wetting shock front and θante ahead of it. The wetting shock front presents a discontinuity of the *WCW*. The rear end of the draining front that was released at *TE* catches up with wetting shock front at time *TI* s that follows from setting *vW* × (*TI* − *TB*) = *cD* × (*TI* − *TE*) as

$$T\_I = \frac{1}{2} \cdot (3 \cdot T\_E - T\_B) \tag{6}$$

Thus, *TI* is exclusively an expression of the pulse duration. The wetting shock front intercepts the draining front at depth:

$$Z\_I = (T\_E - T\_B) \cdot F^2 \cdot \frac{\mathcal{S}}{2 \cdot \eta} \tag{7}$$

m. The two expressions *TI* and *ZI* represent the spatio-temporal scale of a *WCW*. The rear ends of all the other laminae move with decreased celebrities, ultimately leading to the spatio-temporal distribution of the mobile water content from the surface to the wetting shock front as:

$$w(z,t) = L \cdot \left(\frac{\eta}{g}\right)^{1/2} \cdot z^{1/2} \cdot \left(t - T\_E\right)^{-1/2} \tag{8}$$

After *TI* and beyond *ZI* the draining front disappears and *vW*(*z*,*t*) decreases with time and depth, leading to the position of the wetting shock front as:

$$z\_W(t) = \left(\frac{\mathfrak{Z} \cdot V\_{\rm WCW}}{\mathfrak{Z} \cdot L}\right)^{2/3} \cdot \left(\frac{\mathfrak{Z}}{\eta}\right)^{1/3} \cdot (t - T\_E)^{1/3} \tag{9}$$

with the mobile water content of:

$$\left.w(t)\right|\_{z\_W} = \left(\frac{\eta}{\mathcal{S}}\right)^{1/3} \cdot \left(\frac{3 \cdot V\_{\rm WCW}}{2}\right)^{1/3} \cdot \left(t - T\_E\right)^{-1/3} \cdot L^{2/3} \tag{10}$$

while the volume flux density at the wetting shock front becomes:

$$\left.q(t)\right|\_{z\_W} = \frac{V\_{\rm WC}}{2 \cdot (t - T\_E)}\tag{11}$$

where *VWCW* = *qS* × (*TE* − *TB*) m is the total volume of the *WCW*. Figure 2 illustrates the relationships describing a *WCW*.

**Figure 2.** Schematic representation of a water-content wave *WCW*, where the *w(z,t)-*axis represents the mobile water content, *t* and *z* are the axes of time and depth; *wS* indicates the mobile water content that follows from *qS*; *TB* and *TE* s are the beginning and ending times of the water pulse *P(qS*, *TB*, *TE)* that hits the surface at *z* = 0; *TI* and *ZI* indicate time and depth of the wetting front intercepting the draining front. The wetting shock front continues beyond *t* > *TI* and *z* > *ZI* as curve along time and depth. ([44]; with permission from the publisher).

A *WCW* i.e., *w*(*z*,*t*),translates easily into a volume flux density wave according to:

$$\left|\boldsymbol{q}(\boldsymbol{z},t) = \boldsymbol{b} \cdot \boldsymbol{w}(\boldsymbol{z},t)\right|^3 \tag{12}$$

with *b* = *g*/ <sup>3</sup> · *<sup>L</sup>*<sup>2</sup> · <sup>η</sup> m s−1. The approach applies to laminar flow as assessed with the Reynolds number:

$$\text{Re} = \frac{F \cdot v}{\eta} = \frac{F^3 \cdot g}{3 \cdot \eta^2} = \left(\frac{3 \cdot v^3}{g \cdot \eta}\right)^{1/2} \tag{13}$$

*Re* ≤ 1 strictly defines laminar flow; however, depending on the application, *Re* > 1 might be tolerable, yet within an undisclosed range.

The parameters film thickness *F* m and specific contact area *L* m2 m−<sup>3</sup> of a *WCW* due to a pulse *P*(*qS*,*TB*,*TE*) follow from time series of either volumetric water contents θ(*Z*,*t*) or from volume flux density *q*(*Z*,*t*), where *Z* (m) represents a specific depth in the permeable medium. Presumably, *F* and *L* get established shortly after *TB*. Their sizes are due to *qS* and the actual properties of the permeable medium, such as θante and related flow paths supporting *F* and *L*. Predictions of *F* and *L* are currently undisclosed, thus restricting the approach to *a posteriori* data interpretation.

#### **4. Examples**

The section provides examples in support and as illustrations of Nsfa. More, and more elaborated examples are presented in [44].

#### *4.1. Experimental Support of Newton's Shear Flow in Permeable Media*

Constant wetting front velocities, Equation (3): Already [47] observed constant velocities of wetting fronts during finger flow in sandboxes down to the 0.9-m depth. [48] reported the same kind of observations over a depth range of 2 m during infiltration in the Kiel Sand Tank, while [49] found constant wetting front velocities under natural precipitation regimes across the upper 1.5 m of the undisturbed soil in a large weighing lysimeter. With 34 FTDR-probes in three boreholes, [50] followed infiltration in an ancient sand dune in Israel. From their data, [44] deduced two constant wetting front velocities: one of about 5 <sup>×</sup> 10−<sup>6</sup> m s−<sup>1</sup> above the clay-rich layer at the 8-m depth, and one of about <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>6</sup> m s−<sup>1</sup> down to the 21 m depth. The two coefficients of determination (*r*2) of the linear regressions of arrivals vs. depths of the wetting front amounted to 0.83 and 0.94, respectively. Constant wetting front velocities support the early stage of *WCW*s, Equation (3), and imply constant *F*.

Laminar flow: From more than 200 infiltration experiments, [51] reported a range of 10−<sup>5</sup> < *vW* < <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>2</sup> m s<sup>−</sup>1. According to Equation (13), this leads to a range of 10−<sup>4</sup> <sup>&</sup>lt; *Re* <sup>&</sup>lt; 1.6. The upper limit is still close to laminar flow. In General, permissible upper limits of *Re* demand attention in view of the applicability of Nsf.

Atmospheric pressure behind wetting fronts: Atmospheric pressure in the *WCW* is a prerequisite of Nsf-theory. [48] antecedent capillary heads in the range from −1.2 to −0.8 m to collapse to about −0.2 m behind the wetting front. [42] found similar patterns of capillary-head collapses during infiltration into a column of a sandy-loam textured mollic Cambisol. Additionally, acoustic velocities during infiltration in a column of an undisturbed loess soil, [39] led [40] to the same conclusion of capillary head's break-up close to atmospheric pressure after the passing of the wetting front.

Cohesion of Newton's shear flow approach*:* The parameters *F* and *L* suffice to treat infiltration and drainage with Nsf, Equations (1) and (12). In principle, time series of either θ(*Z*,*t*) or *q*(*Z*,*t*) permit calibration of the two parameters: [45] shows their calibration from θ(*Z*,*t*), while [47] provide the procedure from *q*-series. The latter authors and [44] also demonstrate the cohesion of the approach in that derivations from both time series resulted in the same amount of *L*.

#### *4.2. Newton's Shear Flow vis-à-vis Darcy's Law and Forchheimer's Equation*

Darcy's law mutates to the saturated-flow extension of unsaturated vertical shear flow under the following considerations. When replacing kinematic viscosity with the dynamic viscosity, μ = ρ × η, follows from Equation (1) that:

$$
\rho \theta\_{\text{ante}} + w < \varepsilon \text{ and } \Delta p / \Delta z = \rho \times \text{g} \qquad q = \frac{F^3 L}{3\,\mu} \rho \text{ g} \tag{14}
$$

where θante m3m−<sup>3</sup> is the antecedent volumetric water content, ε m3m−<sup>3</sup> is porosity, Δ*p*/Δ*z* Pa m−<sup>1</sup> is the pressure gradient, ρ (= 1000 kg m<sup>−</sup>3) is the density of water, and μ = ρ η Pa s is dynamic viscosity. At saturation we get:

$$\rho g\_{\rm ante} + w = \varepsilon \text{ and } \Delta p / \Delta z = \rho \times \text{g} : \qquad q\_{\rm sat} = \frac{F\_{\rm sat}^3 L\_{\rm sat}}{3\mu} \rho g = K\_{\rm sat} \tag{15}$$

while an external pressure gradient leads to:

$$
\theta\_{\rm anle} + w = \varepsilon \text{ and } \Delta p/\Delta z > \rho \times \text{g}: \qquad q(p) = \frac{F\_{\rm sat}^3 L\_{\rm sat}}{3\mu} \frac{\Delta p}{\Delta z} = q\_{\rm sat} \frac{\Delta p}{\Delta z \rho g} \tag{16}
$$

Darcy's law states that *q* is proportional to Δ*p*/Δ*z* i.e., volume flux density is a linear function of the flow-driving gradient with the proportionality factor *Ksat*. In view of the various dimensionalities of *w* proportional to (*L*1, *F*1), *v* proportional to (*L*0, *F*2), and *q* proportional to (*L*1, *F*3), linearity seems only possible if *Fsat* and *Lsat* remain constant and independent from *p* in the transition from only gravity-driven to pressure-driven shear flow at saturation i.e., in the transition from Equation (15) to Equation (16). The elaboration supports the linearity of Darcy's law, but it is not an independent proof of the law's linearity. As a consequence, *w* = *q*/*v* also remains constant. Further, if θante + *w*

= ε, d*Lsa*t/d*p* = 0 and d*Fsa*t/d*p* = 0 then follows the hypothesis that (*Fsat* × *Lsat*) represent (*F* × *L*)*ma*<sup>x</sup> leading to *Ksat*. However, other combinations of (*F* × *L*) in unsaturated media are feasible that may lead to *qunsat* > *qsat* = *Ksat*. The presumption of *qunsat* > *Ksat* opens an new view on shear flow and the basics of infiltration that are in stark contrast to Richards' capillary flow, where a priori *Ksat* > *K*(θ or ψ). Comparisons of the rates from field infiltration with *Ksat*-measurements in the laboratory vaguely support the presumption of [45] and [52] most recently provide experimental evidence of *qunsat* > *Ksat*.

Atmospheric pressure prevails in a *WCW* in Newton's shear flow. However, parts of the the *WCW* may hit flow paths narrower than *F*, and *p* > *patm* will occur. The difference *p* > *patm* will dissipate almost instantaneously with the acoustic velocity in soil water between 300 and 800 m s−<sup>1</sup> [53] if virtually no water has to be moved over short distances to the next path wider than *F*. However, increasing the distances between the paths that are wider than *F* leads to local pressure build up that may eventually lead to sessile water above a layer of reduced conductivity. Water saturation may occur from relatively short periods, for instance shortly after heavy rains, up to seasons as, for instance, mottles of chemically reduced and oxidized zones in soil profiles may illustrate. Drainage under *p* > *patm* shows as a steeper recession than drainage under *p* = *patm*. [49] provide the details. [54] report θ(*t*)-time series of TDR-measurements immediately above an impervious soil layer that deviates from those of free flow during infiltration and that indicates the evolution of sessile soil water.

Forchheimer (1901) [9] added a quadratic term to Darcy's law that accounts for momentum dissipation under high flow velocities that occur under high pressure gradients i.e.,

$$-\frac{dp}{d\mathbf{x}} = \frac{\eta}{K}\mathbf{v} + \frac{\rho}{k\_2}\mathbf{v}^2\tag{17}$$

where *K* and *k*<sup>2</sup> represent hydraulic conductivity and turbulence, respectivly. Neither of the two conditions apply to Nsf in unsaturated permeable media in the lithosphere near its interface with the atmosphere. Moreover, [54] estimate the relative contribution of kinetic energy in the range of 10−<sup>6</sup> to 10−<sup>4</sup> in comparison with viscous momentum dissipation during flow in soils. Thus, Forchheimer's equation does not apply to Nsf.

#### *4.3. Gravity and Capillarity*

Gravity is the unique driving force in Newton's shear flow in unsaturated permeable media, while capillarity is responsible for the soil water's redistribution and vertical rise from saturated zones. [10] suggested a two-process approach to water flow and storage in partially saturated permeable media. While he recognized capillarity as responsible for the water's rise, and probably also its contribution to water redistribution in soil columns, he left open the mechanism behind infiltration. Here, the focus is on infiltration that is completely gravity-driven and viscosity-controlled, yet allowing for water abstraction due to capillarity from the mobile to the immobile part of the permeable system. Richards' universal equation for flow and storage in unsaturated porous media revolves around capillarity. It prominently, and unnecessarily, reduces the degrees of freedom of flow [39]. Unjustified, it relates the coefficient of momentum dissipation with antecedent flow conditions i.e., *K*(θ or ψ), thus a priori excluding atmospheric pressure in the moving water. The exclusion leads to too slow advancements of wetting as Germann and Hensel (2006) [55] have demonstrated by comparing the results from HYDRUS-2 model [21] performances with observed infiltrations from more than 200 sites. Concentrating on gravity and viscosity liberates infiltration and drainage from the omnipresence of capillarity in soil hydrology with the benefit of avoiding the difficult definitions of non-equilibrium flow and the separation of macropores from the remaining pores. With respect to capillarity, the relative contribution of gravity to flow varies according to *cos* (α), where α◦ is the angle of deviation from the vertical. Thus, at *cos* (0◦) = 1, as in the cases presented above, gravity's contribution is at maximum; it reduces to *cos* (90◦) = *cos* (270◦) = 0, while it completely opposes capillarity at *cos* (180◦) = −1. The juxtaposition illustrates the spatial limitations of Nsf. Moreover, lateral rapid flow requires saturated conditions along a layer with path widths narrower than *F.*

Pressure in the *WCW* is atmospheric while ψ < 0 typically prevails ahead of it. Therefore, water is abstracted from the *WCW* onto *L.* Abstraction is usually completed during short periods. The amount of abstraction shows in the difference between θend, when the *WCW* has ceased, and θante before the passing of the *WCW*. θend is approached with Equation (8) for a particular depth *Z*, while setting the *WCW*'s end at *w*(*Z*,*t*)/*wS* = *u*, where *u* << 1 is an arbitrarily selected threshold that depends on the particular application of Newton's shear flow. Water abstracted from the *WCW* is thus available for redistribution and uptake by plant roots.

#### *4.4. Scales*

The process scale during Nsf with constant *F* and *L* follows from the time and depth of interception, *TI* and *ZI*, Equations (6) and (7) that depend on the duration (*TE* − *TB*) of the input and on *F*. The velocity of the wetting shock front reduces after *TI* and beyond *ZI* according to Equation (9) [48].

The system-related scales of the applicability of Nsf range from a couple of sand grains at the mm-scale to the km-scale. From neutron radiographs taken during infiltration in sand boxes, Hincapié and Germann (2010) [56] calculated fluxes in layers that were about 1 mm thick and volume balances in the cm3-range. Flow across 21 vertical m of an ancient sand dune was already presented above [49], while tracer experiment across crystalline rocks suggests the applicability of Newton's shear flow approach to the km-range. Dubois (1991) [57] injected the tracers uranin and eosin about 1800 m above the Mont Blanc car tunnel that connects Chamonix (France) with Courmayeur (Italy). Within 108 d he detected the tracer fronts in seeps in the car tunnel. This amounts to *vW* <sup>≈</sup> 2 <sup>×</sup> 10−<sup>4</sup> m s−<sup>1</sup> that is well within the range of the *vW*-collection reported above. Moreover, the gentle water seeps in the tunnel revealed that gravity was consumed by viscosity during the long-lasting flow.

The observations of [57] across 1800 m of crystalline rocks of the Mont Blanc massif and the water balance calculations of finger flow in the sand boxes of [56] at the scale of millimeters hint at the spatial tolerance of Newton's shear flow. This may advance the approach to an attractive tool, for instance, for the study of infiltration into groundwater systems.

#### *4.5. Time-Variable Infiltration*

So far input to the surface is considered a single pulse *P*(*qS*, *TB*, *TE*). Time variable input needs to be separated into a series of rectangular pulses each carrying its individual parameters. Under the assumption of the macropore flow restriction i.e., d*L*/d*q* = 0, a series of pulses can be routed as a sequence of kinematic waves according [30], whose mathematical approach models Newton's shear flow correctly if *a* = 3 in Equation (12). [58] provides details of applying characteristics to the multi-pulse infiltration, while [49] lists some results.

Upon observations of infiltrations with varying durations and the rates from 5, to 10, 20, and 40 mm h−<sup>1</sup> into a column of an undisturbed soil [59], postulated the macropore flow restriction of d*L*/d*qS* = 0, meaning that flow occurs along the same paths independently from the input rate. They tested the relationship of *vW* (*qS*) <sup>=</sup> *<sup>b</sup>*1/<sup>3</sup> <sup>×</sup> *qS* <sup>2</sup>/3, Equations (1) to Equation (3), that resulted in an acceptable coefficient of determination of *r*<sup>2</sup> = 0.95. However, that was an exception hardly achieved again. This calls for further investigation of the relationship, most likely under consideration of θante.

#### *4.6. Transport of Tracers and Particles*

Preferential flow in soil hydrology is frequently associated with enhanced and accelerated particle, solute and, primarily, pollutant breakthrough [60]. However, Bogner and Germann (2019) [61] report considerable delays of tracer breakthrough compared with the first arrival of wetting shock fronts at the drain of soil columns with heights of 0.4 m. They referred to the phenomenon as 'pushing out old water' that is well known in catchment hydrology. They statistically explained 81% of the observed delay variations with combinations of *L* and *F* when applying Newton's shear flow to the data. It appears that tracer exchange on large *L* from thin *F* of the *WCW* may be even faster than presumed 'preferential' tracer transport. Under consideration of the mechanistic parameters *F* and *L*, Newton's shear flow provides for a novel tool for the unambiguous investigation of tracer transport and exchange i.e., accelerated as well as decelerated breakthroughs that primarily addresses the tracer's mass balance and secondarily diffusion.

Similar considerations may apply to the transport of particles and microbes. Based on a precursor of Nsf, [62] approached with drag-forces the transport of latex beads and bacteriophages in soils. The three types of electrically uncharged latex beads had dimeters of 0.5, 1.0, and 1.75 μm. Maximal particle concentrations relative to the input suspensions reached about 0.003. While the relationships between relative particle concentrations vs. drag force showed intermediate to strong linearity during the increasing and decreasing limbs of the drain hydrograph, the drag force in the trailing was apparently too small to produce a considerable effect on their translocations. The five bacteriophages H40/l, H4/4, T7, H6/l, and fl carried ζ-potentials in the range of −50.1 to −31.7 mV. Their much weaker relative concentrations vs. drag force showed pronounced, yet undisclosed hysteresis. Recently, [63] called for improving our understanding of the transport of microplastic particles in soils, while [64] found these particles in worm burrows and assumed preferential flow along macropores as the major transport process. Thus, it might be worthwhile to further develop Newton's shear flow application to the particle transport under particular attention of the films thickness *F* and their specific contact areas *L*.

#### *4.7. Ecohydrology*

The early 20th century saw numerous runoff studies in headwater catchments, many of them with the purpose of demonstrating the advantage of forests over open lands in the mitigation of flood and debris flows from steep slopes, mainly to convince governments for subsidizing huge reforestation projects. Burger's (1922) [22] investigations with soil cores as well as his comparative runoff measurements from the 50% forested watershed of Rappengraben vs. the completely forested catchment of Sperbelgraben are well known in forest hydrology. However, not all the compartments from precipitation to runoff got the necessary attention in research they would require for closing the water balance. For instance, the infiltration experiments of [22] were not carried through to assess with them drainage and runoff in the headwater catchments.

Issues of reforestations' remedial effects on poorly permeable clay soils at sites in the Swiss pre-alps have surfaced again because numerous sites require rejuvenation of the more than 100 y old stands. In this context, Lange et al. (2009) [65] were able to demonstrate positive effects of tree root density on rapid infiltration in stagnic soils.

Wiekenkamp et al. (2019) [66] compared in the Wüstebach (BRD) catchment the hydrology of a forest site with the neighboring clear-cut site. They concluded that infiltration via preferential flow paths increased after deforestation. However, their purposeful exclusion of interception may deprive them of further interpreting the difference between the two treatments in rapid infiltration because interception may turn out as major gate controlling the input pulses *P*(*qS*,*TB*,*TE*).

#### **5. Summary and Conclusions**

The first section reports from English chalk rapid contaminant transport about 150 to 700 faster than the slow tracer movements. The second section deals with the evolution of infiltration-drainage concepts since the mid-19th Century. Two lines of thought emerge: On the one hand, Schumacher (1864) [10] suggested a two-processes approach for flow and storage of water in soils, one accounting for infiltration, the other one for capillary flow. On the other hand, a dual-porosity approach follows from the observations of Lawes et al. (1882) [11]. In searching an universal approach for storage and flow of water in partially saturated soils, Buckingham (1907) [12] and Richards (1931) [3] focused on capillarity and the soil hydrological functions *K-*θ*-*ψ. Concepts of non-equilibrium flow emerged upon discovered discrepancies between theory and observations that led later on to dual-porosity approaches. The third section takes off from Schumacher's two-processes approach, exploring fast infiltration and drainage as gravity-driven and viscosity controlled flow, resulting in a water content

wave, *WCW*, that is based on Newton's shear flow i.e., laminar flow in unsaturated permeable media lasting two to ten times longer than the duration of input and occurring under atmospheric pressure. This restriction suggests strong capillary gradients between the *WCW* and the sessile parts. Because of the wide specific contact area, expressed with the factor *L* in the typical range of 500 to 20,000 m<sup>2</sup> m<sup>−</sup>3, and the thin films *F* that are typically between 2 and 50 μm thick, exchanges of water, energy, heat, solutes, and particles are fast and strong. The fourth section elucidates applications of the approach.

On the one side, the universal Richards equation deals with infiltration and redistribution of soil water. It assumes that entire soil moisture θ participates in the flow process, thus inevitably leading to apparent conditions of non-equilibrium, while dual-porosity approaches are hoped to lead away from the uneasy situation. On the other side, Newton's shear flow focuses on viscosity and gravity resulting in adequate time scales during infiltration and associated drainage. Once the wetting shock fronts have ceased to advance, capillarity takes over and soil moisture redistributes towards ψ*-*θ equilibria. This leads to a dual-process approach to infiltration and redistribution, where time and depth of front interception, *TI* and *ZI*, serve as scales for separating the two processes. A further consequence of concentrating flow and transport on Newton's shear flow are the spatio-temporal limits of the processes, expressed with a few multiples of *TI* and *ZI*. Moreover,' pushing out old water' and assuming flow rates under atmospheric pressure exceeding *Ksat* indicate novel aspects associated with Newton's infiltration that were not considerable in previous approaches to preferential flows. Finally, the analytical expressions are amenable to mathematical procedures, such as kinematic wave theory, and their theoretical combinations may lead to new and solid hypotheses calling for experimental testing.

**Funding:** This research did not receive any external funding

**Acknowledgments:** I thank Scholle, Guest Editor of the Special Issue on "Physical and Mathematical Fluid Mechanics", for inviting the manuscript to WATER. I deeply appreciate the positive critiques of three anonymous reviewers that helped considerably to straighten out the manuscript.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
