Investigation of the Flux–Concentration Relation for Horizontal Flow in Soils

Abstract

The objective of the present work is to investigate the flux–concentration (F(Θ)) relation, where Θ is the normalized soil volumetric water content for the case of one-dimensional horizontal flow, subject to constant concentration conditions. More specifically, the possibility of describing F(Θ) by an equation of the form F(Θ) = 1 − (1 − Θ)^p+1 is examined. Parameter p is estimated from curve-fitting of the experimentally obtained λ(Θ) data to an analytic expression of the form (1 − Θ)^p where λ is the well-known Boltzmann transformation λ = xt^−0.5 (x = distance, t = time). The results show that the equation of (1 − Θ)^p form can satisfactorily describe the λ(Θ) relation for the four porous media tested. The proposed F(Θ) function was compared with the limiting F(Θ) function for linear and Green–Ampt soils and to the actual F(Θ) function. From the results, it was shown that the proposed F(Θ) function gave reasonably accurate results in all cases. Moreover, the analytical expression of the soil water diffusivity (D(Θ)) function, as it was obtained by using the equation for λ(Θ) of the form (1 − Θ)^p, appears to be very close to the experimental D(Θ) data (root mean square error (RMSE) = 0.593 m²min⁻¹).

1. Introduction

For a rational and productive irrigation water application, one needs to know well how water infiltrates into the soil, and how the soil moisture profiles evolve. This can be achieved by solving the appropriate Richards’ equation, which is feasible only when major soil hydraulic properties are known, together with the initial and boundary conditions imposed in each case.

To solve the unsaturated flow equation, many procedures have been proposed, e.g., analytical, semi analytical, finite differences, and finite elements. Philip [1] introduced the flux–concentration (F(Θ)) relation in order to develop a quasi-analytical technique for solving the unsaturated flow equation. The technique is rather simple and gives a holistic view of the phenomenon [1,2].

F(Θ) expresses the relation of soil water flux–density with the normalized volumetric soil water content, and can be used in an attempt to elucidate further the flow phenomena in unsaturated soils, both in the horizontal (absorption, as well as desorption) and in the vertical cases. The form of the F for the case of horizontal absorption and the normalized volumetric soil water content is given by the following equations:

(1)$F=\frac{q}{q_0},\mathsf{\Theta }=\frac{\theta -{\theta }_{in}}{{\theta }_0-{\theta }_{in}}$

where θ_in is the initial volumetric soil water content of the soil column, considered uniform throughout the column; θ₀ is the volumetric soil water content at the soil surface (x = 0); q is the soil water flux density at a position x where soil water content is θ; and q₀ is the soil water flux density at x = 0, the soil surface, where water is absorbed under constant concentration conditions with θ = θ₀ = θ_s (θ_s = volumetric soil moisture at saturation) when the soil water pressure head is H = H₀ = 0 cm column of water.

The values of F and Θ lie in the range of 0 to 1. A zero value for F and Θ when q = 0 is far enough from the soil infiltration surface, where θ = θ_in; the value is 1 at the soil surface, where q = q₀ and θ = θ₀ [1]. In general, the F is a function of θ, θ_in, θ₀, and time (t) and it is, in most cases, unknown a priori [3].

For the prediction of the flow regime in porous media, using the flux–concentration relationship, either an iterative procedure should be applied, or accurate estimates of F(Θ) where feasible are used [2].

A number of researchers have proposed analytical expressions for the F(Θ) relationship to solve the flow equation for the case of the horizontal absorption under constant pressure head conditions, applied at the soil surface (x = 0).

Smiles et al. [4] proposed an empirical expression to describe the measured F(Θ) of a fine sand soil:

(2)$F\left(\mathsf{\Theta }\right)=1-{\left(1-\mathsf{\Theta }\right)}^{1.19}.$

Vauclin and Haverkampt [5] proposed the following simple expression

(3)$F\left(\mathsf{\Theta }\right)=\frac{2\mathsf{\Theta }}{\mathsf{\Theta }+1},$

while Evangelides et al. [6] proposed the empirical expression

(4)$F\left(\mathsf{\Theta }\right)=2\mathsf{\Theta }-\left[1- {\left(1-\mathsf{\Theta }\right)}^m\right]$

with the value of parameter m lying in the range 0.718–0.867, and the value m = 0.8 considered as the most appropriate, which is very close to the mean value.

Recently, Ma et al. [7] proposed the expression

(5)$F\left(\mathsf{\Theta }\right)=\mathsf{\Theta }\left[{1+\mathrm{a}}_2\left({1-\mathsf{\Theta }}^{\frac{1}{a_2}}\right)\right]$

with the parameter a₂ being a function of the initial soil water saturation and the soil pore structure index that reflects the shape of the soil water retention curve.

From the abovementioned expressions, those where F is a function of Θ with a parameter for soil properties (Equations (4) and (5)), compared to those where F is a function of Θ without any parameters (Equations (2) and (3)), are more flexible and can simulate F(Θ) relationships with higher accuracy for a wide range of texture soils [7].

Philip [1] investigated the behavior of F for two distinct cases of soils, characterized by their D(θ) functions—i.e., soils, where their D(θ) function is a delta function (the so called Green–Ampt soils, where the moisture profiles advance step-wise), as well as soils whose D(θ) function is constant, independent of moisture content θ (the so called linear soils). In real soils, the F(Θ) relationship will lie in between the Green–Ampt and linear soils [1].

Clothier et al. [8] investigated the possibility of estimating the expression $\lambda$(Θ), where $\lambda$ is the Boltzmann transformation $\lambda \left(\theta \right)=x{t}^{-0.5}$ for each soil from the experimental data of a horizontal absorption experiment, according to the equation below (Philip [9], Table 1, no. 2):

(6)$\lambda \left(\mathsf{\Theta }\right){=\mathsf{\lambda }}_i{\left(1-\mathsf{\Theta }\right)}^p,$

where λ_i is the maximum value of λ(Θ) (i.e., when Θ = 0). Parameter λ_i considered to be a characteristic measure of the wetted region in a horizontal absorption experiment, and is influenced only by θ_in for every soil [10,11,12]; p is a fitting parameter.

Equation (6) has the advantage of an easy estimation of parameters λ_i and p; by a log-transformation, Equation (6) becomes linear, with the slope equal to p and transect equal to logλi when the λ(Θ) is described by Equation (6). This procedure appears preferable to least squares fitting.

It is to be mentioned here that previous similar studies [6] could obtain the F(Θ) relationship using more complex expressions for it. In the present study, we use a simple, mono-parametric F(Θ) expression using a mono-parametric λ(Θ) expression, which adequately describes the experimental data [8].

In this context, the main objectives of this study are (1) tο test the possibility of the mono-parametric λ(Θ) expression to describe the experimental data, as the equation $\lambda \left(\mathsf{\Theta }\right){=\mathsf{\lambda }}_i{\left(1-\mathsf{\Theta }\right)}^p$ through the proper selection of a fitting parameter p; and (2) tο investigate whether the experimental F(Θ) relations are described accurately enough from the analytical, mono-parametric function $F\left(\mathsf{\Theta }\right)=1-{\left(1-\mathsf{\Theta }\right)}^{p+1}$ for a wide range of soil types, as well as finding out the physical meaning of the parameter p and the range of its values.

2. Theory

Applying Darcy’s law and the mass conservation principle for the one-dimensional horizontal infiltration case, one easily gets the partial differential equation

(7)$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial x}\left[D\left(\mathsf{\theta }\right)\frac{\partial \theta }{\partial x}\right],$

where x is the horizontal axis, t the time, and D(θ) the soil water diffusivity. The solution of Equation (7) under the following initial and boundary conditions:

• (a). t = 0, x > 0, θ = θ_in

• (b). t > 0, x = 0, θ = θ₀

• (c). t > 0, x → ∞, θ = θ_in

could be obtained in terms of the Boltzmann transformation $\lambda \left(\theta \right)=x{t}^{-0.5}$, considered a function of θ only as

(8)${{\displaystyle \int }}_{{\theta }_{in}}^{\theta }\lambda \left(\theta \right)d\theta =-2D\left(\theta \right)\frac{d\theta }{d\lambda },$

following the initial and boundary conditions θ = θ_in, λ → ∞, and θ = θ₀, λ = 0.

Introducing the Boltzmann transformation function $\lambda \left(\theta \right)=x{t}^{-0.5}$ in Darcy’s law, one gets

(9)$q=-{\rm K}\left(\theta \right)\frac{dH}{dx}=-D\left(\theta \right)\frac{d\theta }{d\lambda }\frac{\partial \lambda }{\partial x}=-t^{-0.5}D\left(\theta \right)\frac{d\theta }{d\lambda },$

and the combination of Equations (1), (8), and (9) gives the following expression for F

(10)$F=\frac{{{\displaystyle \int }}_{{\theta }_{in}}^{\theta }\lambda \left(\theta \right)d\theta }{{{\displaystyle \int }}_{{\theta }_{in}}^{{\theta }_0}\lambda \left(\theta \right)d\theta },$

From the abovementioned equation, it is obvious that function F, for the case of the horizontal absorption, is independent of time and depends on the values of the volumetric soil water, θ, θ_in and θ₀.

When the D(θ) function is a delta function, the flux concentration relationship F for the horizontal absorption is given by

(11)$F=\mathsf{\Theta },$

while for the case where the D(θ) function is constant, the flux concentration relationship F is given by

(12) $F=\exp \left[-{\left(\mathrm{inverfc}\left(\mathsf{\Theta }\right)\right)}^2\right],$

White et al. [13] have shown that Equation (12) could be reasonably approximated by the expression

(13) $F=\sin \left[\left(\frac{\mathsf{\pi }}{2}\right){\mathsf{\Theta }}^{\frac{\mathsf{\pi }}{4}}\right].$

It is obvious from Equations (8) and (10) that the F(Θ) relationship depends on the form of D(θ).

From Equation (6), the soil sorptivity S_c function could be estimated by analytical expression [8], since

(14) $S_c=\left({\theta }_0 - {\mathsf{\theta }}_{in}\right){{\displaystyle \int }}_0^1{\lambda }_i{\left(1-\mathsf{\Theta }\right)}^p\mathrm{d}\mathsf{\Theta }=\frac{\left({\theta }_0 - {\mathsf{\theta }}_{in}\right){\lambda }_i}{p+1},$

From the combination of Equations (6), (9), and (10), it could be shown that the F(Θ) relationship could be estimated analytically, according to the mono-parametrical analytical expression

(15)$F\left(\mathsf{\Theta }\right)=1-{(1-\mathsf{\Theta })}^{p+1}.$

A similar expression for F(Θ) was presented by Smiles et al. [4] in a horizontal absorption experiment in sand, where the parameter p had a constant value p = 0.19.

Applying the Bruce and Klute [14] method, the soil water diffusivity function D(Θ) can be estimated according to the expression

(16) $D\left(\mathsf{\Theta }\right)=-\frac{{{\displaystyle \int }}_0^{\mathsf{\Theta }}\lambda \left(\mathsf{\Theta }\right)d\mathsf{\Theta }}{2\frac{d\mathsf{\Theta }}{d\mathsf{\lambda }}}.$

In this respect, it is rather easy, by introducing the λ(Θ) expression given in Equation (6), to get an analytical expression for D(Θ) ([8]), as

(17) $D\left(\mathsf{\Theta }\right)=\frac{p\left(p+1\right)S^2\left[{\left(1-\mathsf{\Theta }\right)}^{p-1}-{\left(1-\mathsf{\Theta }\right)}^{2p}\right]}{2{\left({\theta }_0{-\mathsf{\theta }}_{in}\right)}^2}.$

In practice, this means that for soils with a λ(Θ) expression, such as the one given by Equation (6), their diffusivity function D(Θ) should be described analytically by Equation (17), and the flux concentration relationship F(Θ) by Equation (15).

3. Materials and Methods

Horizontal infiltration experiments were conducted (Scheme 1) in three porous media with different soil textures: a sandy loam (SL) (13.2% clay, 8% silt, 78.8% sand), a loam (L) (20% clay, 38% silt, 42% sand), and a silty clay loam (SiCL) (36.5% clay, 52% silt, 11.5% sand). The respective soils were air-dried, ground, and passed through a 2 mm sieve. The soil samples were uniformly packed in a rectangular column 5 cm in width, 2 cm in height, and 50 cm long. The upper soil surface was open to the air, thus securing that the soil air pressure at the soil pores was at atmospheric pressure (Scheme 1). The experimental data from horizontal absorption experiments conducted in sand (S) by Poulovassilis [15] were also used.

The volumetric soil water θ_j(x_j, t*) profile, after a time period t* was elapsed from the beginning of the experiment, was determined by cutting, as quickly as possible, the soil column into small rectangular pieces. Specifically, the soil column was sectioned in 0.01 m increments, and the volumetric water content of each rectangular piece θ_j was determined by using the gravimetric water content and the dry soil bulk density (ρ_φ). The θ_j value corresponds to the center of the soil samples, at distance x_j from the soil infiltration surface.

The initial values of θ and θ_in were those associated to the air-dry soil, and the condition of the soil infiltration surface was that of a constant pressure head H value (H = 0 at x = 0, with x = 0 denoting the soil infiltration surface). The zero value of H was maintained by a Mariotte device (Scheme 1). A thin wire mesh at x = 0 was installed to keep the soil at rest, and also provide the least possible resistance to the soil water entry.

The time duration of these experiments varied according to the soil type—i.e., smaller values for coarse-textured soils and larger values for fine-textured soils. For the sand (S), it was 10.4 min [15], for the SL it was 194 min, for the L it was 251 min, and for the SiCL it was 1630 min. During the experimental process, a continuous monitoring of the water volume entering the column was obtained, thus allowing the cumulative infiltration I(t) experimental curve to be determined. The soil sorptivity S is immediately available from the well-known relationship [16]

(18) $S=\frac{I}{t^{0.5}},$

4. Results and Discussion

In Table 1, some soil characteristics, such as θ₀, θ_in, ρ_φ, λ_i and S (Equation (18)), together with fitting parameter p and the value of S_c (Equation (14)), are shown for each porous media used in this study. It is easily noticeable that the values of S and λ_i depend strongly on the soil type, and they tend to decrease as soils become finer in texture. The same trend for the values S and λ_i for six different soils were presented from McBride and Horton [10].

In Figure 1, the λ(θ) profiles are presented, as these were obtained from the experimental θ(x,t) data and the well-known Boltzmann transformation ($\lambda \left(\theta \right)=x{t}^{-0.5}$), together with the fitted λ(θ) relationship (Equation (6)), after the fitting parameter p was properly selected. From the comparison, it is shown that Equation (6) is suitable for the description of the λ(θ) relationships for the horizontal absorption experiments. In each soil, the p-value was determined as the one where the root mean square error (RMSE) from a series of neighboring p-values was the least.

One could argue, considering that the λ(Θ) relationship is unique for each soil, independent of the duration of the experiment, that the two parameters of Equation (6) (λ_i and p) do not represent simple fitting parameters, but are related to the soil’s hydraulic properties. Some researchers insist that parameter λ_i may be considered as an index of the soil’s hydraulic properties [10,11]. Moreover, Shao and Horton [11] correlated parameter p with the parameter n (p = 1/n) of the equation of van Genuchten [17], which describes the soil moisture retention curve.

In what follows, an investigation of the characteristics of parameter p will be carried out. According to Equation (14), the ratio S/λ_i(θ₀ – θ_in) is equal to (p + 1)⁻¹. For the Green–Ampt soils (D is a Dirac delta-function of θ), the ratio S/λ_i(θ₀ – θ_in) should be unit, thus p would be zero. Similarly, for the linear soils, the expression (p + 1)⁻¹ would be 0.31, and the value of p will be 2.23 [8]. Consequently, the values of p are related to the form of the diffusivity D(θ) function. In order to examine this relationship, Equation (6) is rewritten as in Equation (19):

(19)$\frac{\lambda }{{\lambda }_i}={(1-\mathsf{\Theta })}^p\to \mathsf{\Theta }=1-{\left(\frac{\lambda }{{\lambda }_i}\right)}^{\frac{1}{p}}.$

In Figure 2, the Θ(λ/λ_i) relationship is shown for various values of the parameter p (0 < p < 2.23). From Figure 2, it is shown that all the Θ(λ/λ_i) relations lie in the area with its borders defined by the values of the parameter p (i.e., p $\to 0$; D(θ) Dirac delta function) and p = 2.23 (D(θ) constant). From the investigation of the λ(θ) relations in this study, it is found that p-values fall in the range 0.1 < p < 0.4. Also, Evangelides et al. [6] showed that p-values obtained using data from horizontal infiltration experiments in seven soils fall in the range 0.149 < p < 0.389. In other words, the range of p is narrower than the range $0\le p\le 2.23$.

After the parameter p selection, using also the λ_i, θ₀, and θ_in values for each soil, a re-estimation of sorptivity S_c was performed using Equation (14). From the values of S and S_c presented in Table 1, S_c values are reasonably close to the experimental values of S (S = I/t^0.5) for three out of the four soils examined (absolute values of RE: 1.55% < RE < 3.08%). The relative error value for the SiCL soil appears to be rather high (13.27%). This could be attributed to the long time duration of the experiment (1630 min) and unavoidable soil water evaporation from the upper soil surface. In any case, the overall differences are small, and therefore one may consider that Equation (14) can lead to a quick and reliable way of estimating S from a set of horizontal absorption experimental data.

In Figure 3, the D(θ) data for the loamy soil (L) are shown. Closed circles denote data obtained by the Bruce and Klute [14] method (Equation (16)), while open circles are data obtained from Equation (17). The comparison of the above indicates that Equation (17) could reasonably describe the experimental D(θ) data. One should note that by using the analytical expression (Equation (17)), the problem of differentiating the experimental data $\frac{d\mathsf{\Theta }}{d\lambda },$ in which there is scatter, is overcome. Moreover, the problem of estimating the slope $\frac{d\mathsf{\Theta }}{d\lambda }$ near saturation, where the λ(Θ) relationship is almost parallel to the λ-axis, is overcome by the application of the analytical expression (Equation (17)). Similar results were also obtained for the other soils under present investigation.

In Figure 4, the F(Θ) relationship for the linear Equation (13) and the Green–Ampt soils (Equation (11)) is shown, together with that obtained according to Equation (15) for all soils examined. The experimental F(Θ) points for all soils were obtained by the application of Equation (10) and the measured value of the sorptivity S (Equation (18)). From the results, it is shown that the one-parameter Equation (15) gave practically the same values for the F(Θ) relationships as the ones obtained experimentally, and lies between the two limits (linear and Green–Ampt soils) in all cases that were examined.

It is worth investigating the ability of the equation F(Θ) = 1 − (1 – Θ)^p+1 to describe the upper and lower theoretical boundaries of F(Θ). As has already mentioned, when p = 0 the value of F = Θ—i.e., it converges to the lower limit (D(θ) function is a delta function). However, the values of F for p = 2.23 (resulted from (p + 1)⁻¹ = 0.31 [8]) differs from the values of F resulting from Equation (13), which is the upper theoretical limit of F(Θ). Specifically, it gives higher values of F(Θ) than the theoretical curve. The fitting of the F(Θ) = 1 − (1 – Θ)^p+1 equation to the theoretical curve (Equation (13)) showed that the equation F(Θ) = 1 − (1 – Θ)^2.36 gives very good results of F(Θ) estimation over the entire range of Θ [7]. Therefore, the variation range of the parameter p is between 0 and 1.36, and the resulting shape parameter value for the linear soils ((p + 1)⁻¹ = 0.423) is greater than that presented by Clothier et al. [8].

It can be concluded that F(Θ) = 1 − (1 – Θ)^p+1 seems to be appropriate functional form for describing actual flux–saturation curves of general soils, and its parameter p has a physical meaning, i.e., it represents the shape of the soil moisture profile.

Furthermore, one could observe that from the study of Evangelides et al. (Table 2 [6]), where seven different soils are studied, the fitting parameters m (Equation (4) [6]) and n = p + 1 (Equation (15)) are strongly linearly related, as $n=p+1=-1.483m+2.4375$. If one uses m = 4/5, as Evangelides et al. [6] proposed, then this leads to $\overline{n}=1.251$ which corresponds to our p value 0.25. In this respect, Equation (4) [6], with m = 0.8, becomes equivalent to Equation (15) in the present study, with p = 0.25.

5. Conclusions

The mono-parametric equation of the form (1 – Θ)^p can reliably describe the λ(Θ) relationship after the proper selection of the parameter p for a relatively large range of soils. It is also shown that the analytical expressions of the soil hydraulic diffusivity D(θ) and soil sorptivity S approach the experimental ones well. Moreover, for the case of the D(θ) relationship, there is the advantage of obtaining values near saturation, where the classical methodology of Bruce and Klute might be inadequate.

During the laboratory experiment, a new mono-parametric analytical function was used for F(Θ) evaluation, where the parameter p is estimated from curve-fitting of the experimentally obtained λ(Θ) data to an analytical expression of the form (1 – Θ)^p. Parameter p seems to be strongly related to soil hydraulic properties, and further investigation is needed to find this exact relationship. The analytical F(Θ) relationship, for all soil types investigated, approaches the experimental ones very well, and lies within the limiting F(Θ) function for linear and delta function soils. In addition, the upper and lower limit curves of F(Θ) calculated by the proposed expression were consistent with theoretical curves.

Author Contributions

All authors contributed equally to this work.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Figures, Scheme and Table

Scheme 1. Experimental device.

Figure 1. The λ(θ) relationships, as obtained directly from experimental θ(x,t) data for each soil, and the estimated values according to Equation (6).

Figure 2. The relationship Θ(λ/λi) for different values of the parameter p.

Figure 3. A comparative presentation of the relationship D(θ) as obtained according to the Bruce and Klute [14] method (Equation (16)) and Equation (17) for the loamy (L) soil.

Figure 4. Comparative presentation of the relationship F(Θ) for the linear soils and Green–Ampt soils, as well as the F(Θ) calculated according to Equation (15) and the experimental F(Θ) values described in the text for the four soils examined.

Figure 4. Comparative presentation of the relationship F(Θ) for the linear soils and Green–Ampt soils, as well as the F(Θ) calculated according to Equation (15) and the experimental F(Θ) values described in the text for the four soils examined.

Table 1

The values of the soil characteristics θ₀, θ_in, ρ_φ, and S, together with the values of the parameters λ_i, p, and S_c for each soil examined. The value of the parameter p was estimated from Equation (6). The value of sorptivity S was obtained directly from the experimental I(t) data (Equation (18)), while S_c comes from Equation (14). RE denotes the absolute relative error between actual S and estimated S_c values for each soil examined.

Porous Media	θ0 (m3 m−3)	θin (m3 m−3)	ρφ (g cm−3) 1	λi (cm min−0.5)	p	S (cm min−0.5)	Sc (cm min−0.5)	RE (%)
Sand (S)	0.284	0.06	-	12.210	0.35	1.930	1.980	2.59
Sandy loam (SL)	0.418	0.015	1.41	3.510	0.12	1.300	1.26	3.08
Loam (L)	0.465	0.022	1.12	1.545	0.1	0.644	0.634	1.55
Silty clay loam (SiCL)	0.511	0.028	1.22	0.656	0.18	0.309	0.268	13.27

¹ Dry soil bulk density value for sand soil is not referred to by Poulovassilis [15].