Modelling Unsaturated Flow in Porous Media Using an Improved Picard Iteration Scheme


 The numerical solution of various systems of linear equations describing fluid infiltration uses the Picard iteration (PI). However, because many such systems are ill-conditioned, the solution process often has a poor convergence rate, making it very time-consuming. In this study, a control volume method based on non-uniform nodes is used to discretize the Richards equation, and adaptive relaxation is combined with a multistep preconditioner to improve the convergence rate of PI. The resulting adaptive relaxed PI with multistep preconditioner (MP(m)-ARPI) is used to simulate unsaturated flow in porous media. Three examples are used to verify the proposed schemes. The results show that MP(m)-ARPI can effectively reduce the condition number of the coefficient matrix for the system of linear equations. Compared with conventional PI, MP(m)-ARPI achieves faster convergence, higher computational efficiency, and enhanced robustness. These results demonstrate that improved scheme is an excellent prospect for simulating unsaturated flow in porous media.


Introduction
Many geotechnical engineering problems involve unsaturated seepage, such as the soil slope stability and dam infiltration involved in rainfall or groundwater changes, and the transportation of pollutants in landfill (Morway et al., 2013;Ku et al., 2017;Li et al., 2020). Therefore, effectively simulating and analyzing unsaturated seepage is of great practical significance (Wang et al., 2011;Zang et al., 2016;Wu et al., 2016;2020a).
The Richards equation (RE) provides a basic description of unsaturated infiltration (Richards, 1931), although analytical solutions to the RE are difficult to derive. Srivastava and Yeh (1991) used a Laplace transformation technique to obtain analytical solutions of the transient flows in homogeneous and two-layer soils, while Tracy (2006) proposed two-and three-dimensional analytical solutions of the RE using simple exponential models. Although their actual applicability is limited, these are good references for theoretical analysis (Zambra et al., 2012;Liu et al., 2015;Baron et al., 2017).
(2018) proposed a generalized finite difference method and an adaptive step-size Crank-Nicolson method to obtain more accurate numerical solutions to the RE. Dolejší et al. (2019) proposed an adaptive higher-order space-time discontinuous-Galerkin (hp-STDG) method to obtain highly efficient and accurate solutions to the RE. Svyatskiy and Lipnikov (2017) proposed a second-order-accurate finite volume scheme that solves the RE using high-order upwind algorithms for the relative permeability.
Additionally, other advanced numerical methods for solving variable saturated flow problems exhibit high levels of accuracy, computational efficiency, or ease of implementation under certain conditions. These include the two-level adaptive domain decomposition (Kuraz et al., 2015), finite analytic method (Zhang et al., 2016), and Chebyshev's spectral method (Wu et al., 2020b). A system of linear equations derived from the RE need to be solved iteratively, such as by using the Picard iteration (PI) scheme. Because hydraulic conductivity and water content are nonlinear functions of the pressure head, the numerical results to the RE using PI are relatively unreliable. A modified PI could ensure the mass conservation in the temporal discretization (Celia et al., 1990). Recently, Zha et al. (2017) further improved the algorithms of Celia et al. (1990) to overcome the numerical divergence that occurs when simulating infiltration into extremely dry soil. However, these improvements are limited to the overall mass balance of the RE, and do not improve the convergence rate of the solutions.
As PI has poor convergence and is time-consuming, increasing attention has focused on improving its convergence rate (Durbin et al., 2007;Lott et al., 2012;List and Radu, 2016;Illiano et al., 2020). Introducing a relaxation process can improve the convergence rate of PI (Paniconi and Putti, 1994), but determining the optimal relaxation factor is difficult. Durbin et al. (2007) proposed an adaptive under-relaxation of PI to accelerate the convergence rate for nonlinear groundwater flow problems. Later, De Smedt et al. (2010) developed an unstable manifold correction method for the adaptive relaxation of PI to solve the velocity field in higher-order ice-flow models.
Although the convergence rate of PI can be improved by adaptive relaxation techniques, the numerical results are not very robust. Preprocessing can improve the ill-conditioned nature of the associated systems of linear equations; that is, reducing the condition number of the coefficient matrix results in faster convergence (Benzi, 2002). Recently, Wang and Zhang (2003)

Numerical Procedures
The RE describes multidimensional unsaturated flow (Wu et al., 2013;Duc et al., 2020). The 2D RE is written as: where z and x are the vertical and horizontal coordinates, respectively; h represents the pressure head;  denotes the moisture content; () where  and () Kh can be described using exponential models (Gardner, 1958) where s K is the saturated hydraulic conductivity; s  and r  represent the saturated and residual moisture content, respectively; and  is the fitting parameter.
Furthermore, the control volume method (Patankar, 1980) can be used to discretize Eqs.
(1) and (2). A 1D infiltration model is illustrated in Fig. 1. Firstly, the nonuniform grid coordinates can be expressed as: where i denotes the nonuniform grid nodes; L is the height of the soil layer; and N represents the number of nodes. To ensure the mass balance in the temporal discretization, the following modified iteration format for / t   is proposed (Celia et al., 1990): where C is the specific moisture capacity, defined as ( )  denotes the time step; and the superscripts j and k represent the time step and number of iterations, respectively. Equation (6) has been widely used in commercial software such as HYDRUS (Šimůnek et al., 2009;Zha et al., 2017).
The discretized equation is now derived by integrating Eq.
(2) over the control volume in Fig. 1 and over the time interval from t to tt  (Patankar, 1980). Thus: Substituting Eq. (6) into Eq. (7), we obtain: where: For steady-state seepage, Eq. (8) can be simplified as: According to Eqs. (8) and (10), the matrix format of PI based on a nonuniform grid can be written as: where A is a symmetrical tridiagonal coefficient matrix; h and f are vectors. A system of linear algebraic equations is first derived at each iteration, and then a basic linear iterative method (such as Gaussian elimination, conjugate gradient method, or the Gauss-Seidel method) is used to solve the linear equations (Šimůnek et al., 2009).
As PI determines the computational efficiency, we use the Gauss-Seidel method to solve the system of linear equations at each iteration.
After solving Eqs. (11) and (12) for the first time, the coefficients in these equations are re-evaluated using the first solution, and the new system of linear equations is again solved. The iteration procedure terminates when the iteration error in terms of the l  norm satisfies the following: where  is the tolerance (set to 10 −8 in this paper).

Improved PI Scheme
PI has a first-order convergence rate with poor calculation efficiency. To improve the convergence rate of PI, adaptive relaxation and multistep preconditioner strategies are introduced. First, when k>1 in Eqs. (11) and (12), the following adaptive relaxation method is adopted: where k  is the kth relaxation factor, which is adjusted according to the generalized angle  between the current iteration increment h and the previous iteration increment h as follows:  reflects the trend of convergence in the numerical solution. An acute angle indicates better convergence, and an obtuse angle indicates oscillation. Therefore, when  is acute, k  should be increased appropriately; when  is obtuse, k  should be decreased appropriately. After many preliminary calculations, accurate results are calculated as follows: For each time step, the adaptive relaxed PI (ARPI) can be summarized as follows: (1) Set an initial h0 and tolerance  , and let k=1.
(3) Calculate the kth relaxation factor k  : If k=1, then k  =1; otherwise, the previous iteration increment h and relaxation factor 1 k   are read, and k  is calculated according to Eqs. (15) and (16).
(4) The calculation result of the current iteration step is corrected using Eq. (14).
where D represents the diagonal matrix of A; L and U represent the lower and upper triangular matrices of A, respectively. The adaptive relaxed PI scheme with multistep preconditioner (MP(m)-ARPI) proceeds according to the following algorithm.

Algorithm MP(m)-ARPI.
Given the number of steps 0 m  , the coefficient matrix A, and the filter parameter  : Compute an SGS preconditioner Solve the preconditioned system of linear equations using ARPI.
In this study,  is set to 10 −4 . It is obvious that, when m=0, MP(0)-ARPI completely degenerates to the ARPI scheme. A brief flowchart of MP(m)-ARPI for simulating unsaturated flow is shown in Fig. 2. To verify the calculation efficiency, the speed-up ratio is defined as: where PI T and

Example 1 (homogeneous soil)
To verify the effectiveness of the proposed method, we present an example in which a 1D steady-state, transient unsaturated flow is simulated in homogeneous porous media. The mathematical model is shown in Fig. 1 The initial condition is assumed to be ( ) 0 hz  . The boundary conditions are set as:    (Fig. 7).
The test results indicate that the proposed MP(m)-ARPI effectively simulates 1D steady-state and transient unsaturated flow in homogeneous porous media. The results illustrate that MP(m)-ARPI should use m=2 or 3 to ensure computational efficiency.
When the system of linear equations is more ill-conditioned, MP(3)-ARPI is optimal.

Example 2 (two-layer soils)
In two-layer soils, the coefficient matrix is extremely ill-conditioned because of the different physical parameters between soil layers (Liu et al., 2015;Zhu et al., 2019). To further verify the effectiveness of the proposed scheme, this example simulates the 1D unsaturated flow in a two-layer soil. The mathematical model is shown in Fig. 8, where L1=L2=5 m. The number of nodes N in each layer of soil is set to 80. The parameters and boundary conditions are consistent with those in Example 1. Table 2 presents various parameters for this example, with the hydraulic conductivity for layer 1 set to 10 -1 m/s and that for layer 2 ranging from 10 -2 -10 -9 m/s. Figure 9 shows that the condition number of the coefficient matrix increases as the hydraulic conductivity of the two layers diverges, but the condition number can be greatly reduced using the MP algorithm. In Case 8, in particular, using MP(3) reduces the condition number by an order of 10 7 . The computed results for Cases 2 and 6 are depicted in Figs. 10(a) and 10(b). From Fig. 11(a), it can be seen that the proposed MP(3)-ARPI generally requires the fewest iterations. Additionally, the speed-up ratios of MP(3)-ARPI relative to PI and ARPI are significant ( Fig. 11(b)).
Similar to Example 1, the 1D transient unsaturated flow was examined in two-layer soils. The parameters of the two-layer soils are listed in Table 3. MP(3)-ARPI was used to solve this example, with a total simulation time of 25 h a time step of 0.05 h, and the initial condition   3 , 0 10 h z t    . The computed results are shown in Fig. 12.

Example 3 (2D transient flow)
In this example, we simulate 2D transient flow in homogeneous porous media. In Fig.   13, the length (L) and width (W) of the model are both set to 1 m, and the boundary conditions are constant heads. The top boundary conditions can be expressed as follows: The initial condition is   ana h x z t can be obtained as follows (Tracy, 2006): The total simulation time and time step were set to 1 h and 0.05 h, respectively. Figure   14 shows that the condition number of the coefficient matrix increases as the mesh density increases when solving the 2D transient flow problem. Similarly, the condition number can be reduced significantly by using the MP algorithm. In Fig. 15(a), the average number of iterations required by MP(3)-ARPI is the lowest under different mesh densities. The calculation efficiency of MP(3)-ARPI is higher than that of MP(2)-ARPI when the mesh density is 30×30 (Fig. 15(b)). This numerical result is similar to that of Example 1. Figure 16 shows the numerical results obtained using MP(3)-ARPI.
Compared with the analytical solution (Fig. 17), the accuracy of the numerical solution is of order 10 -4 .

Summary and Conclusions
PI is often used to solve the RE and is easier to implement than Newtonian iteration, as there is no need to compute any derivatives. However, PI has a poor convergence rate and is time-consuming. To improve the convergence rate of PI, we have presented an improved PI scheme that combines adaptive relaxation and multistep preconditioner strategies. The results demonstrate the effectiveness of the improved PI scheme. The following conclusions can be drawn: (1) The proposed scheme effectively simulates unsaturated flow in porous media.
The results show that the proposed MP(m)-ARPI can effectively reduce the condition number of the coefficient matrix derived from the RE. Additionally, compared with conventional PI, CPI, and ARPI, the proposed MP(m)-ARPI offers faster convergence, higher computational efficiency, and enhanced robustness.

No
The solution vector is obtained at time step j.

Yes If No
End Yes Calculate the kth relaxation factor : If k = 1, then = 1; Otherwise, the previous iteration increment and relaxation factor are read, and is calculated according to Eqs. (15) -(16).
The preconditioned system of linear equations is solved, and the current iteration increment is calculated.
The numerical solution is corrected using Eq. (14).          Figure 1