{'fraction_non_alphanumeric': 0.041622334766140356, 'fraction_numerical': 0.02493935396008001, 'mean_word_length': 4.531544256120528, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Localized RBF methods for modeling infiltration using the Kirchhoff-transformed Richards equation

We develop a new approach to solve the nonlinear Richards equation based on the Kirchhoff transformation and localized radial basis function (LRBF) techniques. Our aim is to reduce the nonlinearity of the governing equation and apply LRBF methods for modeling unsaturated flow through heterogeneous soils. In our methodology, we propose special techniques which deal with the heterogeneity of the medium in order to apply the Kirchhoff transformation where we used the Brooks and Corey model for the capillary pressure function and a power-law relation in saturation for the relative permeability function. The new approach allows us to avoid the technical issues encountered in the Kirchhoff transformation due to soil heterogeneity in order to reduce the nonlinearity of the model equation. The resulting Kirchhoff-transformed Richards equation is solved using LRBF methods which have advantages in terms of computational cost since they don't require mesh generation. Furthermore, these LRBF techniques lead to a system with a sparse matrix which allows us to avoid ill-conditioned issues. To validate the developed approach for predicting the dynamics of unsaturated flow in porous media, numerical experiments are performed in one, two, and three-dimensional soils. The numerical results demonstrate the efficiency and accuracy of the proposed techniques for modeling infiltration through heterogeneous soils.


Introduction
Understanding infiltration through soils is of great importance in the fields of agriculture, hydrology, and water resources and environmental management. The modeling of infiltration processes is time-consuming and there is a need in the development of efficient techniques for these processes in the case of heterogeneous soils. The Richards equation [1] describes the dynamic of unsaturated flow through porous medium which is due to the actions of gravity and capillarity. Richards' equation is highly nonlinear because of the largely nonlinear dependencies of both unsaturated hydraulic conductivity and capillary pressure on saturation [2,3,4]. The van Genuchten [4] and Brooks-Corey models [3] are often used for the capillary pressure function of unsaturated soils. In terms of numerical analysis, the Gardner model [2] is important since large class of analytical solutions are available [5,6,7,8]. However, this model has some limitations in practical applications for describing unsaturated flow in soils [9].
The design of efficient numerical techniques for solving the Richards equation is very challenging due to the highly non-linearity of the equation and the technical issues encountered in the numerical treatment of soil heterogeneity. Various classes of approaches have been developed for modeling infiltration in soils such as finite difference methods [10,11,12], finite element methods [13,14,15,16] and finite volume methods [17,18,19,20,21]. While many methods have been developed for modeling unsaturated flow in soils, there is still a need for more efficient techniques to deal with the soil heterogeneity and the nonlinearity of the medium hydraulic properties as functions of saturation [3,4]. Most available approaches used iterative methods, such as Newton and Picard schemes [11,22,23,24,25,26], to linearize the system to be solved. While these iterative algorithms produce accurate results, they are expensive in terms of computational cost and they may have convergence issues for some flow conditions because of the highly non-linearity of the Richards equation and soil heterogeneity [24,27,28,26].
Among the numerical techniques proposed to solve the Richards equation is the Kirchhoff transformation approach [10,29,14,27,30,31]. Several studies have shown the efficiency of this approach because it reduces the non-linearity of the Richards equation [32,27,33,34]. However, this approach is developed for homogeneous soils or particular soil heterogeneity and capillary pressure functions [10,35,27,30,31]. Most of available techniques using the Kirchhoff transformation are limited to Gardner model for the capillary pressure function [36,37,38,39,27,31,40]. For instance, Suk and Park [31] recently developed a new numerical method based on the Kirchhoff transformation and the Gardner model to solve the Richards equation for layered soils. In their approach, the authors used a truncated Taylor series expansion to the Kirchhoff head at the material interface.
We note that the Gardner model is limited for practical applications compared to van Genuchten and Brooks-Corey models which are suitable for the entire range of pressure head [31]. Incorporating both high non-linear models for the capillary pressure [3,4] and soil heterogeneity is still a challenge in solving the Richards equation using the Kirchhoff transformation [31].
This study is a follow-up of the paper [41] in which we developed LRBF techniques for solving Richards equation in homogeneous medium where the Gardner model is used for capillary pressure. The techniques used in [41] which allowed us to linearize the system can not be applied in the case of heterogeneous mediums with other models of capillary pressure such as the Brooks-Corey [3] model used in this study. Here, we develop special techniques which deal with soil heterogeneity in order to apply the Kirchhoff transformation for solving the Richards equation based on the Brooks and Corey model [3] for capillary pressure. In our approach, a power-law relation in saturation for the relative permeability function is used to avoid the technical issues encountered in the Kirchhoff transformation due to the non-linearity of capillary pressure function and the heterogeneity of soils [31].
In our approach, we based on RBF meshless techniques [42,43,44] to solve the resulting system. These techniques don't require mesh generation and are based only on a set of independent points, which makes them advantageous in terms of computational cost. Due to their simplicity to implement, they represent an attractive alternative to the classical methods as a solution method for partial differential equations. Note that there are two versions of RBF methods. The global method [45] and the local one [43]. Several studies have demonstrated the efficiency of the global method [45,45,46] however it suffers from two major drawbacks: the ill-conditioned matrix obtained after the discretization process and the problem of choosing the adequate shape parameter for some RBFs [47,44,48]. To overcome these issues, the local methods were suggested [43,44]. The LRBF methods have advantages in terms of operational memory and calculation procedures where only inversion of sparse matrix are required. On the other hand, the LRBF methods are efficient in solving high-dimensional problems with complex boundaries [47,44,48,41] and are less sensitive to the choice of the shape parameter of RBFs as shown in [43]. These localized meshless methods have been successfully applied to a large variety of problems to solve partial differential equations including the Richards equation [49,33,50,41,34]. In this study, LRBF methods are applied to solve the resulting Kirchhoff-transformed Richards equation for modeling infiltration through soils.
The paper is organized as follows. In Section 2, we introduce the developed numerical approach and the Kirchhoff-transformed Richards equation. The proposed numerical model based on the LRBF method is described in Section 3. In Section 4, numerical simulations are conducted to validate the developed approach for modeling infiltration through soils. Finally, some concluding remarks are provided in Section 5.

Richards' model
We consider the traditional Richards equation describing infiltration through soils [1]: where θ s [L 3 /L 3 ] is the saturated water content, θ r [L 3 /L 3 ] is the residual water content and the parameter φ [−] is given by φ = (θ s − θ r ). Equation (2.2) is highly non-linear due to the nonlinear dependencies of the capillary pressure and relative permeability functions on saturation. Empirical constitutive relationships have been developed for these functions using experiments [2,3,4] and are used in previous studies [11,51,52,41,53].

Capillary pressure and relative permeability functions
Here, we used the Brooks and Corey's model [3] for the capillary pressure function. The saturation is given by:  [54]. The capillary pressure function can be expressed as follows: where the Leverett J-function is given by J(S) = S −1/λ . We propose to use the power-law relation in saturation for the relative permeability k r (S) = S β which can be expressed using capillary pressure as follows: where β > 1.
Let's introduce a reference constanth [L] for the capillary pressure in Equation (2.2) in order to use a dimensionless form of the capillary pressure function. For h h d , we obtain the following expression where we set ω =h/h d : andh can be taken as: where V is the volume of the whole domain. Since ω depends only on space (ω = ω(x)), we obtain: and We substitute Equations (2.6), (2.8) and (2.9) into Equation (2.2), we obtain the new form of the Richards equation associated with the Brooks and Corey model for h h d : where we separate the terms which are dependent on the capillary pressure function and those which depend only on space due to soil heterogeneity. For h > h d , we obtain:

Kirchhoff transformation
In this study, we propose to use the transformation of Kirchhoff which allows us to reduce the nonlinearity of the model equation. The Kirchhoff integral transformation is defined as: where we used the variable σ = h/h. Based on Equation (2.5), the transformation (2.12) can be rewritten as follows: where, we assume that λβ > 1. From Equation (2.13), we can determine the pressure head h as a function of ϕ: . (2.14) In our approach, in the transformation of the model equation, the variation of the intrinsic permeability is assumed dominant in the effect of heterogeneity variability [54,53]. We have for h h d : 16) and the third term of Equation (2.10) becomes: The time derivative in the first term of Equation (2.10) can be rewritten under the same assumption h h d by: Similarly, for h > h d , Equation (2.11) can be written in terms of ϕ: For simplicity, we will use the following parameters: (2.25) Finally, we obtain the new form of the model equation using the Kirchhoff ϕ: where the non-linearity of the original Richards model has been reduced since only the terms E and F are nonlinear and χ depends only on spatial coordinates x due to the heterogeneity of soils.

Approximation methods
In this section, we describe the approximation methods used to solve Equation (2.26). Let ∆t > 0 a time step and t p = p∆t with p ≥ 0 denotes the time level. Temporal discretization of Equation (2.26) using the backward Euler method may be written as: where ϕ p+1 denotes the approximate solution of ϕ at t p+1 , E p+1 and F p+1 are the estimated values of E and F computed using h p+1 which is obtained by substituting ϕ p+1 in Equation (2.14). By applying the Picard iteration scheme to Equation (3.1), we obtain: with m identifies iteration level. The solution is assumed to be known both at time level p and at iteration level m.
where N i denotes the number of interior points and N b denotes the number of points on the boundary For simplicity, we take the following expressions: . By extending the 1D spatial approximation used in [11] to 3D case, we obtain: where the expressions χ i−1/2 are given by: Along the iR are the corresponding values at the right of ϕ i and χ i respectively. Similarly, ϕ iL are the corresponding values at the left of ϕ i and χ i respectively. We then obtain: On the other hand, we approximate the following spatial operators as follows [11]: i−1/2 , (3.9) We then obtain: For simplicity, we use the following expressions: The operator L m is linear for each iteration level m. In addition to initial and boundary conditions, Equation (3.2) may be rewritten as follows: x ∈ Ω. (3.14) ϕ 0 and ϕ Γ are functions associated with the initial and boundary conditions. For each time level p, the linear system (3.14) is solved at each iteration level of Picard m until the following inequality is satisfied at all collocation points: where T ol is the error tolerance.

Approach using local radial basis functions
In this study, we use the local RBF meshfree method [44] which has advantages in terms of operational memory and calculation procedures where only inversions of small size matrices are required. For any point x s ∈Ω, the k -d tree algorithm is used [55] where {α p+1,m+1 i } ns i=1 are constants to be determined and ψ is a RBF for which, in our approach, we chose the exponential function given by ψ(r) = exp(−(cr) 2 ), where r = x s − x i [s] denotes the distance between x s and x i [s] and c > 0 is the shape parameter. According to Equation (3.16), we obtain: where , Υ m is considered as the expansion of Υ m [s] by inserting zeros in the proper position. It follows that: Similarly, for x s ∈ ∂Ω, we apply the linear operator B: Bϕ p+1,m+1 [s] We obtain the following sparse linear system: (3.23) The LRBF approach leads to system of sparse equations (3.23) which allows as to reduce the size of the dense matrices and avoid ill-conditioned problems arising from the global approach [56,57,48,58]. The approximate solutions ϕ p+1,m+1 = {ϕ p+1,m+1 (x i )} N i=1 can be obtained by solving the system (3.23).

Initial and boundary conditions
At t = 0, we assume that h(x, 0) = h 0 (x) which implies in terms of the Kirchhoff variable that ϕ(x, 0) = ϕ 0 (x), where .
The boundary conditions are expressed in terms of the Kirchhoff variable. For Dirichlet conditions, we consider h = g b for z = 0 and z = L which implies that: where g b is given by: For Neumann conditions, we consider −K ∂h ∂x = 0 for x = 0 and x = l 1 which implies that −χ ∂ϕ ∂x = 0. In the same way, we assume that −χ ∂ϕ ∂y = 0 at the lateral boundaries y = 0 and y = l 2 . Therefore, the linear operator corresponding to the boundary conditions is given by:

Numerical experiments
To validate the proposed approach for modeling unsaturated flow through heterogeneous soils, we present numerical solutions of Richards  In Section (4.1), numerical tests are performed using the developed model to simulate flow in unsaturated homogeneous soils. Section (4.2) presents numerical simulations of infiltration through heterogeneous soils. To investigate the accuracy of the developed technique, the results of numerical tests are used to compute the RM SE and L 1 er errors based on the following formulas: where θ(x i ) represents the approximate solution for the water content and θ ref (x i ) represents a reference solution and N is the number of collocation points.

Infiltration in homogeneous soils
In this section, we perform numerical simulations using the numerical method to simulate unsaturated flow in homogeneous soils. We consider different soil samples with L = 1 m. The parameters of these soils are shown in Table 1. The parameters θ 0 and h 0 are the initial

Soil
Type water content and pressure head respectively. In this numerical test, we used N z = 1001 uniform nodes with n s = 3 number of neighboring points and a time step ∆t = 0.0001. The water content profiles are shown in Figure 1 for the considered types of soils. Table 2 illustrates the RM SE, L 1 er errors between the numerical solutions and the reference solutions which are obtained using 1D-Hydrus [59]. We obtain accurate results and the predictions are in good agreement with the reference solutions simulated using 1D-Hydrus.  We compute the evolution of the total mass of water I for each numerical solution:  Figure 2 shows the evolution of the total mass of water for the numerical solutions obtained using the proposed techniques and 1D-solution obtained using Hydrus. The results show the effectiveness of the LRBF method in terms of conservation of mass. In the following, we perform numerical simulations for 3D infiltration problem. We consider a block of soil having the dimensions l 1 = l 2 = 0.3 m and L = 1 m. We consider the same physical parameters as the previous test. The silty clay and clay loam soils are chosen for this numerical test. We set c = 0.6, n s = 7, N x = N y = 90, N z = 300 and ∆t = 0.0001. Figures 3 and 4 show the 3D evolution of saturation (left) for the selected soils. The results on the right side are the x-slices of saturation (x = 0, x = l 1 /4, x = l 1 /2, x = 3l 1 /4, x = l 1 ). The cross sectional average in the vertical direction of the total mass of water of 3D numerical solutions and the 1D-Hydrus reference solutions (l x = l y = 1) are shown in Figure 5. The results confirm the accuracy of the proposed numerical method for infiltration through three-dimensional porous medium.

Sand
Silty clay

Infiltration in heterogeneous soil
In the following sections, we perform numerical tests to study the robustness of the developed numerical method in modeling one-, two-and three-dimensional heterogeneous medium.

Infiltration in 1D-layered soils
In this numerical test, we perform simulations of infiltration using a column of soil (L = 25.5 cm) with three layers. The layered soil consists of a thin surface crust (0.5 cm), a tilled layer (10 cm) and a subsoil layer (15 cm). The hydraulic properties [18] for the layers soil are shown in Table 3. Numerical simulations are performed for two cases using
6 displays the evolution in time of the water content (left) and pressure head (right). The numerical solutions are in good agreement with the 1D-Hydrus simulations. Table 4 presents the RM SE and L 1 er errors between the numerical solutions and the results obtained using 1D-Hydrus for the soil water content. The results confirm the effectiveness of the developed method in terms of accuracy in modeling infiltration in layered soils.

Infiltration in 2D-layered soils
This numerical test is performed for unsaturated flow through 2D layered porous medium. We consider the physical parameters of soils given in Table 3 and l 1 = 5 cm. Figure 7 displays the time-evolution of saturation for h 0 = −1000 cm obtained using the proposed method. The results are obtained using N x = 100, N z = 1001, ∆t = 0.005, n s = 5 and c = 0.6.
We present in Figure 8 the cross sectional average in the vertical direction of the total mass of water of 2D numerical solutions and the 1D-Hydrus reference solutions (l x = 1) where we consider two cases using h 0 = −100 cm and −1000 cm. The results show a good correspondence between the numerical and reference solutions which confirms the accuracy of the proposed method.

Infiltration in 3D-layered soils
Here, we investigate the capability of the developed numerical model in predicting infiltration through three-dimensional layered porous medium. We consider the same hydraulic properties of soils as the previous test. We perform numerical simulations using c = 0.6, n s = 7, N x = N y = 100, N z = 501 and ∆t = 0.001. In Figure 9, we display the 3D evolution of saturation (left) for the considered soils for h 0 = −1000 cm. The results on the right side present the x-slices of saturation (x = 0 , x = l 1 /4, x = l 1 /2, x = 3l 1 /4, x = l 1 ). The cross sectional average in the vertical direction of the total mass of water of 3D numerical solutions and the 1D-Hydrus reference solutions (l x = l y = 1) are shown in Figure 10. The comparison between the results of the total mass shows the accuracy of the developed numerical model for infiltration in three-dimensional layered soils.

Infiltration in layered soil of L-shape form
Here, we used a complex geometry of the interface between soils compared to the previous test. We perform numerical simulations using infiltration problem through layered soil of L-shape [60]. As shown in Figure 11, the computational domain is partitioned into two subdomains which differ in their saturated hydraulic conductivity K s . In [60], the van Genuchten model [4] is used for capillary pressure with the parameter values K s = 0.3319 m/h, θ s = 0.368, θ r = 0.102, α = 3.35 m −1 , n = 2 and m = 0.5, where α = 1/h cap , n and m are empirical parameters satisfy m = 1 − 1/n. In our case, we used the Brooks-Corey model [3] where we approximate λ and h d based on the equivalence between van Genuchten and Brooks-Corey parameters proposed in [61]. The parameters h d and λ are given by [61]: where S x = 0.72 − 0.35 exp(−n 4 ). We used a homogeneous Neumann condition on the vertical sides of the domain (x = 0, 1 m) and a homogeneous Dirichlet on the top and bottom sides (z = 0, 1 m). The initial condition is h(x, z, 0) = −z. We set ε = 0.1,

Water content (θ)
Pressure head (h) Figure 6: Vertical profile of soil water content (left) and pressure head (right). n s = 5, N x = N z = 1000 and ∆t = 0.001. Figure 14 displays the evolution of saturation at times T = 12, 24, and 48 h. We observe good overall agreement between the results of our simulations and those presented in [60,53].

Infiltration in curvilinearly layered soil
In this last numerical test, we perform simulations of infiltration through curvilinearly layered soil. The computational domain is split into two subdomains separated by a curved interface (see Figure 13). The interface equation is given by: ξ(x) = l 2 (0.1 (1 − cos(πx/l 1 )) + 0.45) . (4.6) The hydraulic parameters of the soils [53] are shown in Table 5. As in the previous test, the values of the parameters h d and λ are approximated using Equations (4.4) and (4.5). We  Elevation use l 1 = l 2 = 1 m, ε = 0.8, n s = 5, N x = N z = 1000 and ∆t = 0.001. We consider the initial condition h(x, z, 0) = −z. Homogeneous Dirichlet boundary conditions are imposed on the top and bottom sides of the domain, while homogeneous Neumann conditions are enforced on the two vertical sides. We present in Figure 13 the time evolution of the water content. The results of our simulations are in good overall agreement at different times compared to the results of simulation presented in [18,53], which confirm the robustness of the developed numerical model in predicting infiltration in heterogeneous soils.

Conclusion
In this study, a new approach is developed for modeling unsaturated flow through porous media. The proposed techniques are based on the Kirchhoff transformation, the Brooks and Corey model for the capillary pressure function and a power-law relation for the relative permeability function. The proposed approach allows us to avoid technical issues associated with the use of the Kirchhoff transformation in heterogeneous soils and to reduce the nonlinearity of the model equation. The resulting system is solved based on the LRBF methods which are very effective for solving high-dimensional problems since they don't require mesh generation and have a computational advantage of using reduced memory. The LRBF meshless methods allow us to avoid ill-conditioning problems where a sparse matrix is obtained for the global system. The performance and robustness of the developed numerical model are demonstrated based on comparisons between numerical and reference solutions. Numerical experiments are performed to simulate the infiltration in one, two, and three-dimensional soils. The numerical results show the accuracy of the proposed techniques for modeling infiltration through porous media.