Effective saturation-based weighting for interblock hydraulic conductivity in unsaturated zone soil water flow modelling using one-dimensional vertical finite-difference model

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.

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.

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.

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.

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)).

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 (S_i+1/2, S_i−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.

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

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.

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.

Test case 1 simulated a vertical 1D flow in a homogenous soil profile under constant pressure head boundary conditions with the value of h_z=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 h_z,t=0 = –500 cm at time t = 0. The bottom boundary condition was assigned with a free drainage type boundary condition (h_z=L,t = h_z=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.

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.

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.

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 h_z=0,z=L,t = –1 cm at both the top and bottom nodes for all time steps. The initial condition was assigned as h_z,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 (h_z,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 (h_z,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).

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.

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 h_z=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 (h_z=0 = –15,000 cm) pressure head value, the potential evaporation rate was reduced to a limited evaporation rate. The initial condition was assigned as h_z,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.

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 cm³/cm³.

The three vertical discretizations simulation results with respect to the reference solution were computed in terms of RMSE, MAE and E_MC 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 E_MC 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 (E_MC) 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 E_MC. Although the calculated RMSE, MAE and E_MC 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).

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 h_z,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 h_z,=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 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 E_MC 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.

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.

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.