Abstract

Richards equation is solved for soil water flow modeling in the unsaturated zone continuum. Interblock hydraulic conductivities, while solving for Richards equation, are estimated by some sort of averaging process based on upstream and downstream nodes hydraulic conductivity values. The accuracy of the interblock hydraulic conductivity estimation methods mainly depends on the distance between two adjacent discretized nodes. In general, the accuracy of the numerical solution of the Richards equation decreases as nodal grid discretization increases. Conventional interblock hydraulic conductivity estimation methods are mostly mere approximation approaches while the Darcian-based interblock hydraulic conductivities involve complex calculations and require intensive computation under different flow regimes. Therefore, in this study, we proposed an effective saturation-based weighting approach in the soil hydraulic curve functions for estimating interblock hydraulic conductivity using a one-dimensional vertical finite-difference model which provides a parametric basis for interblock hydraulic conductivity estimation while reducing complexity in the calculation and computational processes. Furthermore, we compared four test case simulation results from different interblock hydraulic conductivity methods with the reference solutions. The comparison results show that the proposed method performance in terms of percentage reduction in root mean square and mean absolute error over other methods compared in this study were 59.5 and 60%, respectively.

INTRODUCTION

Water flows in a variably saturated or unsaturated zone can be modeled by Richards equation (Richards 1931; Hillel 1980). The Richards equation can be solved in three different forms: (i) head-based, (ii) moisture content-based, and (iii) mixed form of head and moisture content-based. The mixed form of Richards equation reduces the mass balance error when compared to head-based form of the Richards equation (Celia et al. 1990). The solution of the Richards equation can be obtained by analytical methods for simplified boundary conditions (Srivastava & Yeh 1991; Basha 1999; Chen et al. 2001; Yuan & Lu 2005; Menziani et al. 2007; Tracy 2007; Zlotnik et al. 2007; Hayek 2014, 2015, 2016; Pugnaghi et al. 2015). However, the complex nonlinear Richards equation limits the analytical solution for complex boundary conditions such as alternate wetting and drying taking place due to rainfall/irrigation and evaporation processes at the soil surface, for example. Thus, various numerical approaches have been adopted by many researchers to solve the highly nonlinear Richards equation under complex and natural boundary conditions (Keese et al. 2005; Kolditz et al. 2008; Essig et al. 2009; Jiménez-Martínez et al. 2009; Šimůnek et al. 2009; Wang et al. 2009; Broda et al. 2011; Carrera-Hernandez et al. 2012; Namin & Boroomand 2012; Mohanasundaram et al. 2013; Herrada et al. 2014; Arrey et al. 2017). Furthermore, a complete review of the numerical solution to Richards equation has been done by Farthing & Ogden (2017) and Zymkiewicz (2013) and the historical background of Richards equation modelling with respect to catchment scale hydrology can be found in Paniconi & Putti (2015).

The soil-water retention curve (SWR) depicts the relationship between moisture content and pressure head variations in the unsaturated soil system. Similarly, the soil hydraulic conductivity curve (SHC) depicts the relationship between hydraulic conductivity and pressure head or water content. As the Richards equation is a complex nonlinear equation, it has to be solved using these additional SWR and SHC functions. The SWR and SHC can be modeled with different parametric-based analytical relationships (Brooks & Corey 1964; Mualem 1976; Haverkamp et al. 1977; van Genuchten 1980). The parameters of these SWR and SHC functions are estimated directly from experimental data points using different optimization algorithms (Maggi 2017; Londra & Kargas 2018; Wang et al. 2018), RETC computer program (van Genuchten et al. 1991) or using neural network-based pedotransfer functions (Rosetta) as adopted in the HYDRUS-1D package (Schaap et al. 2001; Šimůnek et al. 2009).

Interblock hydraulic conductivity or average hydraulic conductivity is computed while solving the Richards equation through an iteration process. There are many methods presently available for calculating interblock hydraulic conductivity values such as arithmetic averaging, geometric averaging, harmonic averaging and upstream mean methods (Haverkamp & Vauclin 1979; Oldenburg & Pruess 1993; Romano et al. 1998; van Dam & Feddes 2000). As these conventional interblock hydraulic conductivity estimation functions are simple approximation-based functions, they tend to deviate from the actual numerical solutions in the unsaturated zone flow problems and especially the deviation error propagates larger as the nodal discretization increases. Therefore, the present study introduces an effective saturation-based weighting concept to calculate interblock hydraulic conductivity values for accurately simulating one-dimensional (1D) vertical flows in the unsaturated soils under a variety of soil-water flow problems with specific emphasis on different grid discretizations.

Traditional interblock hydraulic conductivity estimation functions use some sort of averaging process to estimate interblock hydraulic conductivity values at the interface between upstream and downstream nodes. This is mainly to reduce the computational complexities in the numerical schemes while solving for complex Richards equation. However, the exact interblock hydraulic conductivity can be obtained by assuming a steady-state flow between two adjacent nodes based on Darcy's law (Warrick 1991; Baker 1995, 2000). Some studies approximate the pattern of steady-state solution between adjacent nodes in terms of assigning certain weights to upstream and downstream nodes hydraulic conductivities (Gastó et al. 2002). The concept of upstream mean method evolved based on the condition that the total head (suction/pressure head + elevation head) at the upstream node dominates over the downstream node value (Oldenburg & Pruess 1993). In general, the upstream mean method produces numerical results closer to the actual results under infiltration driven flow conditions and this method even gives a reasonably accurate solution under coarser grid discretization (Szymkiewicz 2009). The integrated mean estimates the approximate interblock hydraulic conductivity at the interface based on integrating the upstream and downstream nodes hydraulic conductivity values with respect to the distance between the respective nodes (Zaidel & Russo 1992). The Darcian-based mean approaches can be adopted for different flow regimes such as infiltration, drainage and capillary rise segments in soil-water flow systems (Szymkiewicz 2009). The Darcian-based mean method is where the interblock hydraulic conductivity is estimated between the range of integrated mean and upstream mean values (Szymkiewicz 2009). The recent study on interblock hydraulic conductivity proposes a new function based on a higher order upwind mean method in which the interblock hydraulic conductivity is approximated with second-order spatial accuracy (An & Noh 2014).

After implementing one of the conventional interblock hydraulic conductivity functions for fine grid discretizations, numerically solved pressure head and water content values reasonably match with the actual numerical solutions. However, when the grid size discretization increases, the accuracy of the numerical solution reduces. This is evident from many studies. For example, arithmetic mean-based interblock hydraulic conductivity overestimates the soil-water fluxes while geometric mean-based interblock hydraulic conductivity underestimates the water fluxes at coarser grid discretization (van Dam & Feddes 2000). However, the Darcian integral method preserves the accurate numerical results which consider Darcy law and assumes the steady-state flux between the adjacent nodes (Warrick 1991; Baker 1995; Szymkiewicz 2009). The comparison of many simple averaging methods against Darcian-based approaches at coarser grid discretization show large variations in the simulated and actual numerical solutions due to highly nonlinear behavior of the hydraulic conductivity – pressure head functions under different soil textures and water regimes (Haverkamp & Vauclin 1979; Gastó et al. 2002; Belfort & Lehmann 2005; Szymkiewicz 2009; Szymkiewicz & Burzyński 2011; Belfort et al. 2013; An & Noh 2014). Nevertheless, the Darcian-mean methods have some limitations in terms of numerous computation processes and unavailability of approximated weighting functions for different soils to compute the interblock hydraulic conductivity values. Even though the Darcian mean methods are reasonably accurate under wide textural classes of homogeneous soil, it complicates the calculation of interblock hydraulic conductivity under heterogeneous layered soil systems (Belfort et al. 2013).

After rigorous application of these standard interblock hydraulic conductivity functions in the homogeneous porous systems, some studies have attempted to assess the numerical accuracy in the heterogeneous layered soil systems (Szymkiewicz & Helmig 2011; An & Noh 2014). Some studies conclude that the arithmetic averaging methods, although they are simple in implementing within the main numerical algorithms, cause huge stability and oscillation problems in coarse grid simulation conditions (Szymkiewicz 2009; Belfort et al. 2013). The Darcian approach-based methods are accurate, although they require more computational efforts to be implemented within the main numerical scheme, they may be challenging to obtain optimum solutions under complex multi-dimensional problems (Belfort et al. 2013). The upwind mean method can be an alternative method when the arithmetic and geometric averaging methods violate the monotonicity conditions while the upwind mean method preserves the monotonicity of the numerical solution for the Richards equation (Belfort et al. 2013).

The standard interblock hydraulic conductivity functions are easy to implement in the numerical schemes for solving Richards equation but the accuracy of the numerical solution decreases when the grid size discretization increases for these methods. Conversely, Darcian-based methods yield accurate results but they are relatively complex functions and require additional supporting information for solving Richards equation under different water flow regimes. Therefore, in the present study, we have proposed an effective saturation-based weighting for calculating interblock hydraulic conductivity which will be an alternative method in terms of producing relatively accurate numerical solutions, even with coarser grid discretizations, and which is simple to implement in the numerical schemes. The computed effective saturation values at the discretized nodes are used as the weighting fractions in the SHC functions for calculating interblock hydraulic conductivities in the proposed approach. The specific objectives of the present work are: (1) to formulate an effective saturation-based weighting in the SHC functions to calculate the interblock hydraulic conductivity for solving Richards equation using 1D vertical finite-difference model, and (2) to compare and assess the accuracy of the proposed effective saturation-based interblock hydraulic conductivity function against arithmetic mean and higher-order mean methods under (i) homogeneous and heterogeneous soils, (ii) different initial and boundary conditions, and (iii) different vertical grid discretizations.

The reference solution resulting from a finely discretized vertical soil column with 1 mm nodal spacing between two adjacent nodes was simulated from HYDRUS 1D model for comparison purposes.

METHODOLOGY

Discretization of the Richards equation

In 1D finite-difference-based numerical modelling methods, the numerical solutions of the dependent variables (pressure head or moisture content) are solved at the grid nodes where the nodes are either equally spaced or non-equally spaced. In the present study, we considered only the equally spaced 1D vertical finite-difference model for simulating the soil water flows. The formulation of a mixed form of Richards equation for solving 1D vertical flow in the variably saturated zone is as follows:
formula
(1)
where θ is moisture content (L3 L–3); h is suction pressure head or pressure head (L); K(h) is unsaturated hydraulic conductivity of the soil which is a function of pressure head (LT–1); z is vertical coordinate assumed positive downward (L); t is time (T).
For simplicity, the sink term is not considered in this study. The finite-difference discretization for fully implicit backward Euler formulation of Equation (1) with a constant nodal discretization is given as follows:
formula
(2)
where n + 1 is new time step; n is the current time step; m + 1 is new iteration level; m is current iteration level; is interblock hydraulic conductivity (LT–1); Δt is time interval (T); Δz is grid discretization (L).

Modified Picard iteration method for solving Richards equation

The numerical solution of Equation (2) can be achieved in an iterative manner with a pre-determined threshold error limit. Celia et al. (1990) proposed a modified Picard iteration method which improved the mass balance property of the numerical solution. The modified Picard iteration technique can be carried out as follows.

Let us consider Equation (2) and expand the term with respect to the term using a truncated Taylor series approximation. The Taylor series approximation for the new time step water content variable for the ith node () can be expressed as follows:
formula
(3)
Furthermore, the increment of the dependent variable from iteration m to m + 1 can be written as follows:
formula
(4)
Rearranging Equation (4) in terms of the dependent variable at new time step (n + 1) is as follows:
formula
(5)
The term in Equation (3) denotes the specific capacity of the soil which can be denoted as . Furthermore second and higher order terms in Equation (3) can be neglected due to smaller values. Substituting Equations (3)–(5) in Equation (2), we obtain:
formula
(6)
Rearranging Equation (6) in terms of new time step – new iteration levels (n + 1, m + 1) in the left hand side and new time step – current iteration levels (n + 1, m) in the right hand side, we obtain:
formula
(7)
Formulating tri-diagonal matrix for the system of linear equations, we obtain:
formula
(8)
where [P]n+1,m is a tri-diagonal matrix containing the coefficients corresponding to i − 1, i, i + 1 nodes; is the right-hand side known vector values of the system of linear equations.
The iteration process continues until the error tolerance limit or absolute error is reached between the two consecutive iterations as follows:
formula
(9)
Furthermore, this absolute error convergence criterion was improved by Huang et al. (1996) and the corresponding formulation is given as follows:
formula
(10)
where γa is the absolute error; Cn+1,m is the specific capacity of the soil (L–1); γθ is the threshold limit for the water content which is usually equal to 0.001.

As this convergence criterion proposed by Huang et al. (1996) significantly improved the computational time, we implemented this convergence criterion in our study for solving the Richards equation. However, the focus of the present study is not in the convergence criteria to improve the computational efficiency of the numerical scheme, rather we implemented a new effective saturation-based interblock hydraulic conductivity estimation function in the numerical scheme to improve the accuracy of the numerical solutions of the Richards equation for a wide range of soil types, boundary conditions, and nodal discretizations.

Parametric soil-water retention functions

The SWR and SHC functions are used as an additional constitutive relationship functions while solving for Richards equation. Many parametric SWR and SHC models were developed in the past with unique fitting parameters (Mualem 1976; Haverkamp et al. 1977; van Genuchten 1980). For simplicity, we considered the van Genuchten–Mualem model in this study (Mualem 1976; van Genuchten 1980). The functional relationships among moisture content, pressure head and hydraulic conductivities based on the van Genuchten–Mualem model are as follows:
formula
(11)
formula
(12)
formula
(13)
formula
(14)
where l is the pore connectivity parameter which is generally taken to be 0.5 as an average value for many soils (Mualem 1976); n is the shape parameter.

Interblock hydraulic conductivity functions

The numerical solution for the dependent variables is solved only at the discretized nodes and the solutions at the interface between the discretized nodes are generally not known from the finite differencing schemes. Therefore, the interblock hydraulic conductivities are estimated based on the adjacent nodes hydraulic conductivity values. Different weighting schemes are used based on the adjacent nodes hydraulic conductivity values for estimating an effective hydraulic conductivity or interblock hydraulic conductivity values (Haverkamp & Vauclin 1979; Celia et al. 1990; Oldenburg & Pruess 1993). Different formulations for computing interblock hydraulic conductivity values are discussed below.

Arithmetic mean

The interblock hydraulic conductivity values in each iteration within each time step can be calculated by a simple arithmetic mean method (KAM) as follows (Celia et al. 1990):
formula
(15)

Geometric mean

The interblock hydraulic conductivities for the upstream and downstream nodes can also be estimated by the geometric mean method (KGM) which is given as follows (Haverkamp & Vauclin 1979):
formula
(16)

Harmonic mean

The harmonic mean (KHM) method of estimating the interblock hydraulic conductivities is as follows (Oldenberg & Pruess 1993; Romano et al. 1998):
formula
(17)

Upwind mean

The upstream mean (KUP) method of estimating the interblock hydraulic conductivity function is as follows (Oldenburg & Pruess 1993):
formula
(18)
where H=h+z is the total head (L)

Higher-order mean method

An & Noh (2014) adopted a scheme based on a higher-order mean approach for estimating interblock hydraulic conductivity in which two different values of hydraulic conductivity at the cell interface was estimated depending on the gradient value in the cell. The formulation for estimating the higher-order mean hydraulic conductivity (An & Noh 2014) is as follows:
formula
(19)
formula
(20)
formula
(21)
where ri= ΔKi+1/2Ki−1/2, ΔKi−1/2 = KiKi−1, ΔKi+1/2 = Ki+1Ki, and ϕ(r) = max (0, min(r,1)) is the slope limiter function.

It is to be noted that if the gradient K is zero, this higher-order upwind mean method becomes a traditional upwind mean method (Equation (18)).

Proposed effective saturation weighting method

A schematic diagram of the proposed method is shown in Figure 1. In this study, we used the well-established parametric-based SHC functions for calculating the interblock hydraulic conductivity values instead of empirical-based weighting functions. Any well-established SHC functions can be used to obtain the interblock hydraulic conductivity values using this proposed approach. For brevity, we considered only the van Genuchten–Mualem model (Equations (11) and (12)) in this study. Equations (11) and (12) denote the SWR and SHC functions respectively. We know that the SWR and SHC functions are used with the Richards equation to obtain numerical solutions at the discretized nodal points. However, the SHC functions (Equation (12)) as such will not be used at the interface location to compute the interblock hydraulic conductivities because the effective saturation value at the interface is not known. Therefore, we used a slightly modified SHC function at the interface between the corresponding adjacent nodes where we replaced the effective saturation term with an equivalent effective saturation term to account for the average value of the effective saturation at the interface. The equivalent effective saturation term at the interface was computed based on the adjacent nodes effective saturation values (Figure 1). The proposed effective saturation-based interblock hydraulic conductivity function is given as follows:
formula
(22)
formula
(23)
formula
(24)
where θi is the water content for ith node at current time level; θr is residual water content of the soil; θs is saturated water content of the soil; Si is the effective saturation at the ith node; Si+1/2, Si−1/2, Ki+1/2, and Ki−1/2 are interface effective saturation and hydraulic conductivity values between i and i + 1, and i and i − 1 nodes respectively.
Figure 1

Schematic diagram of the proposed effective saturation-based weighting method for interblock hydraulic conductivity estimation.

Figure 1

Schematic diagram of the proposed effective saturation-based weighting method for interblock hydraulic conductivity estimation.

The average effective saturation values at the interface are not solved in the finite difference based numerical schemes. So, the equivalent effective saturation values are to be estimated based on adjacent nodes effective saturation values. The equivalent effective saturation values (Si+1/2, Si−1/2,) at the interface either between the ith and i + 1th nodes or ith and i−1th nodes are calculated in such a way that the weights are properly normalized for achieving the stability of the numerical solution (Equation (23)). For example, the following calculation with assumed values confirms the proper estimation of the interblock hydraulic conductivity values at the interface for the numerical stability criteria.

Let us consider a soil column with four equally spaced vertical nodes from top to bottom (i = 1, 2, 3 and 4). The three interface locations in the 4th nodal soil column can be located at , , and positions respectively. Let us consider a sand soil and corresponding van Genuchten model parameters as required in Equation (22) as follows:
formula
Let us also assume the initial effective saturation values at the nodes i= 1, 2, 3 and 4 as follows:
formula

Computing interblock hydraulic conductivities using Equation (22) with the assumed values as mentioned above at the interface Ki+1/2 and Ki−1/2 is as follows:

For node i = 2:
formula
formula
formula
formula
Similarly, for node i = 3:
formula
formula

It can be observed from the above sample calculation that the interblock hydraulic conductivities upstream () and downstream () sides of the ith nodes are calculated properly. For instance, if we consider the interface between nodes 2 and 3, the resulting interblock hydraulic conductivity should be calculated at . The sample calculation results reveal that the downstream side interblock hydraulic conductivity value for the node i = 2 matches precisely with the upstream side interblock hydraulic conductivity value for the node i = 3. This situation can be further visualized whereby the interface location at which the estimated two different interblock hydraulic conductivities ( and ) are nothing but the same. Therefore, the values must be the same. This confirms the stable numerical solutions by the proposed method and thus this approach can be adopted for solving the Richards equation.

Briefly, the significance of the proposed method can be recognized in two aspects: (1) it models the interblock hydraulic conductivity using the specific SHC functions such as the van Genuchten–Mualem (Mualem 1976; van Genuchten 1980) and Brooks–Corey (Brooks & Corey 1964) models which are basically a parametric-based SWR and SHC functions thus controlling the smooth variation of the estimated interblock hydraulic conductivities rather than a simple approximation by the conventional methods, and (2) it directly uses the parametric-based SHC functions to estimate interblock hydraulic conductivity values, thus complex calculations as required by the Darcian-based approaches are not necessary under any flow conditions.

Test case simulations

Effectiveness of the developed interblock hydraulic conductivity methods was assessed by simulating various test case scenarios with different soils, boundary conditions, and grid discretizations. For brevity and clarity, only three methods of calculating interblock hydraulic conductivities, such as arithmetic average, higher-order mean and the proposed effective saturation methods, were compared in each scenario simulated in this study. We specifically chose the higher-order mean method to compare with the proposed method results because the higher-order means method proposed by An & Noh (2014) was the recent study which compared the relevant standard and recently evolved interblock hydraulic conductivity functions and their proposed model results were better than the compared model results.

Accuracy assessment for the numerical solutions

The numerical accuracy of different interblock hydraulic conductivity estimation functions developed in the present study was compared against the reference numerical solution. The reference numerical solution was simulated using the HYDRUS 1D model with finely discretized soil profiles (1 mm) such that all interblock hydraulic conductivity functions yield the same results. Two statistical indices, mean absolute error (MAE) and root mean square error (RMSE), was adopted to assess the performance of the developed interblock hydraulic conductivity functions. The formulations of MAE and RMSE are as follows:
formula
(25)
formula
(26)
where Xref is the reference pressure head or water content values at a specific location and time step (L or L3 L–3); Xi is the modeled pressure head or water content values at the corresponding location and time step (L or L3 L–3); N is the total number of data points.
Furthermore, an additional moisture content error index (EMC) based on simulated and reference solution moisture content values according to An & Noh (2014) was also used to compute the performance of the developed models which is given as follows:
formula
(27)
where θref is the reference solution moisture content (L3 L–3).

TEST CASE SIMULATIONS RESULTS AND DISCUSSION

Test case 1: infiltration under constant head boundary condition in homogenous soil

Test case 1 simulated a vertical 1D flow in a homogenous soil profile under constant pressure head boundary conditions with the value of hz=0,t = –1 cm applied at the soil surface node (z = 0) and at all time steps (t > 0) throughout the simulation time. The total depth of the soil profile (L) was considered as 5 m. The initial condition was assigned with the pressure head value of hz,t=0 = –500 cm at time t = 0. The bottom boundary condition was assigned with a free drainage type boundary condition (hz=L,t = hz=L−1,t). The soil type used in this test case and the corresponding SWR parameters based on Carsel & Parrish (1988) are given in Table 1. Simulations were carried out for four different vertical grid discretizations (Δz) of 5, 10, 25 and 50 cm respectively. The total simulation time considered for this test case was about 2 hours.

Table 1

Soil types and the corresponding soil-water retention parameters (Carsel & Parrish 1988) used in the test cases 1, 2, 3, and 4

Soil-water retention parameters for Mualem-van Genuchten modelTest case 1Test case 2Test case 3Test case 4
Layer 1Layer 2Layer 3
Soil texture Sand Loam Clay loam Loam Silt loam Clay 
θs (cm3/cm30.43 0.43 0.41 0.43 0.45 0.38 
θr (cm3/cm30.045 0.078 0.095 0.078 0.067 0.068 
α (1/cm) 0.145 0.036 0.019 0.036 0.02 0.008 
n 2.68 1.56 1.31 1.56 1.41 1.09 
Ks (cm/hr) 29.7 1.04 0.26 1.04 0.45 0.2 
l 0.5 0.5 0.5 0.5 0.5 0.5 
Soil-water retention parameters for Mualem-van Genuchten modelTest case 1Test case 2Test case 3Test case 4
Layer 1Layer 2Layer 3
Soil texture Sand Loam Clay loam Loam Silt loam Clay 
θs (cm3/cm30.43 0.43 0.41 0.43 0.45 0.38 
θr (cm3/cm30.045 0.078 0.095 0.078 0.067 0.068 
α (1/cm) 0.145 0.036 0.019 0.036 0.02 0.008 
n 2.68 1.56 1.31 1.56 1.41 1.09 
Ks (cm/hr) 29.7 1.04 0.26 1.04 0.45 0.2 
l 0.5 0.5 0.5 0.5 0.5 0.5 

Simulated pressure head profiles of different interblock hydraulic conductivity methods for the vertical discretization of 25 cm is shown in Figure 2. Arithmetic and higher-order mean methods show relatively deviated pressure head profiles, especially in the downstream side of the reference solution when compared to the proposed mean method. This is because arithmetic and higher-order mean methods estimate the interblock hydraulic conductivity by directly using upstream and downstream nodes hydraulic conductivity values thus resulting in overestimation of the interblock hydraulic conductivity which caused deviated results. However, the proposed mean method, unlike other methods, imposed the weights based on actual effective saturation values through the SHC functions for calculating interblock hydraulic conductivity values which resulted in a reasonably accurate estimation of interblock hydraulic conductivities and so the pressure heads and moisture contents.

Figure 2

Pressure head distribution profile after 2 hours of simulation under constant head boundary conditions (hz=0,t= − 1 cm) in sand soil for the grid discretization of 25 cm.

Figure 2

Pressure head distribution profile after 2 hours of simulation under constant head boundary conditions (hz=0,t= − 1 cm) in sand soil for the grid discretization of 25 cm.

The accuracy of the simulated pressure head distribution profiles with respect to the reference solution for different grid discretizations was calculated using Equations (25) and (26). The numerical results for different grid discretizations show that the accuracy of the simulated results from all developed methods deviated more and more from the reference solution as the grid discretizations increased. The calculated model performance indices (RMSE and MAE) for the test case 1 simulations are shown in Figure 3. The proposed method was the relatively least erroneous method when compared to other methods at all grid discretizations. This is due to the fact that the proposed method calculated the interblock hydraulic conductivity based on effective saturation weighting and through proper parametric-based SHC functions thus reasonably estimating the actual interblock hydraulic conductivities at the interface.

Figure 3

Accuracy of the simulated pressure head values with respect to the reference solution in terms of (a) RMSE and, (b) MAE values after 2 hours for test case 1 simulations.

Figure 3

Accuracy of the simulated pressure head values with respect to the reference solution in terms of (a) RMSE and, (b) MAE values after 2 hours for test case 1 simulations.

Test case 2: constant head infiltration and constant head bottom boundary conditions in the homogeneous soil

This test case simulated the flow in homogenous loam soil under constant head top boundary (infiltration) and constant head bottom boundary conditions by imposing the constant pressure head value of hz=0,z=L,t = –1 cm at both the top and bottom nodes for all time steps. The initial condition was assigned as hz,t=0 = –100 cm to all the discretized nodes in the soil domain at time t = 0. The total length of the soil profile was assigned to 100 cm. Similarly, the vertical grid discretizations were varied with 5, 10 and 15 cm. The simulation was continued for up to 2 hours after which the results were compared with the reference solutions. Furthermore, the parameters of SWR and SHC functions of the soil used in this test case are listed in Table 1.

The simulated pressure head distribution profiles from different interblock hydraulic conductivity methods for the grid discretization of 10 cm are shown in Figure 4. This simulation shows the water movement in the soil both due to infiltration (downward movement) from the top and capillary rise from the bottom (upward movement) of the soil profile. It can be observed from Figure 4 that the pressure head profiles from the top and bottom node advanced relatively faster with respect to the reference solution at the end of simulation time (t = 2 h) for arithmetic and higher-order mean methods. This indicates that the arithmetic and higher-order mean methods under initially wet soil profiles (hz,t=0 = –100 cm), as in this case, computed the interblock hydraulic conductivity in such a way that it thrusts the water movement to relatively higher depths than as compared to the previous test case (test case 1) condition where the initial soil profile was considered relatively dry (hz,t=0 = –500 cm). Therefore, both arithmetic and higher-order mean methods under relatively wet initial soils perform poorer than dry initial soils. However, the proposed mean method estimated the interblock hydraulic conductivity through parametric SHC functions and the resulted solution was closely matched with the reference solutions under different initial conditions and at larger grid discretizations as well (Figures 3 and 4).

Figure 4

Pressure head distribution after 2 hours of simulation under constant head top and bottom boundary conditions (hz=0,z=L,t = –1 cm) in loam soil for the grid discretization of 10 cm.

Figure 4

Pressure head distribution after 2 hours of simulation under constant head top and bottom boundary conditions (hz=0,z=L,t = –1 cm) in loam soil for the grid discretization of 10 cm.

The computed RMSE and MAE values for these test case simulations are shown in Figure 5(a) and 5(b) respectively. The arithmetic and higher-order mean methods are on a par with each other while the proposed mean method was the least erroneous method among other methods. Overall, the maximum magnitude of the calculated RMSE and MAE values for all the developed models were within 40 and 25 cm respectively. This is because this test case was simulated for relatively wet soil with the initial pressure head distribution of –100 cm, thus the resulted deviations were calculated within the range of –1 to –100 cm. Nevertheless, the proposed method performed significantly well within the RMSE and MAE values of 15 and 10 cm, respectively.

Figure 5

Calculated (a) RMSE and (b) MAE values after 2 hours for the test case 2 simulations.

Figure 5

Calculated (a) RMSE and (b) MAE values after 2 hours for the test case 2 simulations.

Test case 3: atmospheric boundary condition in homogenous soil

This test simulated the flow under alternate wetting and drying spells of precipitation through the atmospheric boundary at the top node and the free drainage boundary at the bottom node of the soil profile. The specified atmospheric boundary with precipitation and evaporation rates for four consecutive days are listed in Table 2. This atmospheric boundary condition removes the excess water which is resulting from oversaturation as surface runoff from the surface. Similarly, the potential evaporation rate is restricted to a limited evaporation rate whenever the simulated pressure head at the top node exceeds the specified critical pressure head value. For this simulation, the critical pressure head value was assigned with hz=0 = –15,000 cm. It was assumed in this simulation that if the simulated pressure head value at any time during the simulation exceeded the wilting point equivalent (hz=0 = –15,000 cm) pressure head value, the potential evaporation rate was reduced to a limited evaporation rate. The initial condition was assigned as hz,t=0 = –100 cm. Three different vertical discretizations were simulated as 5, 10 and 20 cm respectively. The total time of simulation was approximately 4 days as mentioned in Table 2.

Table 2

Time-varying atmospheric boundary conditions used to simulate the flow in test case 3

DayPrecipitation (cm/days)Evaporation (cm/days)Critical pressure head (cm)
−15,000 
−15,000 
−15,000 
−15,000 
DayPrecipitation (cm/days)Evaporation (cm/days)Critical pressure head (cm)
−15,000 
−15,000 
−15,000 
−15,000 

Figure 6(a) and 6(b) shows the simulated moisture content profiles after the 3rd and 4th days of simulation. The moisture content profiles were simulated for alternate wetting and drying spells of precipitation events in this test case. The simulated reference moisture content profile on the 3rd day shows that the soil profile was completely saturated up to 100 cm below the soil surface due to a precipitation event of 6 cm/day. However, the developed interblock hydraulic conductivity models for 10 cm grid discretization predict the simulated saturated water content profiles at about 60–70 cm below the surface. As we already saw from the other two test cases (test cases 1 and 2), the accuracy of the numerical solution decreased as the vertical grid discretizations increased. Therefore, the simulated moisture content profiles by various interblock hydraulic conductivity methods show relatively larger deviated profiles from the reference solution. In addition to the coarser vertical discretization factor, the length of the simulation time, which in this case was 3 days, was also relatively longer and thus caused the deviation in the simulated moisture content profiles from reference solutions (Figure 6(a)). Similarly, the simulated moisture content profiles after 4 days also show the deviated profiles from the reference profiles (Figure 6(b)). After a precipitation event on the 3rd day, there was a drying spell on the 4th day which resulted in the upward movement of water through an evaporation process from the top layer. As a result, it significantly removed water from the soil profile and brought down the moisture content of the top node from 0.41 to 0.15 cm3/cm3.

Figure 6

Moisture content distribution profiles under atmospheric boundary conditions for the grid discretization of 10 cm after (a) 3 days, and (b) 4 days of the simulation.

Figure 6

Moisture content distribution profiles under atmospheric boundary conditions for the grid discretization of 10 cm after (a) 3 days, and (b) 4 days of the simulation.

The three vertical discretizations simulation results with respect to the reference solution were computed in terms of RMSE, MAE and EMC values after 3 and 4 days of the simulation are shown in Figures 7 and 8, respectively. The performance of the proposed method was relatively better than other methods as the calculated RMSE, MAE and EMC values were relatively lesser when compared to other methods (Figures 7 and 8). An additional error index for the moisture content was also computed for these test case simulations based on the deviation of the simulated and reference-based moisture content values with respect to the reference-based moisture content values. This error index (EMC) was also computed as the lowest for the proposed method when compared to other methods. This shows that the proposed interblock hydraulic conductivity method performs relatively better among all other methods and across all the indices such as RMSE, MAE and EMC. Although the calculated RMSE, MAE and EMC values for 4 days after simulation were on a par with each other for all the methods, the proposed method was nevertheless showing relatively lower values across all indices (Figure 8).

Figure 7

Calculated (a) RMSE, (b) MAE, and (c) EMC values after 3 days for the test case 3 simulations.

Figure 7

Calculated (a) RMSE, (b) MAE, and (c) EMC values after 3 days for the test case 3 simulations.

Figure 8

Calculated (a) RMSE, (b) MAE, and (c) EMC values after 4 days for the test case 3 simulations.

Figure 8

Calculated (a) RMSE, (b) MAE, and (c) EMC values after 4 days for the test case 3 simulations.

Test case 4: heterogeneous layered soil with constant head boundary conditions

Test case 4 simulated the flow in the heterogeneous three-layered soil system under constant head boundary conditions at the top and bottom nodes. Soil type and corresponding SWR parameters of the three soils are listed in Table 1. The depth of each soil layer was 30 cm and thus the total depth of the simulated soil profile was 90 cm. The initial condition was assigned with the pressure head value of hz,t=0 = –500 cm. The top and bottom boundary condition was assigned as a constant head type boundary condition with the pressure head value of hz,=0,t = –1 cm. The total simulation time was about 8 hours and three different vertical discretizations were simulated with 2, 5, and 15 cm respectively.

The simulation results of the heterogenous three-layered soil system are shown in Figure 9. The simulated reference moisture content profiles after 8 hours under constant head boundary conditions show that the profile fully saturated the entire top layer soils and a few cms of the middle layer soils (Figure 9). Conversely, the upward movement of water from the bottom-most node saturated more than 50% depth of the bottom layer soils. It was clearly seen that the arithmetic and higher-order mean methods predict the moisture content profiles away from the reference solution profiles for the vertical discretization of 5 cm (Figure 9). The simulated moisture content profile from the proposed method was relatively better than other methods and closely matched with the reference solution simulated by the HYDRUS 1D model. It has already been shown that the proposed method performs relatively well with the fine-grained soils with constant head-based boundary conditions (infiltration under ponded water) in previous test cases. This test case also demonstrated the soil-water movement under constant head boundary conditions in fine-grained soils arranged in three-layered soil systems and the performance of the proposed method was relatively better when compared to other methods. The simulated moisture content profiles at the interface between two layers (layers 1 and 2) showed oscillated profiles for both arithmetic and higher-order mean methods (Figure 9). This could be due to improper computation of interblock hydraulic conductivities resulting from unique SWR and SHC parameters of two different soils at the layer interface. However, the proposed method computed the numerical solution at the interface in such a way that the moisture content profile was relatively smooth and continuous due to the fact that it adopted parametric-based SHC functions and equivalent effective saturation values to compute the interblock hydraulic conductivities at the interface.

Figure 9

Simulated moisture content distribution profiles after 8 hours under constant head boundary conditions in heterogeneous three-layered soil system for the grid sizes of 5 cm.

Figure 9

Simulated moisture content distribution profiles after 8 hours under constant head boundary conditions in heterogeneous three-layered soil system for the grid sizes of 5 cm.

Figure 10 showing the simulated soil moisture content values after 8 hours with respect to the reference-based solutions were relatively better for the proposed method when compared to other methods. For example, the computed EMC value for the proposed method was reduced by 50 and 30% from the arithmetic mean and higher-order mean methods respectively (Figure 10). It was also observed that the reduction in MAE values was increased as the grid discretization increased for the proposed method when compared to other methods. This was mainly due to the fact that the proposed method simulated moisture content profiles at the interface and was relatively smooth and continuous, thus the deviation from the reference solution was relatively lesser for the proposed method than other methods.

Figure 10

Calculated (a) RMSE, (b) MAE, and (c) EMC values after 8 hours for the test case 4 simulations.

Figure 10

Calculated (a) RMSE, (b) MAE, and (c) EMC values after 8 hours for the test case 4 simulations.

CONCLUSIONS

In the present study an equivalent effective saturation weighting in the SHC function-based interblock hydraulic conductivity has been incorporated in the improved modified Picard iteration scheme for solving Richards equation in the variably saturated soil system. The proposed weighting scheme was demonstrated in terms of the van Genuchten–Mualem SHC model in the present study. The proposed method was tested for its performance in terms of numerical accuracy with four test case examples which include different soil textures, different boundary conditions, different initial conditions, and different grid discretizations.

In general, it can be summarized that the numerical accuracy for all the developed models in this study deviated more and more from the reference solution as the grid discretizations increased. The proposed effective saturation weighting method improved the accuracy of the numerical results over arithmetic and higher-order mean methods relatively in all four unique test case simulations. This was mainly because of the effective saturation-based weighting in the parameteric SHC functions as adopted in the proposed method unlike other compared methods where weighting was mostly based on hydraulic conductivity of the subsequent nodes. The proposed method was considerably effective in improving the numerical solutions under steep boundary condition simulation scenarios such as infiltration in dry soils. Similarly, the proposed method was found effective in simulating water flows in fine-textured soils with relatively lower saturated hydraulic conductivity values.

The numerical simulations in this study were conducted based on the literature-based parameter datasets for the selected soils to demonstrate the performance of the proposed method against other methods. However, the proposed method performance will have to be verified with the field-based experimental results and corresponding customized SWR and SHC model parameters under specific real-world conditions. The numerical results simulated by the proposed method were compared against the recently published methods on interblock hydraulic conductivities. However, the present results were not directly compared with the Darcian-based methods. The present study demonstrates all the test cases in 1D vertical soil-water flow conditions. The impact of the proposed method on two dimensional (2D) and three dimensional (3D) soil-water flow modeling would be necessary for claiming the effectiveness of the proposed approach in 2D and 3D domains respectively.

ACKNOWLEDGEMENTS

We would like to acknowledge and thank the departmental computer facility (DCF), Department of Civil Engineering, Indian Institute of Technology Madras, India and Water Engineering and Management, School of Engineering and Technology, Asian Institute of Technology, Thailand for providing computing facilities to carry out the numerical simulations for this study.

REFERENCES

REFERENCES
Broda
S.
Paniconi
C.
Larocque
M.
2011
Numerical investigation of leakage in sloping aquifers
.
J. Hydrol.
409
(
1–2
),
49
61
.
Brooks
R. H.
Corey
A. T.
1964
Hydraulic Properties of Porous Media. Hydrology Papers, No. 3
.
Colorado State University
,
Fort Collins
.
Carsel
R. F.
Parrish
R. S.
1988
Developing joint probability distributions of soil water retention characteristics
.
Water Resour. Res.
24
,
755
769
.
https://doi.org/10.1029/WR024i005p00755
.
Celia
M. A.
Bouloutas
E. T.
Zarba
R. L.
1990
A general mass-conservative numerical solution for the unsaturated flow equation
.
Water Resour. Res.
26
(
7
),
1486
1496
.
Chen
J.-M.
Tan
Y.-C.
Chen
C.-H.
Parlange
J.-Y.
2001
Analytical solutions for linearized Richards equation with arbitrary time-dependent surface fluxes
.
Water Resour. Res.
37
(
4
),
1091
1093
.
Essig
E. T.
Corradini
C.
Morbidelli
R.
Govindaraju
R. S.
2009
Infiltration and deep flow over sloping surfaces: comparison of numerical and experimental results
.
J. Hydrol.
374
(
1–2
),
30
42
.
Farthing
M. W.
Ogden
F. L.
2017
Numerical solution of Richards’ equation: a review of advances and challenges
.
Soil Sci. Soc. Am. J.
81
,
1257
1269
.
Gastó
J. M.
Grifoll
J.
Cohen
Y.
2002
Estimation of internodal permeabilities for numerical simulation of unsaturated flows
.
Water Resour. Res.
38
(
12
),
1
10
.
Haverkamp
R.
Vauclin
M.
Touma
J.
Wierenga
P. J.
Vachaud
G.
1977
A comparison of numerical simulation models for one-dimensional infiltration
.
Soil Sci. Soc. Am. J.
41
,
285
294
.
Herrada
M. A.
Gutiérrez-Martin
A.
Montanero
J. M.
2014
Modeling infiltration rates in a saturated/unsaturated soil under the free draining condition
.
J. Hydrol.
515
,
10
15
.
Hillel
D.
1980
Fundamentals of Soil Physics
.
Academic
,
New York
, p.
413
.
Jiménez-Martínez
J.
Skaggs
T. H.
van Genuchten
M. T.
Candela
L.
2009
A root zone modelling approach to estimating groundwater recharge from irrigated areas
.
J. Hydrol.
367
(
1–2
),
138
149
.
Keese
K. E.
Scanlon
B. R.
Reedy
R. C.
2005
Assessing controls on diffuse groundwater recharge using unsaturated flow modeling
.
Water Resour. Res.
41
(
6
),
W06010
.
Kolditz
O.
Delfs
J.-O.
Bürger
C.
Beinhorn
M.
Park
C.-H.
2008
Numerical analysis of coupled hydrosystems based on an object-oriented compartment approach
.
J. Hydroinf.
10
,
227
244
.
Mohanasundaram
S.
Suresh Kumar
G.
Narasimhan
B.
2013
Numerical modelling of fluid flow through unsaturated zone using a dual-porosity approach
.
ISH J. Hydraul. Eng.
19
(
2
),
97
110
.
Oldenburg
C. M.
Pruess
K.
1993
On numerical modeling of capillary barriers
.
Water Resour. Res.
29
(
4
),
1045
1056
.
Richards
L. A.
1931
Capillary conduction of liquids through porous mediums
.
Physics (College. Park. Md)
1
,
318
333
.
Romano
N.
Brunone
B.
Santini
A.
1998
Numerical analysis of one-dimensional unsaturated flow in layered soils
.
Adv. Water Resour.
21
(
4
),
315
324
.
Šimůnek
J.
Sejna
M.
Saito
H.
Sakai
M.
van Genuchten
M. T.
2009
The HYDRUS-1D Software Package for Simulating the one-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media
.
Version 4.08. HYDRUS Softw. Ser. 3. Dep. Environ. Sci., Univ. Calif., Riverside. 332
.
van Genuchten
M. T.
Leij
F. J.
Yates
S. R.
1991
RETC Code for Quantifying the Hydraulic Functions of Unsaturated Soils
.
USDA
,
Riverside, California
.
Wang
L.
Huang
C.
Huang
L.
2018
Parameter estimation of the soil water retention curve model with Jaya algorithm
.
Comput. Electron. Agric.
151
,
349
353
.
Zlotnik
V. A.
Wang
T.
Nieber
J. L.
Simunek
J.
2007
Verification of numerical solutions of the Richards equation using a traveling wave solution
.
Adv. Water Resour.
30
(
9
),
1973
1980
.
Zymkiewicz
A.
2013
Modelling Water Flow in Unsaturated Porous Media. GeoPlanet: Earth and Planetary Sciences
.
Springer
,
Berlin, Heidelberg
.