J.
Indones.
Math.
Soc.
Vol.
No.
, pp.
1Ae12.
Implementation of Dual Reciprocity Boundary Element Method to Solve The Mathematical Model of Steady Infiltration in Heterogeneous Soils Vertically Armando Deminto Paa1 .
Imam Solekhudin2O Department of Mathematics.
Universitas Gadjah Mada.
Indonesia, armandodemintopaa@mail.
id, 2 imams@ugm.
Abstract.
The Dual Reciprocity Boundary Element Method (DRBEM) is applied to RichardsAo equation for steady infiltration in unsaturated soil.
The types of soil observed are Guelph Loam (GL) and Pima Clay Loam (PCL).
A mathematical model for infiltration in heterogeneous unsaturated soil is built using RichardsAo equation and the Gardner model, with GL and PCL soil parameters.
The heterogeneous soil types include GL-PCL and PCL-GL, where each soil is assumed to vary linearly in the vertical direction.
This study investigates infiltration in vertically heterogeneous Using DRBEM, the numerical simulation results of hydraulic conductivity are obtained as an approximation of the constructed mathematical model.
Based on the numerical results of the hydraulic conductivity, the soil water potential values can be calculated using the Gardner model.
The results of this study provide the distribution of the water content under steady-state conditions in heterogeneous soils, determined by hydraulic conductivity and soil water potential.
Key words and phrases: RichardsAo equation, hydraulic conductivity, soil water
INTRODUCTION
Numerical methods have been used to solve a mathematical model of steady infiltration, which is formulated based on the RichardsAo equation.
Researchers have widely used the Finite Element Method (FEM).
FEM discretizes the entire domain into smaller finite elements, allowing flexible modeling of complex geometries.
However, this method often requires significant computational resources due to the need for meshing across the entire domain, leading to increased complexity and processing time.
In contrast, the Boundary Element Method (BEM) provides a O Corresponding author 2020 Mathematics Subject Classification: 65M38, 65N38, 74S15.
Received: 05-03-2025, accepted: 13-06-2025.
more efficient alternative by discretizing only the boundary.
Therefore, the number of computational elements is significantly reduced, leading to lower computational costs .
BEM is particularly advantageous for problems dominated by boundary conditions, such as infiltration, where it offers higher accuracy and efficiency.
The problems of steady infiltration in a heterogeneous soil are considered.
This study aims to analyze the hydraulic conductivity (K).
/da.
values, which will then yield the soil water potential (Sp ).
, 3, 4, .
of soil under the steady-state condition.
The soil used in this study is Guelph Loam (GL) and Pima Clay Loam (PCL).
The parameters required for this study are the soil coarseness () and the saturated hydraulic conductivity (Ks ) of GL and PCL.
According to .
, the values of these parameters are presented in Table 1 below.
Table 1.
SoilsAo Parameters in Homogeneous Soil Soils PCL .
Oe1 ) Ks .
/da.
On homogeneous soil, the parameters and Ks in GL and PCL can be used to calculate the values of K using DRBEM, which can then be obtained to find the values of Sp .
Using the values of K and Sp in homogeneous soils, the study of steady infiltration will be applied to a heterogeneous soil layer as a rectangular column with a width of 1 m and a depth of 5 m.
The infiltration process occurs when water enters the soil at a certain rate.
In this study, it is assumed that there is a constant flow on the soil surface, no flow on the vertical sides, and a water table at a depth of 5 m.
The flux on the soil surface is assumed to be a constant value q0 , given by q0 = 0.
75 y 0.
099 m/day .
099 m/day is the saturated hydraulic conductivity of PCL.
The value of q0 in Equation .
is the same as in .
The following figure 1 illustrates the soil column that has been described.
Consider two types of heterogeneous soils.
GL-PCL and PCL-GL.
GL-PCL
denotes a soil column with GL on the surface and PCL at the bottom, specifically at a depth of 5 meters, while the soil in between is assumed to transition linearly from GL to PCL in the vertical direction.
Similarly.
PCL-GL denotes a soil column with PCL on the surface and GL at the bottom, specifically at a depth of 5 meters.
In contrast, the soil in between is assumed to transition linearly from PCL to GL in the vertical direction.
Therefore, the values of and Ks for heterogeneous soils are assumed to follow the form shown in Table 2 below.
Figure 1.
A Column of Soil Table 2.
SoilsAo Parameters in Heterogeneous Soil
Soils GL-PCL
PCL-GL
Oe1 ) Ks .
/da.
4 Oe 5 z 0.
3171 Oe 0.
4 25 z 0.
MAIN RESULTS
Mathematical Formulation.
RichardsAo equation .
, 10, 11, 12, .
is utilized to construct the mathematical model for infiltration in unsaturated heterogeneous soil in this study.
The equation is given by OC
OCSp
OCSp
OCK
OCt OCx OCx OCz OCz OCz where is the volumetric water content.
, 15, 16, .
Flux normal to surface with outward pointing normal n = .
x , nz ) is given by .
OCSp OCSp Oe 1 nz .
F = U n1 V n2 = OeK OCx OCz The hydraulic conductivity function proposed by Gardner .
is given by.
K = Ks eSp where and Ks are the soil parameters related to soil coarseness and saturated hydraulic conductivity, respectively.
The soil coarseness is assumed to be linearly changed along the z-axis formulated as = c dz, .
for any c, d OO R such that > 0.
By Equations .
we have OCSp
1 OCK
OCx
K OCx
OCSp
1 OCK
OCz K OCz Equation .
can be written as 1 OC2K
1 OC2K
OCK
Sp K.
OCt OCx2 OCz 2 OCz and the steady infiltration equation in heterogeneous soil is obtained as we set Equation .
to be OC2K
OC2K
OCK
Oe Sp K = 0.
OCx2 OCz 2 OCz Normal flux .
may be written as 1 OCK bSp K Knz OCn Since we assume that water flows on the soil surface and there is no flow on the vertical sides of the column of soil, we have F =Oe F = Oeq0 for 0 < x < 1 and z = 0, on the surface of soil.
F = 0 for x = 0 and 0 < z < 5.
F = 0 for x = 1 and 0 < z < 5.
Based on the water flux condition, the boundary conditions are obtained as follows OCK = q0 Oe K Oe dSp K, on the surface of soil.
OCn
OCK
= 0, for x = 0 and 0 < z < 5.
OCn
OCK
= 0, for x = 1 and 0 < z < 5.
OCn K = Ks , for 0 < x < 1 and z = 5.
Dual Reciprocity Procedure.
The steady infiltration equation .
can be solved numerically using DRBEM.
By applying the procedure of this method, we obtain Z OC.
, z.
, )
, .
c(, )K(, ) = K.
, .
Oe .
, z.
, )
OCn OCn OCK.
, .
, z.
, ) f1 .
, .
Oe f2 .
, .
, .
dx dz.
OCz Equation .
is the solution of the boundary integral equation of 9, where .
, z.
, )
is the fundamental solution of the two-dimensional Laplace equation, and (, ) lies on smooth part e c(, ) = 2 (, ) OO E.
A numerical solution can be obtained by discretizing the boundary e into line segments and choosing the number of interior points in the E.
Let C .
C .
, .
C (N )
be the line segments on e .
OO C .
O C .
O A A A O C (N ) ) and point .
, b.
) be the midpoint of C .
for i = 1, 2, .
N .
Furthermore, let .
(N .
, b(N .
), .
(N .
, b(N .
), .
(N L) , b(N L) ) be the interior points in E.
An illustration of a typical discretization and interior point is shown in Figure 2.
Figure 2.
A Typical Discretization Using the midpoints and interior points described in Figure 2, we can obtain the numerical solution through DRBEM by recasting the integral Equation .
as c.
K .
= .
[K .
F2 .
Oe KE .
F1 .
AAz Oe f2 AA.
]K .
, .
for m = 1, 2, .
N L.
The notations in .
are c.
= c.
, b.
K m = K.
, b.
KE .
= OCK OCn |.
,b.
) , .
, z.
, ) ds.
, .
, .
OC .
, z.
, ) ds.
, .
C .
OCn AA.
= .
O .
, and AA.
" N L AE.
NX O .
The symbol AE.
is (OCA/OCZ)|.
,b.
) , where A represents the radial basis function centered at the point .
, b.
For mathematical convenience, let be c and let Sp be the average value of the Sp of GL and PCL, denoted by SAp .
Then, .
f1 and f2 are given by .
dSAp c and 2d2 A .
Sp .
The notation .
is defined as .
=c.
, b.
, a.
, b.
) N OCN.
, z.
, b.
) .
, b.
, b.
)F2 OCn .
Oe .
,b.
) and N is a function that satisfying OC2N OC2N 2 = A.
OCx2 OCz By taking to be c and Sp to be Sp , the boundary conditions applied to this numerical approximation are OCK = cq0 Oe cK Oe dSAp K, on the surface of soil.
OCn
OCK
= 0 for x = 0 and 0 < z < 5.
OCn
OCK
= 0 for x = 1 and 0 < z < 5.
OCn K = Ks for 0 < x < 1 and z = 5.
Numerical Simulation.
This subsection presents the numerical solution of steady infiltration in the heterogeneous soil equation constructed in the preceding section using DRBEM.
The collocation points are selected for segments and the interior of the domain.
There are 132 points on the boundary and 2401 points on the interior.
The selection of these collocation points optimizes computational efficiency and facilitates smooth numerical interpolation to visualize the surface plot.
The numerical simulation in this study results in an estimation of K in heterogeneous soil.
Using the obtained values of K, we can determine the corresponding values of Sp .
We focus on analyzing the obtained numerical results of K and Sp in the domainAos interior.
We denote that K .
and Sp as the numerical results of K and Sp respectively, at .
the interior points.
By the obtained K .
Sp can be calculated by .
Sp = .
for p = 1, 2, .
, 2401.
Hydraulic Conductivity in Homogeneous and Heterogeneous Soil.
The numerical solution of K obtained in this study is shown in Figure 3, which represents the numerical values of K in GL and PCL, together with GL-PCL and PCL-GL.
The values of K on the soil surface for the four types of soil are relatively equal, as they are influenced by the constant water flow rate q0 given in Equation .
Overall, the values of K increase from a depth of 3.
4 m to 5 m, which is the surface of the water table.
The numerical simulation yields the values of K, which are relatively uniform in the horizontal direction.
Figure 4 illustrates the distribution of K for homogeneous and heterogeneous soils along the soil depth.
In GL, numerical simulations result in varying K values with every change in depth.
The values of K tend to decrease, starting from the shallow surface at a depth of 0.
1 m about 0.
07417086 down to a depth of 3.
3 m about 0.
This occurs because the constant flow rate q0 causes more water to be retained at the shallow surface.
Then, the values of K increase from a depth of 3.
4 m about 07096157 to a depth of 4.
9 m about 0.
2451368, approaching 0.
3171, the saturated hydraulic conductivity of GL.
The increase in K is influenced by the increasing water content as the soil approaches the water table at a depth of 5 m.
In PCL, the values of K tend to increase as depth increases.
At the shallow surface, at a depth of 0.
1 m, the value of K is about 0.
It rises to about 09579711 as depth reaches 4.
9 m and approaches 0.
099, which is the saturated hydraulic conductivity of PCL.
In GL-PCL, the numerical simulation results in varying K values at every change in depth.
The values of K tend to decrease, starting from the shallow surface at a depth of 0.
1 m about 0.
0739798 to a depth of 3.
3 m about 0.
Then, the values of K increase from a depth of 3.
4 m about 0.
07193379 to a depth of 9 m about 0.
09205779, approaching the saturated conductivity of PCL, which is Figure 3.
Surface Plot of K in Homogeneous and Heterogeneous Soils This result is influenced by GL-PCL characteristics, given by the parameter , which is 3.
4 Oe 52 z.
Similar to the distribution of K in GL and GL-PCL, numerical simulations in PCL-GL result in varying values of K.
The value of K decreases from a depth of 0.
m about 0.
0718888 to a depth of 2.
4 m about 0.
Then, there is a significant increase in the value of K starting from a depth of 2.
5 m, about 0.
06644589 to a depth of 4.
9 m, about 0.
27001055, approaching the saturated conductivity of GL, which is 0.
This result is influenced by PCL-GL characteristics, given by parameter , which is 1.
4 25 z.
Soil Water Potential in Homogeneous and Heterogeneous Soil.
The numerical solution of Sp obtained in this study is shown in Figure 5, which represents the numerical values of Sp in GL and PCL, together with GLPCL and PCL-GL.
Figure 5 represents the Sp in homogeneous and heterogeneous soils during steady vertical infiltration.
The numerical simulation yields the values of Sp , which are relatively equal along the horizontal line.
Figure 6 shows the variability of the value Sp in both homogeneous and heterogeneous soils along the soil depth.
Figure 4.
Hydraulic Conductivity (K) Along The Depth of Soil Figure 6.
Soil Water Potential (Sp ) Along The Depth of Soil Figure 5.
Surface Plot of Sp in Homogeneous and Heterogeneous Soils In GL, the values of Sp vary with every change in depth.
Sp tends to increase .
n larger negative value.
starting from the shallow surface at a depth of 0.
about Oe0.
4273077 to a depth of 3.
3 m about Oe0.
This occurs because there is constant water flow, causing the water content at the shallow surface to increase as depth increases to 3.
3 m.
Subsequently, the values of Sp decrease from a depth of 3.
4 m about Oe0.
44032601 to a depth of 4.
9 m about Oe0.
This decrease occurs as depth increases towards the water table, decreasing Sp .
Unlike in GL.
PCL shows that its Sp values tend to decrease as depth increases.
The values of Sp decrease from a depth of 0.
1 m about Oe0.
20514084 to a depth of 4.
9 m about Oe0.
The Sp of PCL is lower than that of GL.
This is caused by the soil parameter of PCL, which is 1.
4, less than GL with is 3.
The soil type, clay, results in a weaker ability of PCL to retain water than GL.
Furthermore, in GL-PCL, the Sp results in varying values at every change in The Sp tends to increase, starting from the shallow surface at a depth of 0.
m about Oe0.
42903979 to a depth of 1.
9 m about Oe0.
This increase occurs because of a constant water flow, causing the water content at the shallow surface.
In contrast, the Sp tends to decrease from a depth of 2 m about Oe0.
44303507 to a depth of 4.
9 m about Oe0.
The increasing water content influences the decrease as depth increases toward the water table at 5 m.
Similarly to GL and GL-PCL, the calculation in PCL-GL results in varying Sp values at every change in depth.
The Sp tends to increase from a depth of 0.
m about Oe0.
2521646 to a depth of 2.
5 m about Oe0.
This increase occurs due to the constant water flow, causing the water content at the shallow surface to Conversely.
Sp tends to decrease from a depth of 2.
6 m about Oe0.
to a depth of 4.
9 m about Oe0.
The decrease of Sp is influenced by the increasing water content at the specified depths as depth increases towards the water table at a depth of 5 m.
CONCLUDING REMARKS
Steady infiltration in heterogeneous soils has been mathematically modeled and solved using the DRBEM.
The application of DRBEM to the constructed model has resulted in numerical values of hydraulic conductivity, from which the soil water potential values have been determined.
According to the results, the hydraulic conductivity values depend on soil characteristics at any depth level.
shallow soil, the obtained values of hydraulic conductivity are influenced by the fluxes on the surface, which will increase as depth increases to the water table since it contains more water.
The values of soil water potential in heterogeneous soil depend on its characteristics at any depth level.
The more water contained in the soil, the lower the soil water potential values.
Acknowledgement.
This research is partially supported by research grant Penelitian Fundamental Reguler (PFR) 2025, contract number 067/C3/DT.
00/PL/2025.
2412/UN1/DITLIT/Dit-lit/PT.
03/2025.
REFERENCES