Scaling to generalize a single solution of Richards' equation for soil water redistribution

SOILS AND PLANT NUTRITION

Scaling to generalize a single solution of Richards' equation for soil water redistribution

Morteza SadeghiI, ^*; Bijan Ghahraman^I; Kamran Davary^I; Seyed Majid Hasheminia^I; Klaus Reichardt^II

^IFerdowsi University of Mashhad/College of Agriculture - Dept. of Water Engineering - 91779-4897 - Mashhad, Iran

^IIUSP/CENA - Lab. de Física do Solo, C.P. 96 - 13418-900 - Piracicaba, SP - Brasil

ABSTRACT

Using scaling methods, a single solution of Richards' equation (RE) will suffice for numerous specific cases of water flow in unsaturated soils. In this study, a new method is developed to scale RE for the soil water redistribution process. Two similarity conditions are required: similarity in the shape of the soil water content profiles as well as of the water flux density curves. An advantage of this method is that it is not restricted to a specific soil hydraulic model - hence, all such models can be applied to RE. To evaluate the proposed method, various soil textures and initial conditions were considered. After the RE was solved numerically using the HYDRUS-1D model, the solutions were scaled. The scaled soil water content profiles were nearly invariant for medium- and fine-textured soils when the soil profile was not deeply wetted. The textural range of the soils in which the similarity conditions are held decreases as the initial conditions deal with a deeply wetted profile. Thus, the scaling performance was poor in such a condition. This limitation was more pronounced in the coarse-textured soils. Based on the scaling method, a procedure is suggested by which the solution of RE for a specific case can be used to approximate solutions for many other cases. Such a procedure reduces complicated numerical calculations and provides additional opportunities for solving the highly nonlinear RE as in the case of unsaturated water flow in soils.

Keywords: Nielsen-similarity, spatial variability, invariant solutions

Introduction

Scaling methods based on the "similar media" concept (Miller and Miller, 1956) were developed to cope with the spatial variability of soils (Warrick et al., 1977; Sharma et al., 1980; Ahuja and Williams, 1991; Kosugi and Hopmans, 1998; Tuli et al., 2001; Kozak and Ahuja, 2005; Roth, 2008; Sadeghi et al., 2010). Vereecken et al. (2007) comprehensively reviewed the scaling methods developed during the past years.

Scaling has proven its success also as a tool for numerical analyses. Using scaling methods, a single solution of Richards' equation (RE) will suffice for numerous specific cases of unsaturated water flow. Hence, these methods considerably reduce the calculations required for heterogeneous soils (Warrick and Hussen, 1993). So far, various methods for scaling RE have been proposed (Reichardt et al., 1972; Warrick and Amoozegar-Fard, 1979; Warrick et al., 1985; Vogel et al., 1991; Kutilek et al., 1991; Warrick and Hussen, 1993; Neuweiler and Cirpka, 2005). Using specific scaling factors, these methods allow a linear transformation of RE variables to achieve invariant solutions for a set of similar soils. This similarity may be defined based on microscopic-scale geometry (Miller and Miller 1956), shape of soil hydraulic functions (Simmons et al., 1979), or a linear variability concept (Vogel et al., 1991).

Scaling methods can be divided into two classes (Kutilek and Nielsen, 1994). In the first class, the scaling factors are derived to unify the soil hydraulic functions into a single curve and therefore, are invariant quantities for each soil. Using these methods, scaled RE will be invariant for similar soils provided that the scaled boundary and initial conditions are the same. This provision limits the applicability of such methods developed, for example, by Warrick et al. (1985), and Vogel et al. (1991). In the second class, the scaling factors are defined by considering the imposed boundary and/or initial conditions. A clear advantage of these methods over the first class is that the scaled RE is invariant to the boundary and/or initial conditions. However, these methods are limited to specific hydrological processes and/or soil hydraulic models. For example, the method of Warrick and Hussen (1993), developed for infiltration and redistribution, applies only for the Brooks-Corey hydraulic functions.

The objective of this study, following the second class, is to scale RE focused on the redistribution process. Scaling of RE solutions is proposed instead of solving the scaled RE. By this procedure, this new method is not restricted to a specific hydraulic model - hence, all existing and future models can be used for the redistribution solutions.

Theory

Consider a one-dimensional Richards' equation (RE) of the form:

where θ [L³ L^-3] is the volumetric soil water content, K [L T^-1] the unsaturated hydraulic conductivity, h [L] the soil water matric potential head (i.e. absolute value of the soil water pressure head), z the vertical position coordinate below soil surface, and t the time.

Considering a soil water redistribution process without surface evaporation, the following conditions are imposed on the upper (z = 0) and lower boundary (z = L) of the solution domain dealing with zero water flux density and free drainage, respectively:

where q [L T^-1] is the water flux density. The following initial conditions are considered:

where θ_fi and z_fi are the initial values of the wetting front water content θ_f(t) and depth z_f(t), respectively; θ_ι is the initial soil water content below the wetting front (WF). Figure 1 graphically describes the imposed boundary and initial conditions on the solution domain.

A θ-based solution of Eq. (1), considering conditions (2) to (4), yields soil water content profiles (SWCP), θ(z), during redistribution. Considering the shape of the SWCPs during redistribution, two general forms can be found: the SWCP (a) with the θ gradient everywhere positive and with the θ distribution almost uniform above a sharp WF, and (b) with a further wetting below advancing as a step-like profile while the soil desaturates near the surface (Youngs, 1990).

Scaling Method

Following relationships are proposed to scale θ and z:

by which the scaled values of θ_fi and z_fi will become unity and the scaled initial conditions will be invariant as follows:

In Eqs. (5) to (7) θ^* and z^* are the scaled soil water content and depth, respectively. It is assumed that t can be scaled using a constant scaling factor Τ (i.e. t^*=Τt, where t^* is the scaled time) so that the scaled solutions of RE, θ^*(z^* ,t^*), be invariant for a set of cases (i.e. specific soils and initial conditions). To do so, a primary requirement is that the shape of the SWCPs be similar for all the cases. Assume two cases A and B for which this condition is held. Since the scaled initial conditions are invariant, the scaled solutions will be the same for these two cases only if the scaled WF advance velocity (i.e. dz_f^*/dt^*, where z_f^* is the scaled WF depth) is the same at each scaled time:

where subscripts A and B correspond to cases A and B. Eq. (8) gives:

here v_f [T^-1] represents dz_f^*/dt, and v_fi is the initial value of v_f. Based on Eq. (9), v_fi will be the best choice for the time scaling factor, Τ:

To determine v_fi, Darcy's equation is considered in the form of:

Integrating Eq. (11) from the soil surface to the WF depth at t = 0, we have:

which yields:

In Eqs. (12) and (13), h_i [L] is the matric potential head corresponding to θ_i, h_fi [L] and K_fi [LT^-1] are the matric potential head and hydraulic conductivity corresponding to θ_fi, respectively, θ_fi is the initial value of downward flux density at the WF, K_s [LT^-1] is the saturated hydraulic conductivity, and G [L] is WF matric head at t = 0 defined as:

Mass conservation law gives:

Combining (13) and (15) results in:

Using Eq. (16), v_fi can be calculated from the initial conditions (i.e. θ_fi, θ_i and z_fi). The soil hydraulic functions should be known to determine K_fi and h_fi from θ_fi as well as h_i from θ_i (h_fi and h_i are required to calculate G using Eq. (14)). Although there is no restriction on a specific form of the soil hydraulic models, it should be noticed that the selected hydraulic models should be the same for solving RE and calculating v_fi using Eq. (16). Based on the mass conservation law, ∂θ / ∂t = -∂q / ∂z , Eqs. (5), (6), and (10) suggest the following relationship to scale the flux density:

where q^* is the scaled flux density. To scale the redistribution process for a set of soils and initial conditions, q^*(z^*) should be a unified curve. Eq. (17) implies in another condition necessary for the scaling, dealing with the shape similarity of the flux density curves over the scaled soil profile, such that it can be scaled by a linear transformation. This kind of similarity was previously adopted by Simmons et al. (1979) and, as stated by Sposito and Jury (1985), is referred to as "Nielsen-similarity".

Materials and Methods

Richards' Equation was solved using HYDRUS-1D, Version 4 (Simunek et al., 2008) for various soils and initial conditions. To simulate the SWCPs during redistribution, zero water flux and free drainage were set as the upper and lower boundary conditions, respectively. Van Genuchten (1980) hydraulic functions based on Mualem's (1976) model were adopted:

where θ_r and θ_s are soil residual and saturated water contents, respectively, and α, n, and m are empirical parameters with the assumption that m=1-1/n.

Soils of twelve textural classes were considered by applying parameters of Carsel and Parrish (1988) which are the default parameters of HYDRUS for van Genuchten functions (Table 1). Two sets of initial conditions, presented in Table 2, were considered using various combinations of z_fi, h_fi and h_i. Set A (A1 to A6) was considered to evaluate the effect of the total water added to the soil profile, W, while Set B (B1 to B8) was considered to separately evaluate the effects of z_fi, h_fi and h_i. To prevent the divergence of the solutions, the abrupt increase of the matric potential head at the WF was avoided. To do so, the matric potential head gradually increased from h_fi and h_i through four space steps at the WF. It should be noted that HYDRUS works with the pressure head (negative values) instead of the absolute value of the matric head.

SWCPs obtained by HYDRUS were scaled using Eqs. (5), (6), and (10). The value of v_fi was determined from Eq. (16) with G being computed from Eq. (14). The integral of Eq. (14) was approximated by the trapezoidal rule. For convenience of application, it is worth to note that h_fi=0 and h_i à ∞, G approaches the effective capillary drive (H_cM) defined by Morel-Seytoux and Khanji (1974):

When van Genuchten hydraulic models with the assumption of m=1-1/n are applied, H_cM can be approximated as follows (Morel-Seytoux et al., 1996):

with minimal errors (< 2 %) over the range of m from 0.05 to 0.7. Eq. (21) was used to find a proper value for increments of h in the approximation of G using the trapezoidal rule. With 0.1 cm increments of h, approximations of the trapezoidal rule were close to those of Eq. (21).

To quantitatively evaluate the performance of the proposed scaling method, for each set of the scaled SWCPs, two criteria - the mean absolute error of scaling (MAES) and the root mean squared error of scaling of scaling (RMSES) - were defined as:

where θ^*_mi is the arithmetic mean of all θ^* values at each 0.1 increment of z^* at which θ^*_i is located, and N is the total number of the evaluated points. Profile of θ^*_m shows a mean scaled SWCPs, and therefore, MAES and RMSES indicate the deviations of the scaled SWCPs from the mean scaled SWCPs. When all the scaled SWCPs coalesce, these criteria will be equal to zero, suggesting an ideal performance of the scaling method.

Results and Discussion

Figure 2 shows HYDRUS model outputs for the 12 soils of Table 1 and initial conditions of A1 which were scaled by Eqs. (5), (6), and (10). Except for the sand, loamy sand and sandy loam, the remaining nine soils were reasonably well scaled and manifest a nearly unique scaled SWCP. A reason for the undesirable deviations in the three sandy soils is that the SWCPs in these soils with the imposed initial conditions are not similar in shape to those of the other soils. However, it seems that the main reason for the deviations is that, regarding the Nielsen-similarity condition, these soils are not similar to the other soils. To clarify this issue, the similarity condition should be evaluated.

As mentioned earlier, the so-called Nielsen-similarity deals with the shape similarity of the flux density curves during redistribution. It is assumed that the flux curves are described by the following power model (Jury and Horton, 2003):

where a, b, and c are empirical fitting parameters. When time t in Eq. (24) is longer than a day or so, the model shows a linear behavior. Therefore, the equality in the slope of the log-log plot of q(t) verifies the similarity.

We studied the similarity condition in three soils of sand, loam, and clay textures under the conditions corresponding to Figure 2. For these soils, log-log transformed curves of q(t) at z=z_fi (z^*=1) are shown in Figure 3. The slope of the sand is definitely greater than those of the loam and clay which are nearly more equal (Figure 3). Therefore, we consider that the curves for the loam and clay are nearly similar but different from that of the sand.

Figure 4 shows the scaled SWCPs for the 12 soils and initial conditions of A6 which considers a deeply wetted initial profile. A comparison of Figures 2 and 4 indicates that by increasing W (i.e. total water added to the soil profile) the scaling performance decreases. Since dz_f^*/dt at t=0 is used as the scaling factor of time, the scaled SWCPs increasingly diverge as the time increases. The poor performance of the scaling in Figure 4 can be similarly justified by invalidity of the two similarity conditions required for scaling. For the sand, loam, and clay of Figure 4, log-log transformed curves of q(t) at z=z_fi (z^*=1) are shown in Figure 5. The figure shows that the slopes are significantly different indicating that the Nielsen-similarity condition is not held in Figure 4. It can be concluded that by increasing W, the flow properties become more sensitive to soil texture. In other words, the textural range of the soils in which the similarity condition is held decreases as W increases.

The impact of W on the scaling performance, studied in detail by varying the initial conditions from A1 to A6, is shown in Figure 6 for the sand, loam, and clay at t^* = 5.

We conclude that the scaling is more sensitive to W for the very coarse- and very fine- textured soils. The scaled SWCPs for the loam are approximately invariant with respect to the initial conditions. For the sand, when W is relatively small (A1 and A2), the scaled SWCPs show a delay, while the scaling is adequate for larger values of W. On the other hand, the scaled SWCPs in the clay coalesced for all initial conditions except for A6 - the condition corresponding to the greatest value of W.

To separately study the effects of z_fi, h_fi and h_i on the scaling performance, RE was solved by applying set B of the initial conditions in which, for z_fi, h_fi and h_i, extreme values of set A are considered. The scaled solutions are presented in Figure 7 for sand, loam and clay at t^* = 5. The figure indicates that, in very coarse- or fine-textured soils, the scaling performance is more sensitive to z_fi, h_fi and h_i. The scaled SWCPs are approximately invariant for the loam while this is not the case for the sand and clay.

Impacts of z_fi, h_fi and h_i are dependent to each other. Based on the results, h_fi seems to be the most effective parameter on the scaling performance especially for fine-textured soils. Variation of this parameter significantly changes the flow rate and therefore, the shape of the flux curves. Thus, cases with an extreme value of h_fi will be far from the other cases regarding the Nielsen-similarity condition. Also, increasing z_fi decreases the scaling performance (Figure 7). This is the case for the sand when h_fi is large, however, for the clay when h_fi is small. The results confirm that the flow properties are not significantly affected by variation of h_i. In the sand, the impact of h_i is significant only if h_fi is large. This is the case for the clay if simultaneously, h_fi is large and z_fi is small.

Generally, the above discussions suggest that the proposed scaling method can be successfully applied for medium- and fine-textured soils provided that the initial profile is not deeply wetted. Figure 8 shows the scaled SWCPs at three scaled times of 1, 5, and 10 for the nine medium- and fine-textured soils of Table 1 (i.e. from loam to clay) having initial conditions A1 through A5. Even though the figure contains 45 various scenarios for the solution of RE, the scaled SWCPs are nearly invariant with limited scattering around the mean scaled SWCPs (the white points in the figure).

The 45 cases considered in Figure 8 are approximately Nielsen-similar. For these 45 cases, Figure 9a shows the log-log transformed curves of q(t) at z=z_fi (z^*=1). It is obvious that the slopes are approximately equal with an average of -0.955 (standard deviation=0.127). Therefore, the flux curves could be well scaled using Eq. (17), which are presented in Figure 9b. Although the flux values fall in a wide range, Figure 9b indicates that the scaled fluxes coalesced into a unified curve and could be well described using a power model similar to Eq. (24), q^*=0.71(t^*+0.069)^-1.03, with a determination coefficient of 0.995. Recognizing that the two assumed similarity conditions are held for the range of soils and selected initial conditions, we propose that the scaling method can be generalized to numerous other cases in this range leading to approximations of the solutions of RE.

Generalizing a Single Solution

The ability of the scaling method to approximate numerical solutions of RE using a single solution was evaluated. To do so, six medium- and fine-textured soils were selected from the literature: a Beit Netofa Clay (van Genuchten et al., 1980), a Pima Clay Loam (Elmaloglou and Malamos, 2003), and four other soils taken from UNSODA database (Leij et al., 1999) specified by codes 1300, 1370, 3360, and 4030. Some general information of these soils, including texture, taxonomic class and geographical location, are presented in le 3. For each soil, randomly produced initial conditions in the range of A1 to A5 were considered. To determine the van Genuchten parameters of the four soils of UNSODA, van Genuchten hydraulic models, Eqs. (18) and (19), were simultaneously fitted to the measured data using the RETC software (van Genuchten et al., 1991). Table 4 presents the van Genuchten parameters of the soils as well as the initial conditions.

It was assumed that, for these six soils, the scaled solutions of RE are invariant and identical to the mean scaled SWCPs of Figure 6. Subsequently, the mean scaled SWCPs were de-scaled (i.e. converted to the real scale) for the six soils using the Eqs. (5), (6) and (10). Hence, for each soil, SWCPs were approximated during redistribution. To evaluate the accuracy of the approximations, the SWCPs for the same soils and initial conditions were individually simulated using HYDRUS. The simulated and approximated SWCPs were compared using the mean relative error criterion, MRE (i.e. the mean of the absolute errors between the simulated and approximated values relative to the simulated values).

Figure 10 shows the simulated and approximated SWCPs for the six soils of Table 3 at two times corresponding to t^*=1 and t^*=10. With the MRE values never exceeding 5%, the scaling method satisfactorily approximated the numerical solutions of RE for the selected soils and initial conditions.

Conclusion

Scaled soil water content profiles were found to be nearly invariant during scaled redistribution times for medium- to fine-textured soils when the initial profile was not deeply wetted. An advantage of this method is that it is not restricted to a specific soil hydraulic model. A disadvantage is that the method does not adequately scale water content redistribution profiles of sands and other coarse-textured soils wetted partially, or those of fine-textured soils wetted deeply. Textural range of the soils in which the similarity conditions are held decreases as the initial conditions deals with a deeply wetted profile. In such a condition, a classification of soils and initial conditions and then separately scaling of each class may alleviate the problem. The method is promising to reduce complicated numerical calculations and opens a new window to easily obtain approximate solutions of highly nonlinear equation of Richards for water flow in unsaturated soils, within prescribed levels of error.

Acknowledgement

To Prof. D.R. Nielsen and Prof. A.W. Warrick, for their thorough review and helpful comments. The last author thanks CNPq for granting his productivity fellowship.

• Ahuja, L.R.; Williams, R.D. 1991. Scaling water characteristic and hydraulic conductivity based on Gregson-Hector-McGowan approach. Soil Science Society of America Journal 55: 308-319.
• Carsel, R.F.; Parrish, R.S. 1988. Developing joint probability distributions of soil water characteristics. Water Resources Research 24: 755-769.
• Elmaloglou, S.T.; Malamos, N. 2003. A method to estimate soil water movement under a trickle surface line source, with water extraction by roots. Irrigation and Drainage 52: 273-284.
• Jury, W.A.; Horton, R. 2003. Soil Physics. 6ed. John Wiley, New York, NY, USA.
• Kozak, J.A.; Ahuja, L.R. 2005. Scaling of infiltration and redistribution of water across soil textural classes. Soil Science Society of America Journal 69: 816-827
• Kosugi, K.; Hopmans, J.W. 1998. Scaling water retention curves for soils with lognormal pore-size distribution. Soil Science Society of America Journal 62: 1496-1504.
• Kutilek, M.; Nielsen, D.R. 1994. Soil Hydrology. Catena, Germany.
• Kutilek, M.; Zayani, K.; Haverkamp, R.; Parlange, J.Y.; Vachaud, G. 1991. Scaling of Richards' equation under invariant flux boundary conditions. Water Resources Research 27: 2181-2185.
• Leij, F.J.; Alves, W.J.; van Genuchten, M.T.; Williams, J.R. 1999. The UNSODA unsaturated soil hydraulic database. p. 1269-1281. In: van Genuchten, M.T.; Leij, F.J.; Wu. L., eds. Characterization and measurement of the hydraulic properties of unsaturated porous media. University of California, Riverside, CA, USA.
• Miller, E.E.; Miller, R.D. 1956. Physical theory for capillary flow phenomena. Journal of Applied Physics 27: 324-332.
• Morel-Seytoux, H.J.; Khanji, J. 1974. Derivation of an equation of infiltration. Water Resources Research 10: 795-800.
• Morel-Seytoux, H.J.; Meyer, P.D.; Nachabe, M.; Touma, J,; van Genuchten, M.T.; Lenhard, R.J. 1996. Parameter equivalence for the Brooks-Corey and van Genuchten soil characteristics: Preserving the effective capillary drive. Water Resources Research 32: 1251-1258.
• Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research 12: 513-22.
• Neuweiler, I.; Cirpka, O.A. 2005. Homogenization of Richards equation in permeability fields with different connectivities. Water Resources Research 41: W02009, doi:10.1029/2004WR003329.
• Reichardt, K.; Nielsen, D.R.; Biggar, J.W. 1972. Scaling of horizontal infiltration into homogeneous soils. Soil Science Society of America Proceedings 36: 241-245.
• Roth, K. 2008. Scaling of water flow through porous media and soils. European Journal of Soil Science 59: 125-130.
• Sadeghi, M.; Gohardous-Monfared, M.R.; Ghahraman, B. 2010. Scaling of soil hydraulic conductivity function using effective capillary drive. Journal of Water and Soil 24: 189-197 (in Persian).
• Simmons, C.S.; Nielsen, D.R.; Biggar, J.W. 1979. Scaling of field-measured soil-water properties. Hilgardia 47: 77-173.
• Simunek, J.; Sejna, M.; Saito, H.; Sakai, M.; van Genuchten, M.T. 2008. The HYDRUS-1D Software Package for Simulating the One-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media, Version 4.0. University of California, Riverside, CA, USA.
• Sharma, M.L.; Gander, G.A.; Hunt, C.G. 1980. Spatial variability of infiltration in a watershed. Journal of Hydrology 45: 101-122.
• Tuli, A.; Kosugi, K.; Hopmans, J.W. 2001. Simultaneous scaling of soil water retention and unsaturated hydraulic conductivity functions assuming lognormal pore-size distribution. Advance in Water Resources 24: 677-688.
• Van Genuchten, M.T. 1980 A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44: 892-898.
• Van Genuchten, M.T.; Leij, F.J.; Yates, S.R. 1991. The RETC Code for Quantifying the Hydraulic Functions of Unsaturated Soils. Robert S. Kerr Environmental Research Laboratory, Ada, OK, USA.
• Vereecken, H.; Kasteel, R.; Vanderborght, J.; Harter, T. 2007. Upscaling hydraulic properties and soil water flow processes in heterogeneous soils: a review. Vadose Zone Journal 6:1-28.
• Vogel, T.; Cislerova, M.; Hopmans, J.W. 1991. Porous media with linearly hydraulic properties. Water Resources Research 27: 2735-2741.
• Warrick, A.W.; Mullen, G.J.; Nielsen, D.R. 1977. Scaling of field measured hydraulic properties using a similar media concept. Water Resources Research 13: 355-362.
• Warrick, A.W.; Amoozegar-Fard, A. 1979. Infiltration and drainage calculations using spatially scaled hydraulic properties. Water Resources Research 15: 1116-1120.
• Warrick, A.W.; Lomen, D.O.; Yates, S.R. 1985. A generalized solution to infiltration. Soil Science Society of America Journal 49: 34-38.
• Warrick, A.W.; Hussen, A.A. 1993. Scaling of Richards' equation for infiltration and drainage. Soil Science Society of America Journal 57: 15-18.
• Youngs, E.G. 1990. Application of scaling to soil-water movement considering hysteresis. p. 23-37. In: Hillel, D.; Elrick, D.E., eds. Scaling in soil physics: principles and applications. Soil Science Society of America, Madison, WI, USA.

Publication Dates

• Publication in this collection
  03 Oct 2011
• Date of issue
  Oct 2011