Characterization of preferential flow in soils near Zarqa river (Jordan) using in situ tension infiltrometer measurements

Study site

The study area between the KTP and KTD (Fig. 1) contains forage and vegetables, as well as olive and citrus orchards. The farmland is 300–500 m above sea level. The relatively warm weather, characterized by a maximum temperature range of 25 to 46 °C and an average temperature range of 19 °C, enables farmers to undertake 2 or 3 cropping cycles per year. Farmers follow conventional tillage practices, using a moldboard/chisel plow to dig to a depth of 10–15 cm and, in some fields, harrowing the topsoil. Farmers usually till their fields before each cropping cycle, typically during the fall (September/October) or early spring (March/April).

The KTP is used to irrigate the agriculture lands adjacent to the main reach of the ZR. The average annual rainfall depth of less than 150 mm constitutes a small fraction of total crop water requirements. Hence, agricultural practices depend mainly on irrigation water supplied directly from the river or local wells.

Four types of soils exist within the study area (Fig. 2). These are Jerash (JER), Abu Salih (ALI), Aydoun (AYD), and Ramtha (THA). These soils belong to the Lithic and Typic Xerorthents, and Xerochreptic Camborthids subgroups (Table 1). The JER, ALI, and AYD overlay geologic formations that consist mainly of limestone and chalky limestone. Meanwhile, the THA sits on a basalt formation (Al-Mahamid, 2005).

Tension infiltrometer measurements

A standard tension infiltrometer equipped with an infiltration disc that was 20 cm in diameter, a pressure transducer, and a data logger (supplied by ICT International) similar to the device described by Ankeny, Kaspar & Horton (1988) were used. A total of 69 measurements were taken from May through September of 2016. The measurements were made during the summer growing season, at least one month after seedbed preparation and plowing activity. At each location, the upper layer of loose soil and plant residue was removed. The soil surface was then leveled and a thin layer of silica sand was placed under the infiltrometer disc to maintain hydraulic contact between the infiltration disc and the soil surface. The measurements were taken sequentially at three tensions of −6, −3, and −0 cm. The −0 indicates that the tension infiltrometer was adjusted to a slight tension to prevent the free flow of water from the tension disc. The water intake was recorded at intervals of one minute. The measurements were allowed to proceed until a steady infiltration rate had been attained, which took about 20 to 60 min for each tension. Thus, the tension infiltrometer measurements can be interpreted according to the following steady-state flow equation (Wooding, 1968): (1)${\displaystyle Q\left(h_i\right)=\pi r^2{K}_{s}exp\left(\alpha h_i\right)\left(1+\frac{4}{\pi r\alpha }\right)}$ where Q(h_i) is the flow rate at tension h_i(L³∕T), K_s is the saturated hydraulic conductivity (L∕T), α is the sorptive number (1∕L), and r is the radius of the tension infiltrometer disc (L).

Considering two data points at h₁ =  − 6 cm and h₂ =  − 3 cm, and taking K_s as a constant, then dividing the Q(h₂) equation by Q(h₁) and solving for α: (2)${\displaystyle \propto =\frac{ln\left[Q\left(h_2\right)∕Q\left(h_1\right)\right]}{h_2-h_1}.}$ The −6 to −3 tension range was used for the calculation of α instead of the −3 to −0 tension range because the infiltration rate and K_s increase rapidly as the tension approaches zero. This may lead to the overestimation of α values if they were derived from a tension range close to saturation (Clothier & Smettem, 1990; Jarvis & Messing, 1995).

The first measurement location was near the KTP, while additional sites were selected along the westward course of the river (Fig. 2). However, the locations were not necessarily visited sequentially; some back-and-forth occurred due to issues with accessibility and with navigating the terrain of the ZR. The geographic coordinates were recorded using a Garmin Etrex 30 GPS device.

Susceptibility to preferential macropore flow

The physical characterization of preferential flow is generally related to macroporosity, defined as the ratio of the total volume of macropores to the total volume of soil (e.g., Beven & Germann, 1982; Luo, Lin & Schmidt, 2010; Mossadeghi-Björklund et al., 2016). Indeed, this is a beneficial parameter. However, the definition of macropores is rather arbitrary. For example, Beven & Germann (1982) reported that values between 30 to 10,000 µm represented the effective macropore diameter, or the equivalent of −102 to 1.02 cm tension. Jarvis (2007) indicated that pores with diameters larger than 300 µm can conduct rapid non-equilibrium flow. Watson & Luxmoore (1986) defined macropores as pores that conduct flux at tension less than | − 3| cm, or equivalent to a pore radius greater than 500 µm. They found that three-quarters of the flux in tension infiltrometers occurs between 0 and −3 cm tension. Another useful parameter is macroscopic capillary length (λ_c) which is equivalent to the wetting front suction in the Green and Ampt equation. It is also equivalent to the inverse of α determined by Eqn. (2) (Nachabe, 1996).

White & Sully (1987) compared the λ_c for disturbed and undisturbed soils. They found that λ_c values increased drastically in disturbed soils. Furthermore, their results showed that λ_c values of less than 80 mm were considered indicators of preferential flow. For the purpose of this study, two definitions were adopted. First, a macropore is a pore that drains at less than —−3— cm tension. Second, preferential flow is indicated by λ_c values of less than 80 mm.

Accordingly, the macropore flux fraction (Q_macro), defined as the fraction of the flux occurring between h =  − 3 cm and h = 0 cm, was determined by normalizing the difference between I(h =  − 0cm) and I(h =  − 3cm) with respect to I(h =  − 0cm), where I is the infiltration rate (mm/hr) (Watson & Luxmoore, 1986).

Soil sampling

After each measurement, a composite soil sample was collected from a hole excavated beneath the tension disk to a depth of 10 cm. The collected samples were analyzed for the OM and soil texture distribution. The OM was determined using the Walkley–Black (WB) method (Walkley & Black, 1934). The WB method measures the SOC. The OM is calculated by using a standard ratio of 1.72 between the OM and the SOC (i.e., OM = 1.72*SOC). The soil textural classes were determined using the hydrometer method (Bashour & Sayegh, 2007).

Exploratory data analysis

All data sets were tested for normality using the Shapiro–Wilk (SK) and Anderson-Darling (AD) normality tests (Table 2). Furthermore, each data set was compared to quantile–quantile plots (QQ-plots), the probability density function (PDF), and the cumulative distribution function (CDF) of the normal, lognormal, and gamma distributions (Figs. 3–5). Statistical analysis was performed by R language version 3.4.3.

The SK and AD tests returned p-values greater than 0.05 for OM, Q_macro, and sand content. For these data sets, the null hypothesis of normality cannot be rejected at p < 0.05 (Table 2). The QQ-plot (Fig. 3) showed that OM data fits the normal distribution better than it does the lognormal or gamma distributions. This was evidenced by the oscillation of QQ points of the lognormal and gamma distributions around the 1:1 reference line (defined as a straight line of intercept = 0 and slope = 1), while the QQ points of the normal distribution fell adjacent to the reference line except for a few extreme values. The CDF (Fig. 5) and QQ-plots (Fig. 3) of Q_macro and sand content did not include clear indications that distinguished the normal distribution from the lognormal and gamma distributions. However, the PDF (Fig. 4) of the normal distribution fits the peaks, the right and left tails of OM, Q_macro, and the sand content histograms better than the lognormal and gamma distributions.

The normality of the clay data was not confirmed by either SK or AD tests. Further examination of the clay QQ-plots (Fig. 3) showed that normal distribution exhibited the least oscillation around the reference line. Similar to the OM, O_pore, and sand content data, the PDF of the normal distribution reasonably fits the histogram of the clay content data.

The lognormal distribution PDF (Fig. 4) and CDF (Fig. 5) produced the best fit of the λ_c data. In fact, SK and AD (Table 2) confirmed the normal distribution of the log-transformed λ_c data. It was anticipated that the saturated hydraulic conductivity (K_s) and I(h =  − 0cm) are also lognormally distributed. However, normality for the log-transformation of both data sets was not confirmed by either SK or AS, possibly due to the tailing caused by extreme values. Nevertheless, the CDF plots (Fig. 5) of the K_s and I(h =  − 0cm) data suggest a reasonable fit of the lognormal or gamma distributions. Further inspection of the PDF plot (Fig. 4) shows that lognormal distribution properly represents the positively skewed K_s and I(h =  − 0cm) histograms.

Thus, the t-test may be conducted for the OM, sand and clay, Q_macro, and the log-transformed data of K_s and I(h =  − 0cm). However, non-parametric data analysis (Mann–Whitney U test, hereafter MW-test) was also applied to all data sets because the SK and AD test did not confirm the normality of clay content and the log-transformed hydraulic conductivity and I(h =  − 0cm).

In situ tension infiltrometer measurements

The exploratory analysis showed an average λ_c value of 102 mm and a geometric mean of 84 mm, which conforms with the average values that White & Sully (1987) suggested. The λ_c also showed a standard deviation of 75 mm and a positive skewness of 2.8 (Table 3).

Three segments can be distinguished based on the number of λ_c observations below and above the critical value (Fig. 2). The first segment starts near the vicinity of the KTP and then extends approximately 7 kilometers (km) to the west. The second segment extends 3 km farther downstream, then meanders northwest to an approximate distance of 7 km. The third segment is considered the point where the river starts moving in a westward direction. Figure 2 suggests that the first and third sections are dominated by λ_c observations below the critical value.

It is also beneficial to consider observations at the locations where the above critical values of λ_c were observed separately from the rest of the observations. Therefore, the determination of the average values of the observed parameters proceeded with a ranking of the data with respect to λ_c value, then the establishment of two groups by splitting the record at the λ_c critical value (Table 4). This was done for each segment as well as for the entire study area. Thus, the data included in the first group (a total of 35 samples) is associated with the preferential flow, henceforth referred to as PF, as indicated by λ_c of less than 80 mm, specifically, λ_c geometric means of 43.9, 64.6, 54.2, and 53.9 mm for segments 1, 2, 3, and the whole area, respectively (Tables 3 and 4). By contrast, the λ_c values for the second group (a total of 34 samples), or the matrix flow group (MF), are 2 to 3 times greater than that of PF. Furthermore, an unequal variances student’s t-test and an MW-test were performed on the groups representing the entire study area (Table 3). Both tests were not performed on the λ_c values because, by definition, the average λ_c of PF (data associated with λ_c < 80 mm) should be less than the average of MF (data associated with λ_c ≥ 80 mm). The geometric means of λ_c representing the whole study area were 53.9 and 133.71 mm for PF and MF, respectively.

Macropore flux, infiltration rate, and saturated hydraulic conductivity

Considering the data of the first segment, the Q_macro average values of 0.75 and 0.47 were calculated for PF and MF, respectively (Table 4). The indicated 50% difference was also observed in the third segment, where Q_macro values of 0.62 and 0.42 were recorded for PF and MF, respectively. However, an exception could be found in the results for the second segment; contrary to the general expectation, the Q_macro of MF was approximately equal to that of PF evident from Q_macro values of 0.65 and 0.66 representing PF and MF, respectively. Within the whole data set, the Q_macro values of 0.67 and 0.57 for PF and MF, respectively, were significantly different at p < 0.05 by t-test and p < 0.01 by WM-test (Table 3). The geometric means of I(h =  − 0cm) within the whole data set were 42.1 and 15.1 mm/hr for the PF and MF, respectively, while for the K_s they were 12.8 and 2.7 mm/hr for the PF and MF, respectively. Thus, the I(h =  − 0cm) and K_s were generally well associated with the inverse λ_c, as was also attested to by the significant differences at p < 0.01 by a t-test and an MW-test of the lognormally transformed data of I(h =  − 0cm) and K_s (Table 3). The largest differences in I(h =  − 0cm) values were observed within segment 1, which showed geometric means of 92.6 and 12.0 mm/hr for PF and MF, respectively (Table 4). Meanwhile, the lowest difference in I(h =  − 0cm) was observed in segment 2, as indicated by the similar I(h =  − 0cm) values of 22.9 and 17.3 mm/hr for PF and 2, respectively. A similar trend was observed for the K_s values. The K_s values (segment 1/segment 2) were 47.9/7.3 and 2.9/2.9 mm/hr for the PF and MF, respectively. Table 4 also shows that the differences in infiltration rates and K_s values between PF and MF were reduced at tensions of −3 and −6 cm. In segment 2, the average values of I(h =  − 6cm) and K(h =  − 6cm) were similar for both PF and MF. Pearson’s correlations between the hydraulic variables (Table 5) showed a moderate negative correlation of the Q_macro of the PF with I(h =  − 6cm) and K(h =  − 6cm), and moderate positive correlation with I(h =  − 0cm). The Q_macro of the MF was only correlated I(h =  − 0cm). The infiltration rate and hydraulic conductivities values of MF and PF were intercorrelated with each other at all supply tensions. A strong positive correlation was also indicated between α (inverse of λ_c) and I(h =  − 0cm) for the PF and the whole dataset.

Soil Texture

Sandy loam and loamy soils are the dominant soil types in the study area (Table 1). Nevertheless, clay loam soils were found in a few locations. The average clay content values for PF for all segments were above 10% (Table 4). The soils in PF within the first segment demonstrated an average clay content of 15.8% compared to the 8% clay content of the MF. Soils in MF within the third segment were also characterized by low clay content of an average of 5.6%. However, a similar trend was not observed in the second segment, attested to by clay contents of 11.8 and 16.6 for PF and MF, respectively. In fact, the average clay content of PF and MF for the entire dataset was not statistically different (Table 3)

The average sand contents for all sites were greater than 50%. Differences between PF and MF were found in the three segments, however, similar to clay faction, there was no statistical difference between PF and MF for the entire data set.

Organic matter

The results (Tables 3 and 4) showed that OM average values of PF were larger than that of MF. The OM values of the PF/ MF were 2.9/1.4, 2.3/1.9, 1.7/1.6 and 2.2/1.8 for segments 1, 2, 3 and entire study area, respectively. In addition, the student’s t-test and MW-test (Table 3) confirmed that the difference between the two data groups (whole data set) is statistically significant at p < 0.01. In general, the OM demonstrated a decreasing trend with increasing λ_c values (Fig. 6). Thus, indicating a negative association between OM and λ_c.

Soil Organic Carbon/clay fraction (SOC/clay)

The scatter plot of SOC versus clay content (Fig. 7) showed that the PF data are close to the 1:10 (SOC/clay) saturation line, whereas most of the MF data are below the saturation line except for soils with a clay content greater than 8%. This suggests that the SOC/clay ratios that were greater than 10% within MF were the result of low clay content. The SW and AD tests performed for SOC/clay for the PF and MF data sets returned p < 0.01, indicating that the normality assumption for the SOC/clay data cannot be accepted. Therefore, the MW-test was used to test the significance of the SOC/clay of the PF and MF data set. The average SOC/clay ratio of the PF (14%) was larger than that of the MF (13.3%). The MW-test produced a p-value of 0.047, which is close to the critical p-value level of 0.05 (Fig. 8). Further investigation of the SOC/clay ratio was carried out through the exclusion of soils with clay content below 8%. The averages of the newly amended SOC/clay ratio were 9.8 and 7.5% for the PF and MF, respectively. The MW-test indicated that the differences between the amended groups were statistically significant at p < 0.05 (p = 0.021).

Macropore flux

Macropores are major pathways for preferential flow and transport. Thus, high Q_macro values are expected at sites that are susceptible to preferential flow. However, this is anticipated due to the empirical correlation between the flow capacities at the supply tension ranges of 0 to −3 cm and −3 to −6 cm, measured by Q_macro and λ_c, respectively (e.g., Jarvis et al., 2013). The significant difference between the Q _macro values of PF and MF reflects the presence or absence of preferential flow as deduced from λ_c values (Table 3). The test of significance also indicated that high I(h =  − 0cm) values were associated with PF, which can be inferred from the considerable difference between I(h =  − 0cm) of PF and MF within segments 1 and 3. However, the influence of preferential flow on I(h =  − 0cm) was not clearly observed in segment 2, where the observed average I(h =  − 0cm) values were 27.0 and 25.8 mm/hr for PF and 2, respectively. These values are considered to be the upper end of the I(h =  − 0cm) range reported by Brouwer et al. (1988). Thus, some aspects of the preferential flow and transport may exist in soils located within the second segment, even when the preferential flow process is not fully apparent.

Also, the relatively low K_s values of segment 2 (as compared to those of segments 1 and 3) can facilitate preferential flow due to the build-up of saturation in the soil as a result of the slow movement of soil water, which subsequently induces macro-pore flow (Koestel et al., 2013). Furthermore, Larsbo, Koestel & Jarvis (2014) showed that macropore flow is not necessarily related to large macroporosities. Preferential flow can occur in soils of low macroporosity through the few preserved macropores. This explains the proximity of I(h =  − 0cm) for PF and MF within segment 2 despite the disparity in λ_c values (A.M. of 65.6 and 136.1 for PF and MF, respectively). In fact, the reduction of macroporosity increases the severity of preferential flow because the elimination of some macropores and the decrease of near-saturated hydraulic conductivity can quickly bring the soil water content to saturation, subsequently activating the remaining macropores that at saturation can conduct a larger amount of water than the soil matrix, thus dominating the flow in the soil profile. In this study, the K(h =  − 6cm) and I(h =  − 6cm) values were negatively correlated with the Q_macro values of the PF data set. Thus, high K(h =  − 6cm) can minimize or prevent the activation of soil macropores by keeping the soil water content near saturation. However, high irrigation application rates that exceed the intake rate of the soil matrix or intense rainfall events will probably lead to soil saturation and preferential macropore flow. Rahbeh, Srinivasan & Mohtar (2019) indicated that 20 mm of daily rainfall was enough to trigger preferential flow in a small irrigated site near the KTP.

Role of organic matter and soil texture

Soil analysis of OM and soil texture are independent of the in situ tension infiltrometer measurements. Therefore, the results of the soil textual analysis and OM were examined for possible association with the preferential flow. The difference between the average clay contents of PF and MF was insignificant (Table 3). Therefore, a direct correlation between clay content and preferential flow cannot be inferred from the results. However, it is worth noting that, in this study, preferential flow was associated with clay content of between 11–16%. This is in agreement with Iversen et al. (2012) and Ghafoor et al. (2013), who reported a positive correlation between preferential transport and fine-textured soils. Also, it aligns with the threshold of 8% clay for macropore formation as suggested by Koestel, Moeys & Jarvis (2012).

However, the positive association of clay content with preferential flow was not observed in the second segment, attested to by clay content of 11.8 and 16.6 for PF and MF, respectively. Apparently, low clay content may impede preferential flow. However, high clay content does not necessarily create favorable conditions for preferential flow (Koestel, Moeys & Jarvis, 2012).

Preferential flow is also influenced by OM (Dexter et al., 2008). The role of OM as a cementing agent is well known (e.g., Hillel, 2004). OM improves the structure and physical properties of soils. It facilitates the formation of soil macropores, thus increasing the infiltrability of the soil as well as improving the aeration and redistribution of soil water. These effects are favorable from an irrigation point of view, but may also increase the potential of preferential flow, thus leaving the rechargeable groundwater aquifer susceptible to environmental risks. The positive association between OM and preferential flow is indicated by the significant difference in OM between PF and MF (Table 3), and also by the decreasing OM trend with increasing λ_c values (Fig. 6).

The differences between the SOC/clay ratios of soils associated with λ_c < 80 mm and λ_c ≥ 80 mm were marginally significant. Therefore, SOC/clay may not be a sufficient indicator of susceptibility to preferential flow (Fig. 8). The lack of statistical significance may be attributed to high SOC/clay ratios resulting from clay content below the 8% minimum threshold required for soil aggregation and the formation of macropores (Koestel, Moeys & Jarvis, 2012), thus negating the influence of OM. The exclusion of soils with clay content below 8% increased the statistical difference in SOC/clay ratio between PF and MF (Fig. 8). Thus, the combination of an SOC/clay ratio close to or greater than 10% and clay content above 8% may be used as a criterion for indicating susceptibility to preferential flow.

Apparently, for the types of soils along the ZR, the macropore formation favors fully complexed clay content or soils with an SOC/clay ratio slightly greater than 10%, provided that the soil has a clay content of 8%, which is the minimum required for aggregation. These soils are generally more susceptible to preferential flow than are soils of low clay content. Therefore, the observation that the preferential flow was influenced interdependently by OM and clay content aligns with the results obtained by Dexter et al. (2008), Zaffar & Lu (2015), and Larsbo et al. (2016).

It should be noted that previous studies indicated a positive correlation between preferential flow and clay content, whereas OM was either secondary in importance to clay content or essential to the stability of the soil structure but did not directly correlate with macropore formation (Karup et al., 2016). Furthermore, Ghafoor et al. (2013) and Larsbo et al. (2016) observed that high OM reduced preferential transport in soils that were rich in clay. The lack of a correlation between OM and preferential flow may also be attributed to the lack of variation in OM or SOC/clay (Koestel et al., 2013; Karup et al., 2016). However, the more important role of OM in this study, as compared to its role in previous studies, can also be attributed to climate differences (Jarvis et al., 2013). Most studies of preferential flow and transport were conducted in temperate to humid climates with a mean annual precipitation range between 300 and 1,100 mm (Jorda et al., 2015). By contrast, the ZR is an arid region characterized by an annual mean rainfall of less than 150 mm. Jorda et al. (2015) indicated that the temperature and precipitation norms are key controls of near-saturated hydraulic conductivity. Furthermore, a regional analysis using soil data compiled from across the USA (Hirmas et al., 2018) showed greater macropore formation in dry regions than in wet regions. This was indicated by the inverse relationship between macroporosity and precipitation variables that included an average magnitude of single rainfall events and mean annual precipitation. They inferred that the decline of macroporosity in humid regions may be attributed to intense rainfall events that caused soil macropores to shift into smaller pores.

Climatic conditions may also influence local levels of OM and clay content. Clay content is proportional to annual rainfall. Therefore, soils in temperate to humid regions are usually richer in clay than are soils in arid regions (Jenny, 1941). Furthermore, high clay content enhances the storage of SOC (Hartati & Sudarmadji, 2016). In this study, clay content and OM were less than 20 and 2%, respectively, in most locations. These values are generally lower than the range of values reported from temperate regions. For example, Ghafoor et al. (2013) reported a range of clay content between 5 and 40% and an OM range between 1.5 and 4%. Also, Larsbo et al. (2016) collected undisturbed soil columns from a site characterized by high clay content (64–70%) and a variable OM range of between 4 and 15%. Larsbo et al. (2016) explained that high OM induced the increased formation of smaller macropores with a diameter of 200–600 µm, which resulted in the reduction of preferential transport. Nemes, Rawls & Pachepsky (2005) showed that high clay content drives the pedotransfer functions into an inverse correlation between OM and K_s; however, for low OM and clay content, K_s is positively correlated with OM. Rawls, Nemes & Pachepsky (2004) suggested that K_s increases with an increase of OM in clay content of less than 30%. At a higher clay content, the relation between OM and K_s becomes more complex. Therefore, it is plausible that at low OM and clay content (e.g., 10–20% clay content and 0.5–2% SOC), the preferential flow was positively associated with OM due to the formation of macropores when SOC/clay approached or exceeded the 1:10 saturation level.

However, macropores are also a reflection of the activity of soil fauna and plant roots (Beven & Germann, 1982). Soil OM can be an outcome rather than the cause of preferential flow. In fact, the preferential flow pathway can alter the distribution of OM (Hagedorn et al., 2015). Cheng et al. (2014) observed higher OM in the preferential flow pathways than in the soil matrix, which they attributed to the presence of decaying roots in the macropores. The ZR valley is intensely cultivated. Farmers run up to three cultivation cycles per year. The intense irrigated agriculture activity increases the lateral and vertical distribution of roots, and also enhances the growth of worms and insects. The steady biological activity in the soil subsurface generates new pathways for express water and contaminant movement and, simultaneously, the accumulation of new organic matter in the soil.