Maria Lobo*
| Petrus D. L. Moensaku | Waltrudis A. Asa | Robertus Guntur
© 2025 The authors. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).
OPEN ACCESS
Efficient modeling of soil water infiltration is essential for sustainable furrow irrigation, particularly in semi-arid regions with limited water resources. This study aims to simulate transient infiltration in trapezoidal irrigation channels with layered Pima clay loam soil, where boundary heterogeneity and hydraulic complexity make numerical modeling challenging. The Dual Reciprocity Boundary Element Method (DRBEM) was applied to the Richards equation, which was transformed into a modified Helmholtz form and solved numerically using MATLAB. Simulations were performed at four soil depths beneath the furrow base (0.15, 0.25, 0.35, and 0.45 cm), where θ represents the water content at each depth, and these fixed values represent the spatial observation points of the of the matric flux potential $\phi_{\mathrm{m}}$. Results show that increasing boundary discretization improves numerical accuracy and stability up to an optimal range of N = 160 - 200, while a slight non-monotonic fluctuation at N = 320 is attributed to numerical over-refinement rather than a physical anomaly. The findings confirm that DRBEM accurately captures the hydraulic potential gradients in fine-textured soils and efficiently models unsaturated flow dynamics under semi-arid conditions. This study provides the first application of DRBEM for transient infiltration modeling in trapezoidal furrow irrigation with layered Pima clay loam soil, offering a robust and computationally efficient framework for improving irrigation management and water use efficiency.
DRBEM, furrow irrigation, numerical modeling, Pima clay loam soil, water infiltration
Efficient irrigation water management is a global priority for ensuring food security and agricultural sustainability, particularly in semi-arid and drought-prone regions where precipitation is limited, and evapotranspiration is high. Agriculture accounts for more than 70% of global freshwater withdrawals, yet inefficiencies in irrigation practices result in substantial water losses, groundwater depletion, and declining crop productivity [1-3]. As climatic variability intensifies, improving infiltration modeling and water distribution efficiency has become a cornerstone in addressing agricultural vulnerability and ensuring long-term sustainability [4].
At the national level, Indonesia faces comparable challenges in optimizing irrigation performance. Although surface irrigation remains the most widely used system due to its simplicity and low cost, it often suffers from substantial conveyance losses through infiltration and evaporation. The absence of real-time monitoring, coupled with limited technical capacity in rural areas, exacerbates inefficiencies in water allocation and soil moisture retention [5, 6]. These challenges are particularly evident in regions characterized by heterogeneous soils and inadequate canal maintenance, where infiltration behavior is difficult to predict and control.
At the regional scale, East Nusa Tenggara (NTT) exemplifies these issues. The province experiences a semi-arid climate with short rainy seasons, prolonged droughts, and stratified clay loam soils that complicate infiltration dynamics [7-9]. Farmers in NTT predominantly use trapezoidal furrow irrigation channels, which, while easy to construct, often lead to uneven infiltration, lateral seepage losses, and suboptimal water distribution. These limitations underscore the urgent need for site-specific, physically based infiltration models capable of capturing the complex interaction between soil layering, geometry, and boundary conditions under semi-arid environments.
Globally, numerous integrated irrigation management approaches—such as regulated deficit irrigation, rainwater harvesting, and soil-moisture-based automation—have significantly enhanced water-use efficiency in comparable semi-arid contexts [10-14]. However, the adoption of such precision technologies in NTT remains limited due to economic and infrastructural constraints. Consequently, numerical modeling represents a viable and cost-effective alternative to improve irrigation design, allowing researchers to simulate transient soil–water interactions and optimize infiltration performance [15, 16].
Infiltration, defined as the downward and lateral movement of water from the surface into the subsurface soil profile, has traditionally been represented through empirical models such as the Kostiakov and Philip equations [17, 18]. While computationally simple, these models assume homogeneous soil texture and steady-state flow, which limits their applicability in layered and anisotropic soils [19, 20]. To overcome these limitations, semi-analytical approaches—notably the modified Green-AMPT equation—introduce additional parameters for hydraulic conductivity contrasts between soil layers [21, 22]. However, these formulations remain inadequate under transient and spatially variable conditions typical of semi-arid environments [23, 24].
To address the shortcomings of empirical and semi-analytical formulations, numerical models have been increasingly employed to capture complex soil–water interactions. Among these, the Dual Reciprocity Boundary Element Method (DRBEM) offers distinct advantages by converting domain integrals into boundary-only computations using Radial Basis Functions (RBFs). This transformation substantially reduces computational cost while maintaining high numerical precision [25-29]. DRBEM has proven effective in solving nonlinear, time-dependent infiltration problems in heterogeneous soils, and when integrated with field soil-moisture data, it significantly improves infiltration prediction accuracy [4, 30].
Recent advances have expanded the DRBEM to simulate diverse geometries, soil textures, and boundary conditions, demonstrating its flexibility in hydrological modeling [31-35]. However, most studies remain idealized, rarely addressing trapezoidal channels, layered anisotropic soils, or transient boundary effects typical of semi-arid irrigation. Comparative assessments of DRBEM with conventional numerical methods are also limited [36].
This gap underscores the need for models capable of simulating transient infiltration from trapezoidal surface channels into stratified soils, where lateral seepage and heterogeneous conductivity significantly affect water movement. This study proposes a novel DRBEM-based transient infiltration model that integrates trapezoidal geometry, layered Pima clay loam soil, and semi-arid conditions. To the authors’ knowledge, it represents the first application of DRBEM in this context, offering both methodological and practical contributions to irrigation modeling.
The model is developed from Richards’ equation, transformed into the modified Helmholtz equation, and used to examine the effects of boundary discretization and soil depth on accuracy and computational efficiency. By leveraging DRBEM’s boundary-only formulation, this study provides a robust and efficient framework for simulating transient infiltration under realistic semi-arid conditions.
2.1 Comparative analysis of modeling approaches
To contextualize the current study, Table 1 compares DRBEM with conventional numerical and empirical methods regarding computational cost, accuracy, and flexibility.
Table 1. Comparison of infiltration modeling approaches
|
Method |
Key Features |
Computational Cost |
Accuracy |
Flexibility (Geometry/Boundary) |
Main Limitations |
References |
|
Empirical Models (Kostiakov, Philip) |
Simple equations derived from observed data |
Very low |
Moderate (steady-state only) |
Low – assumes homogeneity |
Cannot handle transient or layered conditions |
[17, 18] |
|
Finite Difference Method (FDM) |
Domain discretization using grid-based schemes |
High |
High (for simple geometries) |
Low – limited to rectangular domains |
Inefficient for irregular or open boundaries |
[37-39] |
|
Finite Element Method (FEM) |
Domain discretization with triangular or quadrilateral meshes |
High |
Very high |
Moderate – flexible for complex geometries |
Requires extensive meshing and high computational effort |
[36] |
|
Boundary Element Method (BEM) |
Boundary-only discretization; integral formulation |
Moderate |
High |
High for simple boundaries |
Limited to problems with known fundamental solutions |
[40-43] |
|
DRBEM |
Converts domain integrals using RBFs |
Low–moderate |
Very high |
Very high – handles complex and transient boundaries efficiently |
Still underutilized for real-field soil infiltration modeling |
[25, 26, 31] |
2.2 BEM
The BEM is a numerical technique for solving problems involving partial differential equations (PDEs). The method works by transforming the fundamental solution of Laplace’s equation [44, 45] into boundary integral equations.
BEM converts the governing differential equation into an integral equation, and the method involves discretizing the boundaries into interconnected segments to approximate the boundary integral equations [46]. This method can be applied to a variety of equations, including Helmholtz’s Equation, Diffusion-Convection Equation, Laplace’s Equation, Potential and Viscous Flow Equations, and Electrostatic and Electrodynamic Linear. The boundary integral equation is then discretized into several elements, either in the form of cells or segments. Since BEM discretizes the boundary directly from the governing differential equation, it requires fewer elements than other methods. As a result, the linear algebraic system’s matrix is smaller, enabling faster computational solutions [18].
2.3 DRBEM
The DRBEM is an extension of the BEM. While the Helmholtz equation requires a solution, its fundamental solution is difficult to obtain and is non-singular. To address this, DRBEM was developed to solve PDEs for which fundamental solutions are challenging to find. DRBEM modifies BEM to handle such cases effectively, providing a method to solve these complex equations.
2.4 Boundary integral equation
The Laplace solution of Laplace’s equation is
$\Phi(x, y ; \xi, \eta)=\frac{1}{4 \pi} \ln \left[(x, \xi)^2 \mid(y, \eta)^2\right]$ (1)
Domain modification is performed for $(\xi, \eta) \in R$ and $(\xi, \eta) \in C$, which is analogous to the boundary element method yielding the equation.
$\begin{gathered}\lambda(\xi, \eta) \phi(\xi, \eta)=\int_C(\phi(x, y)) \frac{\partial \Phi(x, y ; \xi, \eta)}{\partial n}-\phi(x, y ; \xi, \eta) \frac{\partial \phi(x, y)}{\partial n} d S+\iint_R \Phi(x, y ; \xi, \eta)(g(x, y))-k^2 \phi(x, y) d x d y\end{gathered}$ (2)
where,
$\lambda(\xi, \eta)=\left\{\begin{array}{c}0, \text { if }(\xi, \eta) \in R \cup C \\ \frac{1}{2}, \text { if }(\xi, \eta) \text { on } C \\ 1, \text { if }(\xi, \eta) \text { on } R\end{array}\right.$
Eq. (2) represents a boundary integral equation for the Helmholtz Equation.
2.4 Domain integral approach
In Eq. (2), the domain integral over R is still present. To address this, the domain integral is approximated using a linear combination of RBFs specifically g(x, y-k 2 \phi(x, y) is approximated by the sum RBFs centered at collocation points:
$g(x, y)-k^2 \phi(x, y) \approx \sum_{m=1}^M \beta^{(m)} \rho\left(x, y ; a^{(b)}, b^{(m)}\right)$ (3)
$\begin{aligned} \beta^{(m)}\left(\int_C \phi(x, y)\right. & \left.\frac{\partial \Phi(x, y ; \xi, \eta)}{\partial n}-\Phi(x, y ; \xi, \eta) \frac{\partial \phi(x, y)}{\partial n}\right) d s \approx \sum_{k=1}^N \int_{C^{(k)}}\left(\bar{\phi}^{(k)} \frac{\partial \Phi}{\partial n}(x, y ; \xi, \eta)-\Phi(x, y ; \xi, \eta) \bar{p}^{(k)}\right) d s =\sum_{k=1}^N\left[\bar{\phi}^{(k)} f_2^{(k)}(\xi, \eta)-\bar{p}^{(k)} f_1^{(k)}(\xi, \eta)\right]\end{aligned}$ (4)
This expression in (4) is used to compute the $\phi(a, b),(a, b) \in C \cup R$. Thus, we have:
$\lambda(a, b) \phi(a, b)=\sum_{j=1}^{N+L} \mu^{(j)}\left[g\left(\bar{x}^{(j)}, \bar{y}^{(j)}\right)-k^2 \bar{\phi}^{(j)}\right]+\sum_{k=1}^N g g\left[\bar{\phi}^{(k)} f_2^{(k)}(a, b)-\bar{p}^{(k)} f_1^{(k)}(a, b)\right]$ (5)
In this section, we will discuss the results obtained for the problem of infiltration of furrow irrigation water with the reciprocity boundary element method for the Pima clay loam soil type.
This study utilizes a numerical simulation approach to address the complex problem of transient water infiltration in unsaturated soils. The methodology revolves around solving Richards’ equation, which is widely accepted as the governing model for unsaturated soil water movement. However, due to its strong non-linearity, simplification is often required for practical computation. Recent research shows that linearization techniques, when paired with advanced numerical methods such as the DRBEM, can enhance computational efficiency and maintain solution accuracy under transient conditions [31]. These improvements allow for better modeling of the dynamic interactions between soil texture, moisture, and infiltration rate in semi-arid irrigation systems. Surface irrigation is one of the oldest and most commonly used water application techniques, where water flows over the land surface, with some infiltrating into the soil and the rest either running off or being stored.
Taken together, these studies provide a strong foundation for applying DRBEM in solving complex infiltration problems, especially where traditional analytical or empirical methods fall short. However, there remains a need for context-specific validation and refinement of these models in semi-arid environments like NTT, where real-world variables—such as clay-loam stratification and irregular flow—require more robust, efficient, and adaptive numerical solutions. If embankments are constructed, the water will infiltrate the soil. The furrow irrigation system is designed to meet agricultural needs in regions with low or inconsistent rainfall, ensuring that rows of crops receive sufficient water. On dry soil, pressure is used to draw water from small pores. Suction potential or capillary potential, denoted as $\psi$, is closely related to suction pressure, $p$, where $p=\rho g \psi$. In this equation, $\rho$ represents the density of water, $g$ is the gravitational acceleration, and $\psi$ is negative, indicating the suction pressure. When the soil is saturated, $\psi=0$. Through capillarity action, plants can collect water up to a height, $H$ [47]. The mathematical formulation adopted in this study aligns with theoretical assumptions regarding exponential functions of hydraulic conductivity in relation to soil water potential, as previously detailed [48].
This study implements the DRBEM to solve the modified Helmholtz equation derived from the linearized Richards’ equation. DRBEM transforms domain integrals into boundary integrals using RBFs, allowing the governing equations to be solved through boundary-only discretization.
The general DRBEM formulation is expressed as:
$\nabla^2 \phi-\lambda^2 \phi=f(x, z)$,
where, $f(x, z)$ represents the non-homogeneous term approximated using a series expansion of RBFs:
$f(x, z) \approx \sum_{m=1}^M \beta^{(m)} \rho^{(m)}(r)$,
and $r=\sqrt{\left(x-x_m\right)^2+\left(z-z_m\right)^2}$ is the radial distance from each collocation point.
In this study, the multiquadric (MQ) radial basis function was selected, defined as:
$\rho^{(m)}(r)=\sqrt{r^2+c^2}$,
where, c is a shape parameter optimized to balance accuracy and numerical stability. The MQ-RBF was chosen because it provides smooth approximations and stable interpolation, even with irregular node distributions—beneficial for modeling infiltration in trapezoidal domains with curved or sloped boundaries. The coefficients $\beta^{(m)}$ were determined by solving the system of linear equations from collocation at boundary and interior points, implemented in MATLAB R2017. The algorithm follows the classical DRBEM scheme proposed by for non-linear transient problems [49-51].
In this section, we will discuss the results obtained for the problem of infiltration of furrow irrigation water with the Reciprocity Boundary Element Method for the Pima clay loam soil type.
4.1 Governing equation
DRBEM is an effective and stable approach for solving the modified Helmholtz equation, supporting its use in modeling the transformed Richards equation in this study [27]. We obtain
$\begin{aligned} & \frac{\partial \theta}{\partial T}=-\nabla\left(K(\theta)\left[\frac{\partial \psi}{\partial X} i+\left(\frac{\partial \psi}{\partial z}-1\right)\right] j\right) = \frac{\partial}{\partial X}\left(K(\theta) \frac{\partial \psi}{\partial X}\right)+\frac{\partial}{\partial z}\left(K(\theta) \frac{\partial \psi}{\partial z}\right)-\frac{\partial K(\theta)}{\partial z}\end{aligned}$ (6)
Richard’s equation is transformed into an equation in the form linear differential equation using the Kirchhoff transform. Instance:
$K=K_0 e^{\alpha \Psi}, \alpha>0 ; \theta=\int_{-\infty}^{\Psi} K(s) d s$ (7)
K is the hydraulic conductivity, represented in an exponential model of the suction potential $\psi$ [31] with the dimensions [length]. Those of K are [length] [time] ${ }^{-1}$ and those of $[\theta]$ are $[ { length }]^2[{ time }]^{-1}$. $\alpha$ is a measure of gravity relative to the capillarity particular porous material [52].
When $\theta$ forms Matric Flux Potential (MFP), so as
$\frac{\partial \theta}{\partial X}=\frac{\partial}{\partial X} \int_{-\infty}^\psi K(s) d s=K \frac{\partial \psi}{\partial X}$ (8)
later on,
$\frac{\partial \theta}{\partial Z}=\frac{\partial}{\partial Z} \int_{-\infty}^\psi K(s) d s=K \frac{\partial \psi}{\partial X}$ (9)
Step down toward Z, obtain
$\frac{\partial K}{\partial Z}=\frac{\partial(\alpha \theta)}{\partial \theta} \frac{\partial \theta}{\partial Z}=\alpha \frac{\partial \theta}{\partial Z}$ (10)
Substitute Eqs. (8)-(10) into Eq. (7) to get
$\frac{\partial \theta}{\partial T}=\frac{\partial^2 \theta}{\partial X^2}+\frac{\partial^2 \theta}{\partial Z^2}-\alpha \frac{\partial \theta}{\partial Z}$ (11)
Due to constant water infiltration conditions, $\left(\frac{\partial \theta}{\partial T}=0\right)$ then generates Eq. (11) is a modified Helmholtz equation that shows water infiltration in a saturated irrigation channel Eq. (11) is what will be solved. When getting the value of $\phi$, transform it back into
$\frac{\partial^2 \theta}{\partial X^2}+\frac{\partial^2 \theta}{\partial Z^2}=\alpha \frac{\partial \theta}{\partial z}$ (12)
obtained
$\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial z^2}=\phi$ (13)
where, $\phi$ is a transformed potential function introduced for numerical stability and solvability within the DRBEM framework. Eq. (13) is the modified Helmholtz equation which shows water infiltration in a saturated irrigation channel. Eq. (13) will be resolved. When getting $\phi$ value, transform it back to the initial form. Eq. (13) was solved numerically using DRBEM, yielding the potential distribution $\phi$(X,Z) across the infiltration domain. However, to relate the mathematical solution to physical soil parameters, a two-step transformation was applied. First, the suction potential $\psi$ was reconstructed from $\phi$ using the logarithmic relationship derived from the exponential hydraulic conductivity model:
$\psi=\frac{1}{\alpha} \ln \left(\frac{\alpha \phi e^z v_0 L}{\pi k_0}\right)$ (14)
Next, the corresponding volumetric matric flux potential ($\theta$) was computed from $\psi$ using the Van Genuchten soil-water retention model:
$\theta=\left(\frac{1}{1+(\alpha y)^n}\right)^m\left(\theta_s-\theta_r\right)+\theta$ (15)
where, $\theta_s$ and $\theta_x$ represent the saturated and residual matric flux potentials, respectively, and $\alpha, n, m$ are soil-specific empirical parameters obtained from laboratory calibration for Pima clay loam. In the numerical implementation, $\phi$ values were computed at all boundary and internal nodes using the DRBEM algorithm in MATLAB. These $\phi$ values were then transformed into $\psi$ using Eq. (14), and subsequently into $\theta$ using Eq. (15). The resulting $\theta$ values represent the spatial distribution of soil moisture at different infiltration depths and positions along the furrow channel. Physically, higher $\phi$ corresponds to a greater suction potential ($\psi$) and therefore a lower matric flux potential ($\theta$). Conversely, smaller $\phi$ values indicate zones of higher saturation. This transformation chain $\phi \rightarrow \psi \rightarrow \theta$ ensures that the numerical potential obtained from the Modified Helmholtz Equation can be directly interpreted in terms of soil moisture dynamics.
4.2 Boundary conditions
The shape of the irrigation channel used is trapezoidal because this shape is usually used by farmers. Due to the nature of symmetry, where the coordinates are OXYZ, where O is the center of the channel, and OZ represents the appropriate depth for the area to be analyzed is 0≤X≤L+D and Z>0, where the domains are a semi-infinite region 0≤X≤L+D and Z>0 which is assumed to be R [53].
Then the domain limitation is done by taking z = c, where c is a positive real number on the boundary condition because solving the Boundary Condition Problem (BCP) with a semi-infinite domain is difficult to obtain because this known domain is a semi-infinite region (semi-infinite) R, so the line segments can be defined with $C_1, C_2, C_3, C_4, C_5$, where,
(1) C1: Upper Infiltration Surface
This segment represents the water-soil interface at the base of the trapezoidal furrow where infiltration begins.
$0 \leq x \leq \frac{\alpha L}{\pi}, z=0$
The term $\frac{\alpha L}{\pi}$ corresponds to the horizontal projection of the wetted length, derived from the normalized transformation of the real channel width using the DRBEM’s dimensionless mapping parameter $\alpha$.
(2) C2: Right Trapezoidal Wall
This section extends along the sloping side of the channel, where lateral infiltration occurs.
$\frac{\alpha L}{\pi} \leq x \leq \frac{\alpha}{2}(L+D), z=0$
The slope angle relates to the physical trapezoidal geometry through the ratio $\frac{L}{D}=\tan (\beta)$, where, $\beta$ is the side-wall angle. This configuration ensures continuous hydraulic potential along the wall, allowing for proper flux computation.
(3) C3: Right Vertical Boundary
Represents the lateral limit of the infiltration domain, where the hydraulic gradient becomes negligible.
$x=\frac{\alpha}{2}(L+D), 0 \leq z \leq c$
This boundary is treated as hydraulically closed, assuming no lateral flux beyond this distance.
(4) C4: Left Symmetric Boundary
The left wall of the trapezoidal channel is defined as the mirror image of the right wall relative to the vertical axis at x = 0. The lower limit of this boundary, z = $\frac{3 \alpha L}{4 \pi}$, arises from the geometric projection of the left trapezoidal wall when mapping the physical channel into DRBEM’s dimensionless coordinate system.
The derivation proceeds as follows: Let the physical depth of the channel be h, and the side slope ratio be $\tan (\beta)=\frac{L}{h}$.
The mapping between the physical coordinates (x,z) and the transformed DRBEM coordinates ($x^{\prime}, z^{\prime}$) is given by:
$x^{\prime}=\frac{\alpha x}{\pi}, \quad z^{\prime}=\frac{\alpha z}{\pi}$.
At the left wall, the intersection of the inclined wall and the lower truncation plane (z = c) satisfies:
$z^{\prime}=\frac{3}{4} \frac{\alpha L}{\pi}, z=\frac{3 \alpha L}{4 \pi}$
This scaling ensures the angular projection of the left boundary corresponds to 75% of the wetted perimeter’s vertical projection, a geometric approximation validated in previous DRBEM infiltration models [54].
Hence, the segment is defined as:
$x=0, \frac{3 \alpha L}{4 \pi} \leq z \leq c$.
(5) C5: Bottom Truncation Boundary
This boundary closes the domain at finite depth and represents the lower limit of water penetration.
$0 \leq x \leq \frac{\alpha}{2}(L+D), z=c$.
It assumes negligible hydraulic gradient $\partial \phi / \partial z=0$ beyond this depth, implying no downward flux.
Suppose $C=C_1 \cup C_2 \cup C_3 \cup C_4 \cup C_5$ so that it can be obtained BPC for water infiltration in channel irrigation channels in a dimensionless variable with an $R$ domain that is closed and limited by a $C$ curve, namely
$\frac{\partial^2 \phi(x, z)}{\partial x^2}+\frac{\partial^2 \phi(x, z)}{\partial z^2}=\phi(x, z)$
With boundary conditions,
$\begin{gathered}\frac{\partial \phi}{\partial n}=\phi n_2+\frac{2 \pi}{\alpha L} e^{-z}, \text { for } C_1 \\ \frac{\partial \phi}{\partial n}=-\phi, \text { for } C_2 \\ \frac{\partial \phi}{\partial n}=0, \text { for } C_3 \\ \frac{\partial \phi}{\partial n}=0, \text { for } C_4 \\ \frac{\partial \phi}{\partial n}=-\phi, \text { for } C_5\end{gathered}$
Therefore, the water infiltration model can be formed in channel irrigation which has a BCP with a governing equation of the form Modified Helmholtz [54].
4.3 Boundary conditions
DRBEM will be used to solve water infiltration in furrow irrigation, according to the form of the equation that has been given for this problem, namely the modified Helmholtz equation of the form
$\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial z^2}=\phi$
Substitute the boundary conditions into the boundary integral equation that has been given, namely
$\begin{gathered}\lambda(\xi, \eta) \phi(\xi, \eta)=\int_{C_1}\left(\phi(x, z) \frac{\partial \Phi(x, z ; \xi, \eta)}{\partial n}-\Phi(x, z ; \xi, \eta)\left(\frac{2 \pi}{\alpha L} e^{-z}+n_2^{(k)} \phi(x, z)\right)\right) d(x, z) \\ +\int_{C_2 \cup C_5}\left(\phi(x, z) \frac{\partial \Phi(x, z ; \xi, \eta)}{\partial n}-\Phi(x, z ; \xi, \eta)(-\phi(x, z)) d s(x, z)+\int_{C_3 \cup C_4}\left(\phi(x, z) \frac{\partial \Phi(x, z ; \xi, \eta)}{\partial n}\right) d s(x, z)\right) \\ +\iint_R \Phi(x, z ; \xi, \eta) \phi(x, z) d x d z\end{gathered}$ (16)
where, $\Phi(x, z ; \xi, \eta)$ is the fundamental solution of the twodimensional Laplace equation, then it can be formed in the Linear Equation System (LES) by discretizing the domain boundaries in a finite number of line segments with several interior points, then with the value $\phi\left(a^{(k)}, b^{(k)}\right)$ you can find the value $\phi(a, b)$, where, $a, b$ is any point on $R \cup C$ with:
$\begin{aligned} \lambda(a, b) \phi(a, b)= & \sum_{k=1}^N \phi^{(k)}\left[f_2^{(k)}(a, b)-v^{(k)} f_1^{(k)}(a, b)\right] -\sum_{k=1}^N f l^{(k)} e^z f_1^{(k)}(a, b) +\sum_{j=1}^{N+L} \mu\left(a, b ; a^{(j)}, b^{(j)}\right) \phi^{(j)}\end{aligned}$ (17)
Then DRBEM is implemented in MATLAB R2017 to solve the problem of furrow infiltration of irrigation water in trapezoidal for Pima clay loam soil type [55], with detailed component values required as follows.
At this stage, the process using MATLAB for the Pima clay loam soil type. Following the conceptual correction that θ represents the matric flux potential (θ) rather than volumetric water content, the results presented in Tables 1-5 were reinterpreted in terms of hydraulic energy distribution and infiltration dynamics within the Pima clay loam profile. The θMFP values express the integrated hydraulic energy driving unsaturated flow, with higher magnitudes indicating zones of higher potential infiltration flux.
Table 1 presents the computed matric flux potential (θ) values for Pima clay loam using a coarse boundary discretization (N = 40, M = 225). The θ ranges between 17.30 and 17.38 cm²/s, indicating a relatively high and spatially non-uniform hydraulic potential across the domain. The variation along both the x- and z-directions suggests the presence of unresolved boundary gradients, a typical behavior in low-resolution DRBEM meshes. These results demonstrate that coarse discretization leads to overestimated infiltration potentials and limited numerical stability, emphasizing the need for finer boundary resolution to capture the true infiltration energy distribution.
Table 1. Theta value on Pima clay loam with N = 40 and M = 225
|
(x, z) |
θ |
(x, z) |
θ |
|
0.15, 2 |
17.37752245 |
0.35, 2 |
17.37744347 |
|
0.15, 4 |
17.36971442 |
0.35, 4 |
17.36963549 |
|
0.15, 6 |
17.36192342 |
0.35, 6 |
17.36184452 |
|
0.15, 8 |
17.35414914 |
0.35, 8 |
17.35407027 |
|
0.15, 10 |
17.34639119 |
0.35, 10 |
17.34631232 |
|
0.15, 12 |
17.33864902 |
0.35, 12 |
17.33857013 |
|
0.15, 14 |
17.33092171 |
0.35, 14 |
17.33084275 |
|
0.15, 16 |
17.32320764 |
0.35, 16 |
17.32312847 |
|
0.15, 18 |
17.31550375 |
0.35, 18 |
17.31542399 |
|
0.15, 20 |
17.30780385 |
0.35, 20 |
17.30772224 |
|
0.25, 2 |
17.37748295 |
0.45, 2 |
17.37740398 |
|
0.25, 4 |
17.36967495 |
0.45, 4 |
17.36959603 |
|
0.25, 6 |
17.36188397 |
0.45, 6 |
17.36180508 |
|
0.25, 8 |
17.35414914 |
0.45, 8 |
17.35403084 |
|
0.25, 10 |
17.34635175 |
0.45, 10 |
17.3462729 |
|
0.25, 12 |
17.33860957 |
0.45, 12 |
17.33853071 |
|
0.25, 14 |
17.33088222 |
0.45, 14 |
17.33080329 |
|
0.25, 16 |
17.32316805 |
0.45, 16 |
17.32308892 |
|
0.25, 18 |
17.31546386 |
0.45, 18 |
17.31538417 |
|
0.25, 20 |
17.30776301 |
0.45, 20 |
17.30768155 |
Table 2. Theta value on Pima clay loam with N = 80 and M = 225
|
(x, z) |
θ |
(x, z) |
θ |
|
0.15, 2 |
16.89771715 |
0.35, 2 |
16.89768803 |
|
0.15, 4 |
16.89093642 |
0.35, 4 |
16.8909073 |
|
0.15, 6 |
16.88416968 |
0.35, 6 |
16.88414055 |
|
0.15, 8 |
16.87741686 |
0.35, 8 |
16.87738772 |
|
0.15, 10 |
16.87067788 |
0.35, 10 |
16.87064873 |
|
0.15, 12 |
16.86395267 |
0.35, 12 |
16.86392351 |
|
0.15, 14 |
16.85724111 |
0.35, 14 |
16.85721193 |
|
0.15, 16 |
16.85054302 |
0.35, 16 |
16.85051383 |
|
0.15, 18 |
16.84385805 |
0.35, 18 |
16.84382883 |
|
0.15, 20 |
16.83718515 |
0.35, 20 |
16.83715585 |
|
0.25, 2 |
16.89770259 |
0.45, 2 |
16.89767348 |
|
0.25, 4 |
16.89092186 |
0.45, 4 |
16.89089274 |
|
0.25, 6 |
16.88415512 |
0.45, 6 |
16.88412599 |
|
0.25, 8 |
16.87740229 |
0.45, 8 |
16.87737315 |
|
0.25, 10 |
16.87066331 |
0.45, 10 |
16.87063416 |
|
0.25, 12 |
16.86393809 |
0.45, 12 |
16.86390893 |
|
0.25, 14 |
16.85722652 |
0.45, 14 |
16.85719734 |
|
0.25, 16 |
16.85052843 |
0.45, 16 |
16.85049924 |
|
0.25, 18 |
16.84384344 |
0.45, 18 |
16.84381423 |
|
0.25, 20 |
16.8371705 |
0.45, 20 |
16.83714122 |
Table 2 presents the calculated matric flux potential (θ) values for Pima clay loam using an increased boundary discretization of N = 80. The θMFPvalues range between 16.84 and 16.89 cm2/s, showing a noticeable reduction compared with the N = 40 configuration. This decrease reflects improved numerical stability and spatial resolution, resulting in smoother and more uniform hydraulic energy distribution throughout the infiltration domain. The more consistent vertical (z-direction) gradients of indicate enhanced accuracy in capturing infiltration potential variations across the soil depth up to 20 cm. These results confirm that a moderate increase in boundary nodes significantly improves the precision of DRBEM simulations without introducing computational instability.
Table 3 displays the matric flux potential () values for Pima clay loam obtained using a finer boundary discretization (N = 160, M = 225). The θ values range from 16.17 to 16.23 cm²/s, indicating further numerical convergence compared with the previous cases (N = 40 and N = 80). The reduction in θMFPvariability across both the x- and z-directions demonstrates that the model achieves higher spatial precision and smoother hydraulic gradients at this resolution. The nearly uniform distribution of θ suggests that the DRBEM solution begins to approach steady-state behavior, reflecting a stable infiltration potential throughout the simulated soil depth. This result highlights that increasing the boundary nodes up to N = 160 provides an optimal balance between computational cost and numerical accuracy.
Table 3. Theta value on Pima clay loam with N = 160 and M = 225
|
(x, z) |
θ |
(x, z) |
θ |
|
0.15, 2 |
16.22589704 |
0.35, 2 |
16.22589603 |
|
0.15, 4 |
16.22037337 |
0.35, 4 |
16.22037236 |
|
0.15, 6 |
16.2148595 |
0.35, 6 |
16.21485849 |
|
0.15, 8 |
16.2093554 |
0.35, 8 |
16.20935439 |
|
0.15, 10 |
16.20386104 |
0.35, 10 |
16.20386003 |
|
0.15, 12 |
16.19837638 |
0.35, 12 |
16.19837537 |
|
0.15, 14 |
16.1929014 |
0.35, 14 |
16.19290038 |
|
0.15, 16 |
16.18743606 |
0.35, 16 |
16.18743504 |
|
0.15, 18 |
16.18198032 |
0.35, 18 |
16.18197931 |
|
0.15, 20 |
16.17653416 |
0.35, 20 |
16.17653315 |
|
0.25, 2 |
16.22589653 |
0.45, 2 |
16.225895502 |
|
0.25, 4 |
16.22037287 |
0.45, 4 |
16.22037185 |
|
0.25, 6 |
16.214859 |
0.45, 6 |
16.21485798 |
|
0.25, 8 |
16.2093549 |
0.45, 8 |
16.20935388 |
|
0.25, 10 |
16.20386053 |
0.45, 10 |
16.20385952 |
|
0.25, 12 |
16.19837587 |
0.45, 12 |
16.19837486 |
|
0.25, 14 |
16.19290089 |
0.45, 14 |
16.19289988 |
|
0.25, 16 |
16.18743555 |
0.45, 16 |
16.18743454 |
|
0.25, 18 |
16.18197982 |
0.45, 18 |
16.1819788 |
|
0.25, 20 |
16.17653366 |
0.45, 20 |
16.17653264 |
Table 4. Theta value on Pima clay loam with N = 200 and M = 225
|
(x, z) |
θ |
(x, z) |
θ |
|
0.15, 2 |
16.05406431 |
0.35, 2 |
16.05406413 |
|
0.15, 4 |
16.04883805 |
0.35, 4 |
16.04883786 |
|
0.15, 6 |
16.04362067 |
0.35, 6 |
16.04362049 |
|
0.15, 8 |
16.03841215 |
0.35, 8 |
16.03841197 |
|
0.15, 10 |
16.03321247 |
0.35, 10 |
16.03321229 |
|
0.15, 12 |
16.02802159 |
0.35, 12 |
16.0280214 |
|
0.15, 14 |
16.02283948 |
0.35, 14 |
16.0228393 |
|
0.15, 16 |
16.01766613 |
0.35, 16 |
16.01766595 |
|
0.15, 18 |
16.01250151 |
0.35, 18 |
16.01250133 |
|
0.15, 20 |
16.00734558 |
0.35, 20 |
16.0073454 |
|
0.25, 2 |
16.05406422 |
0.45, 2 |
16.05406403 |
|
0.25, 4 |
16.04883796 |
0.45, 4 |
16.04883777 |
|
0.25, 6 |
16.04362058 |
0.45, 6 |
16.0436204 |
|
0.25, 8 |
16.03841206 |
0.45, 8 |
16.03841188 |
|
0.25, 10 |
16.03321238 |
0.45, 10 |
16.03321219 |
|
0.25, 12 |
16.0280215 |
0.45, 12 |
16.02802131 |
|
0.25, 14 |
16.02283939 |
0.45, 14 |
16.02283921 |
|
0.25, 16 |
16.01766604 |
0.45, 16 |
16.01766586 |
|
0.25, 18 |
16.01250142 |
0.45, 18 |
16.01250124 |
|
0.25, 20 |
16.007345`49 |
0.45, 20 |
16.05406403 |
Table 4 reports the matric flux potential (θ) values computed for Pima clay loam with N = 200 boundary nodes and M = 225 interior points. The θ values range between 16.02 and 16.05 cm²/s, showing a consistent downward trend compared with the coarser meshes (N = 80 and N = 160). This indicates that the solution has achieved near-convergence, with minimal variation across the spatial domain. The uniform θMFPdistribution across both lateral (x) and vertical (z) directions suggests that the DRBEM formulation effectively stabilizes under this discretization density. These results confirm that N = 200 represents an optimal resolution, providing accurate and computationally efficient infiltration modeling performance without numerical oscillations.
Table 5. Theta value on Pima clay loam with N = 320 and M = 225
|
(x, z) |
θ |
(x, z) |
θ |
|
0.15, 2 |
16.6388 |
0.35, 2 |
16.6388 |
|
0.15, 4 |
16.6325 |
0.35, 4 |
16.6324 |
|
0.15, 6 |
16.6261 |
0.35, 6 |
16.6261 |
|
0.15, 8 |
16.6198 |
0.35, 8 |
16.6198 |
|
0.15, 10 |
16.6135 |
0.35, 10 |
16.6135 |
|
0.15, 12 |
16.6072 |
0.35, 12 |
16.6072 |
|
0.15, 14 |
16.6010 |
0.35, 14 |
16.6010 |
|
0.15, 16 |
16.5947 |
0.35, 16 |
16.5947 |
|
0.15, 18 |
16.5885 |
0.35, 18 |
16.5884 |
|
0.15, 20 |
16.5822 |
0.35, 20 |
16.5822 |
|
0.25, 2 |
16.6388 |
0.45, 2 |
16.6388 |
|
0.25, 4 |
16.6325 |
0.45, 4 |
16.6324 |
|
0.25, 6 |
16.6261 |
0.45, 6 |
16.6261 |
|
0.25, 8 |
16.6198 |
0.45, 8 |
16.6198 |
|
0.25, 10 |
16.6135 |
0.45, 10 |
16.6135 |
|
0.25, 12 |
16.6072 |
0.45, 12 |
16.6072 |
|
0.25, 14 |
16.6010 |
0.45, 14 |
16.6010 |
|
0.25, 16 |
16.5947 |
0.45, 16 |
16.5947 |
|
0.25, 18 |
16.5884 |
0.45, 18 |
16.5884 |
|
0.25, 20 |
16.5822 |
0.45, 20 |
16.5822 |
Table 5 presents the results for N = 320, where θ values vary between 16.58 and 16.64 cm²/s. Interestingly, this configuration shows a slight increase in θMFP compared with N = 200, indicating a minor numerical fluctuation or over-refinement effect. Such an anomaly suggests that excessively dense boundary discretization may amplify rounding or interpolation errors within the DRBEM framework, particularly when the radial basis function parameters are not adaptively tuned. Despite this deviation, the θMFP distribution remains generally smooth, implying overall model stability. However, the non-monotonic convergence between N = 200 and N = 320 highlights the practical limitation of over-discretization, reinforcing that N = 200 provides the most balanced trade-off between accuracy and computational cost.
Figures 1-4 compare the θ values in Pima clay loam with the corresponding water content obtained from varying numbers of boundary points (N = 40, 80, 160, 200; M = 225). The results, presented graphically along the evaluation points X = 0.15, 0.25, 0.35, and 0.45 cm, show the variation of water content with soil depth (Z = 0–20 cm), where the θ value obtained shows the value of matric flux potential or humidity (unit or length square per unit time (cm)). Based on the fluctuations shown by each graph, from Figure 4, the graph of theta value at X = 0.15 to the graph of theta value at X = 0.45, it can be seen that there is a decrease in the theta value for each theta value at different N. This behavior when observed for each graph has similarities. The values are obtained from Tables 1-5, which can then be presented in both the distance from the point X (cm) and the depth indicated by Z cm 0–20 cm.
The pattern of decreasing the matric flux potential value at each evaluation point X(cm) = 0.15, 0.25, 0.35, 0.45 (cm) is shown in Figures 5-10.
Figure 1. The value of theta at X = 0.15 with a value of N different from M = 225 on Pima clay loam
Figure 2. The value of theta at X = 0.25 with a value of N different from M = 225 on Pima clay loam
Figure 3. The value of theta at X = 0.35 with a value of N different from M = 225 on Pima clay loam
Figure 4. The value of theta at X = 0.45 with a value of N different from M = 225 on Pima clay loam
Figure 5. The value of theta at N = 40 with M = 225 on Pima clay loam
Figure 6. The value of theta at N = 80 with M = 225 on Pima clay loam
Figure 7. The value of theta at N = 160 with M = 225 on Pima clay loam
Figure 8. The value of theta at N = 200 with M = 225 on Pima clay loam
Figure 9. The value of theta at N = 320 with M = 225 on Pima clay loam
Figure 10. The value of theta with N = 40 - 320, with a value of X different from M = 225
That can be noticed based on Figures 5-10 graphs of theta value a N = 40, 80, 160, 200, 320 with M = 225 on Pima clay loam which shows the change in value at points X(cm) = 0.15 (cm), 0.25 (cm), 0.35(cm) 0.45(cm) for each of these points. It is noted that there is a pattern of decreasing the value of matric flux potential or humanity as the depth increases Z (cm), which leads to the point of convergence.
5.1 Influence of boundary point discretization (N) on θ values
The analysis across Tables 1-5 demonstrates that the number of boundary discretization points (N) plays a central role in determining the accuracy and stability of the simulated matric flux potential (θ) for Pima clay loam under furrow irrigation. At coarse resolutions (e.g., N = 40), θ values were clearly overestimated, indicating that the model was unable to capture steep hydraulic gradients near the wetting front. This behavior is consistent with findings [56], which illustrated that insufficient boundary refinement in BEM applications can lead to numerical diffusion and inflated transient flow estimates. Similar limitations have been reported in other boundary-based modeling fields, where coarse meshes distort stress or flux gradients and reduce local solution fidelity [57].
As the boundary resolution increased, θ values progressively decreased and stabilized, suggesting enhanced representation of infiltration dynamics and reduced interpolation errors within the DRBEM framework. This trend agrees with the results [58, 59], who demonstrated that finer boundary discretization improves DRBEM accuracy in modeling nonlinear soil–water interactions. Comparable improvements with boundary refinement have been observed in computational acoustics [33], thin-structure mechanics [60], and piezoelectric systems [34], underscoring that the accuracy of boundary-based numerical methods is highly sensitive to boundary density across different scientific disciplines.
However, at the highest refinement level (N = 320), a slight increase in θ values was observed compared with N = 160 and N = 200, breaking the expected monotonic convergence trend. This deviation is not attributed to a physical inconsistency in the infiltration process but rather to numerical over-refinement. When the boundary mesh becomes excessively dense, the influence matrix in DRBEM may suffer from ill-conditioning, amplifying rounding and interpolation errors associated with the RBF approximation. As a result, minor oscillations appear in the computed θ values, even though the underlying hydraulic behavior remains physically stable. This effect has been similarly documented in advanced BEM applications, where excessive node densities reduce numerical robustness rather than enhance it [43, 61, 62]. Hence, the anomaly observed at N = 320 is interpreted as a numerical artifact rather than a modeling error.
From a practical perspective, the results indicate that a boundary discretization within the range of N = 160 - 200 offers the most effective balance between computational efficiency and numerical precision. At this level, the model achieves sufficient spatial resolution to capture infiltration gradients without introducing unnecessary computational overhead. This finding supports the use of DRBEM as a computationally efficient and physically consistent tool for modeling infiltration processes in fine-textured, semi-arid soils, where hydraulic properties are often heterogeneous and spatially variable [63].
5.2 Numerical stabilization and optimal N range
The numerical response of the model to increasing boundary discretization indicates a phase of numerical stabilization, where further refinement of the boundary structure no longer yields significant variations in θ values. This behavior is characteristic of BEM formulations, including the DRBEM, in which an equilibrium state is reached between element density and solution accuracy [25, 64-68]. Once this equilibrium is attained, additional boundary nodes contribute only marginal improvements. Similar stabilization effects have been reported in isogeometric BEM simulations for fluid–structure interaction [57] and in dual transformation techniques [69], where refinement beyond a certain spatial threshold produced diminishing computational returns.
Interestingly, at N = 320, the model exhibits a nonlinear response, with a slight increase in θ compared to N = 200, breaking the earlier smooth stabilization pattern. This fluctuation does not imply a physical inconsistency but rather a numerical artifact arising from over-refinement. Excessive boundary density tends to produce an ill-conditioned influence matrix, amplifying round-off and interpolation errors inherent to the RBF approximation [70-76]. Similar effects have been documented in large-scale acoustic scattering and peridynamic fracture analyses using BEM, where overly fine meshes added computational cost without proportional gains in accuracy [33].
From a practical standpoint, this finding highlights the trade-off between precision and computational efficiency, particularly in arid-region irrigation modeling, where multiple water management scenarios must be simulated rapidly [21, 77]. Physically, the stabilization of θ reflects the resolution limit of the soil-water system in Pima clay loam, where infiltration behavior is governed by pore-scale heterogeneity and capillary control. Once the mesh resolution surpasses the representative scale of hydraulic conductivity variation, further discretization fails to add meaningful predictive value [63, 78]. Comparable trends were observed in transient heat conduction studies employing radial integration BEM, in which mesh refinement beyond the optimal range produced identical thermal fields but considerably higher computational demands.
Considering both numerical and physical constraints, an optimal operational range of N = 160 - 200 is recommended for DRBEM-based infiltration modeling. This range provides a balanced compromise between accuracy, numerical stability, and computational efficiency. The recommendation is consistent with efficiency principles established in classical BEM literature [79] and particularly advantageous for large-scale irrigation networks in heterogeneous soils, where excessive refinement can render simulations computationally prohibitive [18].
5.3 Physical interpretation for Pima clay loam
Pima clay loam is a fine-textured soil characterized by a high clay fraction, low hydraulic conductivity, and strong capillary suction, all of which fundamentally control its infiltration dynamics. These soil properties create steep gradients in matric flux potential (θ) near the wetting front, where the transition from saturated to unsaturated flow occurs rapidly. Within this region, water movement is predominantly governed by capillary rather than gravitational forces, resulting in slow, non-uniform infiltration that is highly sensitive to pore connectivity and local suction potential.
The spatial pattern of θ obtained from the DRBEM simulation reveals that higher θ values are concentrated near the channel base, where water first enters the soil, and the hydraulic head is largest. As the distance from the infiltration source increases—both laterally and vertically—θ values gradually decline, illustrating the attenuation of hydraulic energy due to the combined effects of pore friction, air entrapment, and reduced hydraulic conductivity. This behavior reflects the actual physical mechanism of infiltration in clay-rich soils, where strong capillary forces resist rapid percolation, leading to shallow wetting fronts and slower moisture redistribution. The numerical outcome, therefore, captures the characteristic infiltration signature of Pima clay loam: an initial rapid intake followed by a marked decline in infiltration rate as the soil approaches hydraulic equilibrium.
Physically, this pattern signifies that most of the infiltrated water remains concentrated near the furrow base, contributing primarily to lateral rather than vertical moisture spread. Such lateral dominance explains the frequent observation of shallow wetting profiles in furrow-irrigated systems under semi-arid conditions [48, 80]. The DRBEM results are consistent with this mechanism, demonstrating that once the capillary barrier forms within the fine-textured layer, the infiltration flux transitions from a predominantly vertical to a quasi-horizontal direction. This shift has important implications for irrigation efficiency: beyond a certain infiltration depth, additional application time or water volume yields diminishing returns, as most water redistributes laterally rather than infiltrating deeper into the root zone.
The simulation also emphasizes the nonlinear coupling between suction potential and hydraulic conductivity. In Pima clay loam, even small reductions in θ correspond to large decreases in K(ψ), reinforcing the soil’s low permeability under unsaturated conditions. This explains the sharp curvature of infiltration profiles observed in the results and highlights the critical importance of representing θ accurately to predict field-scale infiltration rates. Such insights are valuable for irrigation design, as they indicate that optimizing furrow spacing and irrigation duration requires models capable of capturing this nonlinear soil-water interaction [63, 78].
Beyond numerical convergence, the DRBEM outcomes illustrate a physically meaningful trend: finer boundary discretization allows the model to reproduce the actual hydraulic structure of infiltration fronts. The refinement from N = 80 to N = 200 improves the resolution of local suction gradients, thereby yielding more realistic hydraulic potential distributions that align with expected soil behavior. Conversely, excessive refinement (N = 320) introduces small oscillations unrelated to the infiltration physics—an indication of numerical sensitivity rather than a deviation in soil behavior. The underlying infiltration process remains stable and consistent with the physical constraints of capillarity and hydraulic resistance inherent in fine-textured soils.
These physical interpretations clarify that infiltration in Pima clay loam is not a uniform diffusion process but an energy-limited phenomenon governed by pore-scale interactions and capillary retention. The DRBEM framework effectively resolves these dynamics by concentrating computational resolution along boundaries where the most significant hydraulic gradients occur. This boundary-based representation provides an efficient and physically grounded method to simulate unsaturated infiltration without relying on excessive domain discretization typical of finite element or finite difference approaches.
In summary, the DRBEM results bridge numerical modeling and soil physics by linking θ variations directly to infiltration mechanisms within the Pima clay loam profile. The analysis reveals that accurate prediction of infiltration behavior requires capturing both the capillary-controlled hydraulic gradients near the furrow and the rapid decline in conductivity with depth. By representing these nonlinearities explicitly, DRBEM provides not only a mathematically stable but also a physically coherent framework for simulating infiltration in fine-textured, semi-arid soils.
The physical interpretation of the DRBEM results for Pima clay loam offers valuable insights for improving the design and operation of surface irrigation systems in semi-arid agricultural environments. The model outcomes indicate that infiltration in fine-textured soils is primarily limited by capillary retention and low hydraulic conductivity, which confines water redistribution to shallow depths. This suggests that traditional long-duration irrigation practices may result in inefficient water use, as a significant portion of the applied water remains near the surface rather than reaching the deeper root zone.
From a hydraulic management perspective, the pronounced lateral movement of water observed in the simulation implies that furrow spacing and inflow rate must be optimized to ensure uniform wetting between adjacent furrows. Excessively wide spacing can lead to dry inter-row zones, while overly narrow spacing increases water loss through surface evaporation and runoff. Model-based calibration using DRBEM can therefore assist in determining the optimal geometric configuration of furrow systems, balancing infiltration depth, lateral spread, and irrigation uniformity.
The strong dependence of infiltration on capillary suction also underscores the importance of soil preparation and surface management. Practices such as soil leveling, maintaining fine aggregates, and reducing crust formation can improve the hydraulic connectivity of surface pores, enhancing initial infiltration rates without altering the intrinsic hydraulic conductivity of the soil. Conversely, compacted or crusted layers at the surface can amplify capillary resistance, delaying infiltration and promoting runoff even under optimal irrigation flow conditions.
From a modeling standpoint, the DRBEM framework provides a practical advantage for adaptive irrigation scheduling. By efficiently resolving transient infiltration dynamics with minimal computational demand, the method can be integrated into decision-support systems that predict infiltration behavior under varying soil moisture and climatic conditions. Such integration would enable farmers and irrigation managers to determine the precise duration and volume of irrigation required for specific field conditions, minimizing losses due to deep percolation and surface evaporation.
Furthermore, the energy-based interpretation of θ offers a physically consistent foundation for coupling infiltration models with evapotranspiration and root-water uptake processes. This coupling could yield a more holistic representation of the soil–plant–atmosphere continuum, improving predictive capabilities under fluctuating climatic and crop conditions. When applied to semi-arid regions such as NTT, where water availability is highly seasonal, such predictive precision is essential for sustainable irrigation management and long-term soil conservation.
In practical terms, the findings advocate for the use of DRBEM as a design and diagnostic tool for furrow irrigation systems. Its ability to capture complex infiltration behavior in layered, fine-textured soils allow for scenario testing—such as varying furrow geometries, soil textures, or infiltration durations—without requiring extensive field trials. This capability can significantly reduce design costs while improving the resilience and efficiency of irrigation networks in water-limited environments.
5.4 Spatial distribution of θ across depth and distance
The spatial analysis from Figures 1-4 demonstrates that θ decreases with both depth (z) and horizontal distance (x) from the irrigation channel. Near the source (x = 0.15), θ values are highest across all depths, reflecting rapid saturation in the upper layers. This pattern mirrors physical infiltration behavior, where gravitational and capillary forces are strongest near the water source. Similar proximity effects have been documented in BEM simulations of local stress concentrations in fracture mechanics [32] and in near-field acoustic propagation [33], where field variable magnitudes decay sharply with distance from the source.
As x increases to 0.25 and 0.35, θ values decline steadily, especially at greater depths, indicating reduced lateral movement of water. This lateral attenuation effect is consistent with empirical furrow irrigation studies in fine-textured soils, where water infiltration efficiency drops sharply with distance from the furrow [81]. From a modeling perspective, the DRBEM captures this decline by accurately resolving the coupled lateral–vertical infiltration flow, similar to how isogeometric BEM captures coupled structural–fluid responses in engineering contexts [57].
At the farthest distance examined (x = 0.45), θ values are lowest across all depths, with moisture penetration minimal beyond the shallow layers. This is indicative of significant water loss to both evaporation and preferential flow near the furrow, leaving distal areas under-irrigated. Such behavior underscores the importance of optimal furrow spacing to ensure uniform water distribution—an insight supported by DRBEM-based infiltration studies for heterogeneous and layered soils.
The combined depth-distance patterns highlight the DRBEM’s strength in resolving multi-dimensional infiltration dynamics in anisotropic and heterogeneous soils. This mirrors BEM’s versatility in capturing multi-scale field variations in applications as diverse as thermal conduction for irrigation management. These results emphasize the need to integrate spatially explicit infiltration models into planning tools, enabling optimized scheduling and spacing to maximize water use efficiency in water-scarce environments.
The convergence behavior observed in θ values across increasing N aligns with established theoretical principles of the DRBEM, which emphasize the relationship between boundary discretization and the accuracy of numerical solutions to partial differential equations. Theoretically, finer boundary discretization improves the approximation of boundary integrals, reduces interpolation errors, and enhances the stability of solutions [45, 48]. This is consistent with the mathematical foundations of BEM in other contexts, such as isogeometric modeling and peridynamic formulations, where increasing boundary element resolution leads to rapid improvements in accuracy up to a point of diminishing returns. The current results reinforce these theoretical principles, demonstrating that beyond N ≈ 200, improvements in model precision are negligible, which is in line with convergence thresholds documented in multi-physics simulations.
From an empirical standpoint, the observed trends in θ distribution with depth and horizontal distance mirror field measurements in furrow irrigation systems on fine-textured soils. Studies in both homogeneous [30] and heterogeneous soils [63] consistently report that moisture content decreases sharply with both increasing depth and distance from the water source, due to gravitational drainage, capillary effects, and lateral flow resistance. The DRBEM results in this study reproduce these field patterns with high fidelity, suggesting that the model is capable of capturing both vertical infiltration dynamics and the lateral redistribution of water in a realistic manner. This empirical validation is critical because it bridges the gap between numerical model predictions and practical irrigation performance assessments in real-world agricultural settings.
These findings have direct practical implications for irrigation management in arid and semi-arid regions. The optimal N range of 160–320 identified here allows for accurate yet computationally efficient simulations, enabling iterative scenario testing for different furrow spacings, irrigation durations, and soil types. Field-based irrigation optimization studies, such as studies [31, 81], highlight the importance of coupling accurate infiltration modeling with operational decision-making to reduce water losses and improve uniformity of moisture distribution. The DRBEM framework, when parameterized within the optimal discretization range, can serve as a predictive tool to determine irrigation schedules that minimize evaporation and percolation losses, particularly in soils like Pima clay loam with slow infiltration rates and high water-holding capacities.
Theoretically, the broader implication of this study is that the convergence patterns observed are not unique to infiltration problems but rather represent a generalizable numerical characteristic of the DRBEM and related Boundary Element Method (BEM) formulations. These principles can be extended to various environmental and engineering systems governed by diffusion–advection mechanisms, including contaminant transport, heat conduction, and groundwater recharge modeling. This notion aligns with the development of advanced fuzzy-stochastic numerical frameworks, such as those employing dual-Wiener processes to address uncertainty and randomness in differential equation systems [82]. Similarly, comparative numerical analyses of torpedo-shaped and cubic symmetrical autonomous underwater vehicles under variable marine hydrodynamic conditions [83] highlight the adaptability of numerical modeling strategies when confronted with nonlinear and environment-dependent behaviors. These studies collectively reinforce the methodological flexibility and robustness that DRBEM offers when simulating complex, spatially heterogeneous infiltration domains. Empirically, the present work demonstrates the capability of DRBEM to generate spatially explicit predictions of soil moisture distribution that align with measured field behavior, thereby validating its potential as a reliable decision-support tool for precision irrigation management. The synthesis of theoretical rigor and empirical validation established herein provides a sound foundation for the broader application of DRBEM-based infiltration modeling in both research and practical water resource management contexts.
This study employed the DRBEM to simulate transient water infiltration from trapezoidal furrow channels into Pima clay loam soil, emphasizing the effects of boundary discretization density, numerical convergence, and spatial distribution of the matric flux potential (θ). The principal conclusions are summarized as follows.
First, boundary discretization strongly influenced model precision. Increasing the number of boundary nodes (N) from 40 to 200 improved numerical stability and reduced θMFPvalues, indicating smoother hydraulic potential gradients and enhanced representation of infiltration dynamics. However, when N was further increased to 320, a slight rise in θMFP values was observed, deviating from the previously monotonic convergence trend. This anomaly likely arises from over-refinement, where excessive boundary discretization amplifies interpolation and rounding errors inherent to the RBF approximation. Thus, N = 160-200 was identified as the optimal range, balancing computational cost and numerical accuracy.
Second, the refined DRBEM configuration successfully represented the infiltration characteristics of fine-textured soils such as Pima clay loam, which exhibit slow, non-linear water movement governed by strong capillary suction. Spatially, θ decreased consistently with both soil depth and lateral distance from the water source, confirming the method’s ability to reproduce realistic unsaturated flow patterns under semi-arid soil conditions.
The results collectively suggest that DRBEM, when applied with moderate boundary refinement, offers a computationally efficient and physically consistent framework for infiltration modeling in semi-arid agricultural systems. Its boundary-only formulation provides flexibility in handling irregular geometries and variable boundary conditions, making it a promising numerical tool for irrigation design and optimization.
Despite the encouraging outcomes of this study, several limitations should be acknowledged. The simulations assumed layer-wise homogeneity within the Pima clay loam profile, whereas real field conditions typically display vertical and horizontal variability in hydraulic conductivity, which can significantly influence infiltration behavior. Furthermore, the boundary conditions were considered static and time-invariant, excluding the effects of variable inflow rates, evaporation, and transpiration that occur under natural irrigation cycles. Another key limitation is the absence of experimental or field validation, meaning that the simulated matric flux potential values were not directly compared with measured infiltration data. Additionally, the RBF used in the DRBEM formulation was fixed rather than adaptively optimized, which may have contributed to the minor numerical fluctuation observed at N = 320. Acknowledging these constraints provides a constructive foundation for refining the model and guiding future improvements in both numerical implementation and field applicability.
Future studies should extend the DRBEM framework to incorporate dynamic boundary conditions, root-water uptake, and temporal soil moisture variation to enhance physical realism. Further, coupling DRBEM with climate-responsive evapotranspiration models and validating its predictions through field-scale monitoring will be critical to establishing its practical applicability. Exploring multi-layered and heterogeneous soil structures with spatially variable parameters will also help generalize the method for broader agricultural use. Ultimately, integrating DRBEM into decision-support systems for irrigation management could enable more accurate water allocation, reduce conveyance losses, and improve sustainability in water-scarce regions.
The authors gratefully acknowledge the academic guidance and institutional support provided by the Department of Mathematics, Faculty of Science and Engineering, Universitas Nusa Cendana. The constructive feedback from internal reviewers and colleagues contributed significantly to the development of this research.
This research was entirely self-funded, with no financial support received from any institution or external organization. All simulations and analyses were independently conducted using secondary data, hypothetical and theoretical models, ensuring academic integrity and compliance with ethical research standards.
Maria Lobo was responsible for the overall manuscript preparation, including organizing the research framework, compiling results, and ensuring the paper’s coherence and academic structure.
Petrus Moensaku conducted the numerical coding and computational analysis, implementing the Dual Reciprocity Boundary Element Method (DRBEM) for modeling water infiltration in the Pima clay loam soil.
Waltrudis Asa contributed to the drafting of the manuscript, particularly in developing the introduction, literature review, and discussion of numerical results.
Robertus Guntur carried out proofreading and language editing, ensuring clarity, consistency, and adherence to academic writing standards.
[1] Kay, M., Bunning, S., Burke, J., Boerger, V., Bojic, D., Bosc, P.M., Clark, M., Dale, D., England, M., Hoogeveen, J. (2022). The state of the world’s land and water resources for food and agriculture 2021. Systems at Breaking Point. https://doi.org/10.4060/cb9910en
[2] Kummu, M., de Moel, H., Porkka, M., Siebert, S., Varis, O., Ward, P.J. (2012). Lost food, wasted resources: Global food supply chain losses and their impacts on freshwater, cropland, and fertiliser use. Science of the Total Environment, 438: 477-489. https://doi.org/10.1016/j.scitotenv.2012.08.092
[3] Mishra, R.K. (2023). Fresh water availability and its global challenge. British Journal of Multidisciplinary and Advanced Studies, 4(3): 1-78. https://doi.org/10.37745/bjmas.2022.0207
[4] Salahou, M.K., Jiao, X., Lü, H., Guo, W. (2020). An improved approach to estimating the infiltration characteristics in surface irrigation systems. PloS One, 15(6): e0234480. https://doi.org/10.1371/journal.pone.0234480
[5] Das, S., Pardeshi, S.D. (2018). Integration of different influencing factors in GIS to delineate groundwater potential areas using IF and FR techniques: A study of Pravara basin, Maharashtra, India. Applied Water Science, 8(7): 197. https://doi.org/10.1007/s13201-018-0848-x
[6] Gulhane, A.R., Giri, G.K., Khambalkar, S.V. (2018). Scanning Electron Microscopy (SEM) of seed infected with seed borne fungi. International Journal of Current Microbiology and Applied Sciences, 7(7): 4096-4101. https://doi.org/10.20546/ijcmas.2018.707.476
[7] Jamal, A., Cai, X., Qiao, X., Garcia, L., Wang, J., Amori, A., Yang, H. (2023). Real-time irrigation scheduling based on weather forecasts, field observations, and human-machine interactions. Water Resources Research, 59(12): e2023WR035810. https://doi.org/10.1029/2023WR035810
[8] Lal, I. (2024). In-field optimization of soil water parameters for irrigation scheduling. Master thesis, University of Nebraska.
[9] Seyar, M.H., Ahamed, T. (2024). Optimization of soil-based irrigation scheduling through the integration of machine learning, remote sensing, and soil moisture sensor technology. In IoT and AI in Agriculture. Springer, Singapore. https://doi.org/10.1007/978-981-97-1263-2_18
[10] Ahmed, Z., Gui, D., Murtaza, G., Yunfei, L., Ali, S. (2023). An overview of smart irrigation management for improving water productivity under climate change in drylands. Agronomy, 13(8): 2113. https://doi.org/10.3390/agronomy13082113
[11] Ray, S., Majumder, S. (2024). Water management in agriculture: Innovations for efficient irrigation. Modern Agronomy, 169-185.
[12] Touil, S., Richa, A., Fizir, M., Argente Garcia, J.E., Skarmeta Gomez, A.F. (2022). A review on smart irrigation management strategies and their effect on water savings and crop yield. Irrigation and Drainage, 71(5): 1396-1416. https://doi.org/10.1002/ird.2735
[13] Meriç, M.K. (2025). Implementation of a wireless sensor network for irrigation management in drip irrigation systems. Scientific Reports, 15(1): 14157. https://doi.org/10.1038/s41598-025-97303-w
[14] Salem, H.M., Schott, L.R., Ali, A.M. (2025). Real-time soil moisture monitoring using wireless sensor networks for precision agriculture. In Machine Learning for Environmental Monitoring in Wireless Sensor Networks, pp. 145-172. https://doi.org/10.4018/979-8-3693-3940-4.ch007
[15] Carletti, A. (2017). Trial of protocols and techniques for integrated groundwater management in arid and semi-arid regions to combat drought and desertification. https://iris.unica.it/handle/11584/249607.
[16] Mishra, A., Verbist, K. (2021). Addressing Water Security: Climate Impacts and Adaptation Responses in Africa, Asia, Latin America and the Caribbean: Accomplishment Report. UNESCO Publishing.
[17] Chen, Y., Ma, G., Wang, H., Li, T., Wang, Y. (2019). Application of carbon dioxide as working fluid in geothermal development considering a complex fractured system. Energy Conversion and Management, 180: 1055-1067. https://doi.org/10.1016/j.enconman.2018.11.046
[18] Nurhasanah, A., Manaqib, M., Fauziah, I. (2020). Analysis infiltration waters in various forms of irrigation channels by using Dual Reciprocity Boundary Element Method. Jurnal Matematika MANTIK, 6(1): 52-65. https://doi.org/10.15642/mantik.2020.6.1.52-65
[19] Isidoro, D., Grattan, S.R. (2011). Predicting soil salinity in response to different irrigation practices, soil types and rainfall scenarios. Irrigation Science, 29(3): 197-211. https://doi.org/10.1007/s00271-010-0223-7
[20] Nciizah, A.D., Wakindiki, I.I.C. (2015). Soil sealing and crusting effects on infiltration rate: A critical review of shortfalls in prediction models and solutions. Archives of Agronomy and Soil Science, 61(9): 1211-1230. https://doi.org/10.1080/03650340.2014.998203
[21] Inna, S., Manaqib, M., Samudra, V.D., Erhandi, R. (2024). An analysis of infiltration in furrow irrigation channels with root water uptake. Journal of Applied Mathematics, 2024(1): 1869349. https://doi.org/10.1155/2024/1869349
[22] Zhang, M., Shi, W. (2020). Compositional balance should be considered in soil particle-size fractions mapping using hybrid interpolators. Hydrology and Earth System Sciences Discussions, 2021: 1-25. https://doi.org/10.5194/hess-2020-384
[23] Crompton, O., Sytsma, A., Thompson, S. (2019). Emulation of the Saint Venant equations enables rapid and accurate predictions of infiltration and overland flow velocity on spatially heterogeneous surfaces. Water Resources Research, 55(8): 7108-7129. https://doi.org/10.1029/2019WR025146
[24] Zhu, J., Mohanty, B.P. (2002). Upscaling of soil hydraulic properties for steady state evaporation and infiltration. Water Resources Research, 38(9): 11-17. https://doi.org/10.1029/2001WR000704
[25] Al-Bayati, S.A. (2018). Boundary element analysis for convection-diffusion-reaction problems combining dual reciprocity and radial integration methods. PhD theses, Brunel University London.
[26] Aydin, C., Tezer-Sezgin, M. (2019). The DRBEM solution of Cauchy MHD duct flow with a slipping and variably conducting wall using the well-posed iterations. International Journal of Optimization and Control: Theories and Applications, 9(3): 76-85. https://doi.org/10.11121/IJOCTA.01.2019.00677
[27] Balista, T.G., Loeffler, C.F., Lara, L., Mansur, W.J. (2023). Comparisons between direct interpolation and reciprocity techniques of the boundary element method for solving two-dimensional Helmholtz problems. Engineering Computations, 40(9/10): 2841-2861. https://doi.org/10.1108/EC-06-2023-0290
[28] Katsikadelis, J.T. (2016). The Boundary Element Method for Engineers and Scientists: Theory and Applications. Academic Press.
[29] Liu, Y.J., Mukherjee, S., Nishimura, N., Schanz, M., Ye, W., Sutradhar, A., Pan, E., Dumont, N.A., Frangi, A., Saez, A. (2011). Recent advances and emerging applications of the boundary element method. Applied Mechanics Reviews, 64(3): 30802. https://doi.org/10.1115/1.4005491
[30] Inayah, N., Manaqib, M., Majid, W.N. (2021). Furrow irrigation infiltration in various soil types using Dual Reciprocity Boundary Element Method. AIP Conference Proceedings, 2329(1): 40011. https://doi.org/10.1063/5.0042682
[31] Munadi, M., Rokhman, M.S., Oktaviani, D.N., Ahmadi, A., Jannah, H.R. (2024). Matric flux potential in time independent infiltration problems from a single triangular and a trapezoidal irrigation channel. Jurnal Teori Dan Aplikasi Matematika, 8(1): 285. https://doi.org/10.31764/jtam.v8i1.17033
[32] Andrade, H.C., Trevelyan, J., Leonel, E.D. (2023). Direct evaluation of stress intensity factors and T-stress for bimaterial interface cracks using the extended isogeometric boundary element method. Theoretical and Applied Fracture Mechanics, 127: 104091. https://doi.org/10.1016/j.tafmec.2023.104091
[33] Jiayuan, G., Junying, A.N., Li, M.A., Haiting, X.U. (2022). Numerical quadrature for singular and near-singular integrals of boundary element method and its applications in large-scale acoustic problems. Acta Acustica, 41(5): 768-775.
[34] Gu, Y., Sun, L. (2021). Electroelastic analysis of two-dimensional ultrathin layered piezoelectric films by an advanced boundary element method. International Journal for Numerical Methods in Engineering, 122(11): 2653-2671. https://doi.org/10.1002/nme.6635
[35] Sun, S., Liu, Z., Cheng, A.H.D. (2021). Indirect boundary element method for modelling 2-D poroelastic wave diffraction by cavities and cracks in half space. International Journal for Numerical and Analytical Methods in Geomechanics, 45(14): 2048-2077. https://doi.org/10.1002/nag.3255
[36] Tanaka, M., Herein, A.N. (2025). A boundary element method applied to the elastic bending problem of stiffened plates. WIT Transactions on Modelling and Simulation, 19: 1-10. https://doi.org/10.2495/BE970201
[37] Zou, X., D’Antino, T., Sneed, L.H. (2021). Investigation of the bond behavior of the fiber reinforced composite-concrete interface using the finite difference method (FDM). Composite Structures, 278: 114643. https://doi.org/10.1016/j.compstruct.2021.114643
[38] Ricciu, R., Congedo, P.M., Baccoli, R., Mastino, C.C., Baglivo, C., Buonomo, P.L.M., Popolano, G. (2024). An overview on Finite Different Method (FDM) in multilayer thermal diffusion problem. Journal of Physics: Conference Series, 2685(1): 12054. https://doi.org/10.1088/1742-6596/2685/1/012054
[39] Liszka, T., Orkisz, J. (1980). The finite difference method at arbitrary irregular grids and its application in applied mechanics. Computers and Structures, 11(1-2): 83-95. https://doi.org/10.1016/0045-7949(80)90149-2
[40] Zabala, I., Kelly, T.E., Henriques, J.C.C., Peña-Sanchez, Y., Penalba, M., Blanco, J.M. (2026). Boundary element methods in action: A performance assessment for ocean engineering applications. Ocean Engineering, 343: 123168. https://doi.org/10.1016/j.oceaneng.2025.123168
[41] Fahmy, M.A., Almutlg, A. (2025). Boundary element method solution of a fractional bioheat equation for memory-driven heat transfer in biological tissues. Fractal and Fractional, 9(9): 565. https://doi.org/10.3390/fractalfract9090565
[42] dos Santos, G.A.R., Loeffler, C.F., Bulcao, A., de Oliveira Castro Lara, L. (2025). The direct interpolation boundary element method for solving acoustic wave problems in the time domain. Computational and Applied Mathematics, 44(2): 68. https://doi.org/10.1007/s40314-024-03023-8
[43] Akalin-Acar, Z., Gençer, N.G. (2004). An advanced boundary element method (BEM) implementation for the forward problem of electromagnetic source imaging. Physics in Medicine and Biology, 49(21): 5011-5028. https://doi.org/10.1088/0031-9155/49/21/012
[44] Manaqib, M. (2018). Penyelesaian masalah syarat batas persamaan Helmholtz menggunakan dual reciprocity boundary element method. Jurnal “LOG!K@,” 8(2): 115-132. https://journal.uinjkt.ac.id/logika/article/view/9353.
[45] Solekhudin, I. (2018). A numerical method for time-dependent infiltration from periodic trapezoidal channels with different types of root-water uptake. IAENG International Journal of Applied Mathematics, 48(1): 84-89.
[46] Megasari, M. (2021). Dual Reciprocity Boundary Element Method untuk menyelesaikan masalah infiltrasi stasioner pada saluran datar periodik. Journal of Mathematics Computations and Statistics, 4(1): 30. https://doi.org/10.35580/jmathcos.v4i1.20447
[47] Wang, J., Wei, Z., Ren, C., Cai, B., Wang, Z., Cai, C., Lei, H. (2025). Ficus aerial root-inspired photothermal hygroscopic cellulose-based aerogel with capillarity-hydrated hydrogen bonds water retention structure for atmospheric water harvesting and power generation. Chemical Engineering Journal, 512: 162724. https://doi.org/10.1016/j.cej.2025.162724
[48] Clements, D.L., Lobo, M., Widana, N. (2007). A hypersingular boundary integral equation for a class of problems concerning infiltration from periodic channels. Electronic Journal of Boundary Elements, 5(1): 1-16. https://doi.org/10.14713/ejbe.v5i1.779
[49] Gümgüm, S., Wrobel, L.C. (2017). DRBEM formulation for transient Stokes flow with slip boundary condition. Engineering Analysis with Boundary Elements, 75: 65-78. https://doi.org/10.1016/j.enganabound.2016.12.003
[50] Fahmy, M.A. (2024). A time-stepping DRBEM for nonlinear fractional sub-diffusion bio-heat ultrasonic wave propagation problems during electromagnetic radiation. Journal of Umm Al-Qura University for Applied Sciences, 1-7. https://doi.org/10.1007/s43994-024-00178-2
[51] Al-Bayati, S.A., Wrobel, L.C. (2020). Transient convection-diffusion-reaction problems with variable velocity field by means of Drbem with different radial basis functions. In Computational and Analytic Methods in Science and Engineering, pp. 21-43. https://doi.org/10.1007/978-3-030-48186-5_2
[52] Philip, J. (1984). Steady infiltration from spherical cavities. Soil Science Society of America Journal, 48(4): 724-729. https://doi.org/10.2136/sssaj1984.03615995004800040006x
[53] Lobo, M. (2008). Boundary element method for the solution of problem a class of infiltration. Doctor Thesis, University of Adeleide.
[54] Manaqib, M. (2017). Pemodelan matematika infiltrasi air pada saluran irigasi alur. Jurnal Matematika MANTIK, 3(1): 25-31. https://doi.org/10.15642/mantik.2017.3.1.25-31
[55] Majid, W.N. (2020). Perbandingan karakteristik infiltrasi air irigasi alur dalam berbagai jenis tanah menggunakan Dual Reciprocity Boundary Element Method. Bachelor's thesis, Fakultas Sains dan Teknologi Universitas Islam Negeri Syarif Hidayatullah Jakarta.
[56] Clements, D.L., Lobo, M. (2010). A BEM for time dependent infiltration from an irrigation channel. Engineering Analysis with Boundary Elements, 34(12): 1100-1104. https://doi.org/10.1016/j.enganabound.2010.06.013
[57] Chouliaras, S.P., Kaklis, P.D., Kostas, K.V., Ginnis, A.I., Politis, C.G. (2021). An isogeometric boundary element method for 3D lifting flows using T-splines. Computer Methods in Applied Mechanics and Engineering, 373: 113556. https://doi.org/10.1016/j.cma.2020.113556
[58] Solekhudin, I., Ang, K.C. (2013). A DRBEM for time-dependent infiltration from periodic irrigation channels in a homogeneous soil. Electronic Journal of Boundary Elements, 11(1): 1-12. https://doi.org/10.14713/ejbe.v11i1.1810
[59] Munadi. (2021). Dual Reciprocity Boundary Element Methods (DRBEM) pada Masalah Infiltrasi Stasioner dari suatu Saluran Irigasi Tunggal. In Seminar Nasional Pendidikan Matematika Ahmad Dahlan 2018, pp. 643-652.
[60] Huang, R., Xie, G., Zhong, Y., Geng, H., Li, H., Wang, L. (2022). Boundary element analysis of thin structures using a dual transformation method for weakly singular boundary integrals. Computers & Mathematics with Applications, 113: 198-213. https://doi.org/10.1016/j.camwa.2022.03.014
[61] Mivehchi, A., Harris, J.C., Grilli, S.T., Dahl, J.M., O’Reilly, C.M., Kuznetsov, K., Janssen, C.F. (2017). A hybrid solver based on efficient BEM-potential and LBM-NS models: Recent BEM developments and applications to naval hydrodynamics. In Proceedings of the International Offshore and Polar Engineering Conference, pp. 721-728.
[62] Łach, Ł., Svyetlichnyy, D. (2025). Advances in numerical modeling for heat transfer and thermal management: A review of computational approaches and environmental impacts. Energies, 18(5): 1302. https://doi.org/10.3390/en18051302
[63] Paa, A.D., Solekhudin, I. (2025). Implementation of Dual Reciprocity Boundary Element Method to solve the mathematical model of steady infiltration in heterogeneous soils vertically. Journal of the Indonesian Mathematical Society, 31(3): 1964. https://doi.org/10.22342/jims.v31i3.1964
[64] Zhang, M., Ayala, L.F. (2020). The dual-reciprocity boundary element method solution for gas recovery from unconventional reservoirs with discrete fracture networks. SPE Journal, 25(6): 2898-2914. https://doi.org/10.2118/201250-PA
[65] Samanta, A. (2022). Dual reciprocity boundary element analysis of Helmholtz equations with variable coefficients. Master thesis, Jadavpur University.
[66] Pinheiro, V.P., Loeffler, C.F., Mansur, W.J. (2022). Boundary element method solution of stationary advective–diffusive problems: A comparison between the direct interpolation and dual reciprocity technique. Engineering Analysis with Boundary Elements, 142: 39-51. https://doi.org/10.1016/j.enganabound.2022.05.003
[67] Gümgüm, S. (2010). The Dual Reciprocity Boundary Element Method solution of fluid flow problems. Doctoral dissertation, Middle East Technical University, Turkey.
[68] Ghadiali, S.N., Halpern, D., Gaver III, D.P. (2001). A dual-reciprocity boundary element method for evaluating bulk convective transport of surfactant in free-surface flows. Journal of Computational Physics, 171(2): 534-559. https://doi.org/10.1006/jcph.2001.6792
[69] Zhang, S., Zheng, M., Tang, Y., Zang, R., Zhang, X., Huang, X., Chen, Y., Yamauchi, Y., Kaskel, S., Pang, H. (2022). Understanding synthesis–structure–performance correlations of nanoarchitectured activated carbons for electrochemical applications and carbon capture. Advanced Functional Materials, 32(40): 202204714. https://doi.org/10.1002/adfm.202204714
[70] Yang, Y., Liu, H., Yang, S. (2025). Stock price prediction with variational mode decomposition, ecosystem-based optimization, and radial basis function models: Korea composite index insights. Computational Economics, 1-33. https://doi.org/10.1007/s10614-025-11137-2
[71] Wu, Y., Wang, H., Zhang, B., Du, K.L. (2012). Using radial basis function networks for function approximation and classification. International Scholarly Research Notices, 2012(1): 324194. https://doi.org/10.5402/2012/324194
[72] Wang, Z., Liu, Y., Lang, J., Li, Y. (2025). Successive radial basis function approximation for finite-horizon nonlinear H∞ control problems. IEEE Control Systems Letters, 9: 1045-1050. https://doi.org/10.1109/LCSYS.2025.3579408
[73] Mai-Duy, N., Tran-Cong, T. (2003). Approximation of function and its derivatives using radial basis function networks. Applied Mathematical Modelling, 27(3): 197-220. https://doi.org/10.1016/S0307-904X(02)00101-4
[74] Fornberg, B., Driscoll, T.A., Wright, G., Charles, R. (2002). Observations on the behavior of radial basis function approximations near boundaries. Computers & Mathematics with Applications, 43(3-5): 473-490. https://doi.org/10.1016/S0898-1221(01)00299-1
[75] Deng, R., Liu, C.Z., Zhang, Y., Qin, L. (2025). Finance informed neural networks with residual gain attention based radial basis functions. In 2025 IEEE 20th Conference on Industrial Electronics and Applications (ICIEA), Yantai, China, pp. 1-6. https://doi.org/10.1109/ICIEA65512.2025.11149257
[76] Alomairy, R., Cao, Q., Ltaief, H., Keyes, D. (2026). Analyzing data sparsity in global and compact support radial basis functions for 3D unstructured mesh deformation. In Asynchronous Many-Task Systems and Applications, Lecture Notes in Computer Science. https://doi.org/10.1007/978-3-031-97196-9_9
[77] Torres, E.O., López Secanell, I. (2024). ArtEFactum project: Student perceptions in the teacher training Master’s Program on a transdisciplinary project in physical education and technology. Retos, 61: 334-343. https://doi.org/10.47197/retos.v61.108360
[78] Munadi, Solekhudin, I., Sumardi, S., Zulijanto, A. (2017). A Dual Reciprocity Boundary Element Method for water infiltration problems in a single flat and trapezoidal irrigation channels. In Proceedings of 1st Ahmad Dahlan International Conference on Mathematics and Mathematics Education Universitas Ahmad Dahlan, Yogyakarta, pp. 94-100.
[79] Tanaka, M., Sladek, V., Sladek, J. (1994). Regularization techniques applied to boundary element methods. Applied Mechanics Reviews, 47(10): 457-499. https://doi.org/10.1115/1.3111062
[80] Vereecken, H., Weihermüller, L., Assouline, S., Šimůnek, J., et al. (2019). Infiltration from the pedon to global grid scales: An overview and outlook for land surface modeling. Vadose Zone Journal, 18(1): 1-53. https://doi.org/10.2136/vzj2018.10.0191
[81] Irene, Y., Manaqib, M., Ramadhanty, V.W., Affriani, A.R. (2024). An analysis of water infiltration in furrow irrigation channels with plants in various types of soil in the special region of Yogyakarta using Dual Reciprocity Boundary Element Method. Jurnal Teori Dan Aplikasi Matematika, 8(3): 780-799.
[82] Saadoon, N.Q., Abduirazzaq, N.T., Hussain, A.S., Mohammed, N.A., Az-Zo’bi, E.A., Tashtoush, M.A. (2025). Numerical solutions for fuzzy stochastic ordinary differential equations using Heun’s method with a dual-Wiener process framework. Mathematical Modelling of Engineering Problems, 12(3): 763-773. https://doi.org/10.18280/mmep.120303
[83] Wahyuadnyana, K.D., Indriawati, K., Darwito, P.A., Aufa, A.N., Tnunay, H. (2023). Comparative numerical analysis of torpedo-shaped and cubic symmetrical autonomous underwater vehicles in the context of Indonesian marine environments. Mathematical Modelling of Engineering Problems, 10(6): 1917-1926. https://doi.org/10.18280/mmep.100601