Fixed point results with applications to nonlinear fractional differential equations

1.   Introduction

Groundwater is one of the important water resources in the world. Groundwater is reserved in aquifers with many different sources, like block mountain recharge, direct infiltration of meteoric water, diffuse infiltration, and so on. All the water recharge will result in the variation of groundwater level and flow discharge in the aquifer. In general, it is not easy to observe and realize the distribution of groundwater flow in aquifers. The interaction between surface recharge and groundwater flow is also difficult to understand. Therefore, it is helpful to simulate the variations of groundwater level and flow discharge at any time and any location by using a mathematical model based on the hydraulic theory of groundwater under different hypothetical scenarios.

In many early studies, groundwater problems are usually investigated in homogeneous and isotropic aquifers. Anderson ^[1] developed an analytical solution and used it to examine the interaction of groundwater and surface water in streams, developing explicit analytical solutions for complex recharge situations. Lee et al. ^[2] used the principle of mass conservation combined with base flow records and steady base flow analysis to estimate the long-term average annual groundwater recharge in Taiwan. Their proposed model did not require complex hydrogeological data, vegetation cover, and detailed land use information. The resulting contours of mean annual runoff, groundwater recharge, and recharge rate fields closely match the topographic distribution of Taiwan. Ke ^[3] proposed that a combination of soil and water assessment tools and modular three-dimensional groundwater flow model can be employed to solve the multi-aquifer groundwater flow problem of the Jhuoshuei alluvial fan in Taiwan. The model determined monthly recharge rates for several aquifers so that daily water flows could be properly simulated. The results showed that the mixed model performed better than each single model and presented a spatial and temporal distribution of groundwater, providing insights for authorities to manage groundwater resources. Kong et al. ^[4] studied groundwater level variations in unconfined aquifers under the influence of finite-thickness vadose zones. A semi-analytical solution considering the effect of land cover was derived to calculate seepage. Due to the presence of vadose zones, the average height of the groundwater table and the amount of exchange between surface water and groundwater were reduced. Zomlot et al. ^[5] estimated average long-term groundwater recharge in aquifers using a spatially distributed water balance model. How to accurately assess the dominant variables of the groundwater recharge is critical in their study. Finally, multiple linear regression analysis was performed to evaluate the effect of watershed characteristics on recharge. Águila et al. ^[6] estimated groundwater recharge from discrete precipitation events in unconfined aquifers via the observed groundwater level records in shallow wells. Such recharge estimates are prone to uncertainty when the recharge response is not instantaneous and groundwater seepage is present. Numerical analysis of uncertainties explains the effect of non-instantaneous recharge on groundwater level changes and the effect of water level changes in rivers adjacent to unconfined aquifers. Mahdavi ^[7] analyzed the dynamic changes of groundwater mounds in aquifers under time-varying recharge. The calculated water mounds seemed to agree well with the numerical results obtained by a finite element method. Furthermore, the factors affecting each component of the groundwater budget were identified with the help of a sensitivity analysis. Kar et al. modeled hydraulic head using geological and hydrological data such as alluvium depth, hydraulic conductivity, piezometric measurements, and pumping well information and pointed out that large groundwater withdrawals increase the risk of land subsidence ^[8]. The comprehensive results could be applied to the safe exploitation of groundwater resources. Sedghi and Zhan ^[9] presented an analytical solution for groundwater dynamics in aquifers under time-varying recharge. Considering the lateral boundary recharge from a constant water head, the groundwater level change was greatly affected.

Recently, Pastore et al. ^[10] analyzed the dynamics of groundwater level changes in aquifers due to intermittent rainfall supplies. The developed model analyzed immediate response of groundwater under rainfall events and was validated by using time-series rainfall data and groundwater levels obtained from monitoring stations. Xin et al. ^[11] proposed a piecewise linear approximation based on a conceptual model of riparian aquifer systems to describe river level fluctuations. The mathematical model was verified by the observational data of the river water level in the study area. Wu and Hsieh developed a complete analytical solution to simulate the impact of any type of recharge distribution on groundwater flow in sloping unconfined aquifers ^[12]. The proposed analytical solution was verified by previous research results, and the practicability of the analytical solution was validated by the data of the groundwater station in Taichung City, Taiwan, in 2012 and 2013. Zheng et al. ^[13] used the Fourier transform method to derive an analytical solution for a horizontal two-dimensional (2D) groundwater flow problem under Dupuit-Forchheimer assumptions. A numerical model further verified the accuracy of the analytical solution, showing that for larger aquifer thicknesses, nonlinear effects can be neglected. Their analytical solution extended the theoretical understanding of groundwater dynamics under rainfall. Hassan et al. ^[14] stated that under the influence of global warming and climate change, rising temperatures and fluctuating rainfall will exacerbate the shortage of different water resources and the deterioration of water quality. In their study, the impact of future global warming and climate change on the natural recharge of groundwater in desert unconfined aquifers was assessed using the modeling tool WetSpass. Tao et al. ^[15] presented integrated models of rainfall, interception, and infiltration under vegetation cover. Their study introduced a novel simulation approach to investigate the relationship between rainfall and recharge that enhances our knowledge about rainfall production and concentration over complex slope conditions. Hsieh and Wu derived an analytical solution to the nonlinear Boussinesq equation by a perturbation method ^[16]. The groundwater level in an unconfined aquifer adjacent to a river was affected by rainfall recharge and river water level simultaneously. They also considered Horton's infiltration law to estimate the groundwater recharge for the simulation of groundwater flow.

On the other hand, recent research about groundwater flow is focused on aquifer heterogeneity and anisotropy. Heterogeneous sloping aquifers are widely found in alluvial plain areas, in which the storage of groundwater is generally large. Because the distribution of surface recharge varies with time and space, it is feasible to apply mathematical models to predict or estimate the spatiotemporal variation of groundwater. Serrano ^[17] proposed a new analytical model for an unsteady restricted radial flow in aquifers with heterogeneity. In the analysis, the Theis solution was used as the initial condition, and the decomposition method was adopted instead of using the complicated perturbation method. The mean and root mean square error of the heterogeneous hydraulic conductivity were estimated. Scheibe and Yabusaki ^[18] presented synthetic hydraulic conductivity fields for numerical simulations of 3D subsurface flow and transport processes considering geological changes in aquifers. The results showed that groundwater transport is mainly affected by hydraulic conductivity. Meier et al. ^[19] developed a numerical model to study the effect of aquifer heterogeneity on the relationship between specific yield and hydraulic conductivity. The hydraulic tests in heterogeneous media showed the underestimation of the geometric means of hydraulic conductivities obtained in short- and medium-term tests. Winter and Tartakovsky ^[20] improved the estimation of hydraulic head by introducing a model of water flow through heterogeneous composite porous media. In composite media, the hydraulic conductivity is given by a random distribution. Hemker and Bakker ^[21] derived analytical solutions for groundwater flow in heterogeneous aquifers and found the streamlines are straight for isotropic aquifers, whereas the streamlines appear in a helical shape for anisotropic aquifers. Further, if the aquifer is multi-layered and heterogeneous, multiple groups of groundwater eddies will be generated, and adjacent eddies rotate in opposite directions. Sanchez-Vila et al. ^[22] reported that heterogeneity is the most prominent feature of hydrogeology and estimated the representative hydraulic conductivity of an aquifer. Numerical results have shown that traditional hydraulic tests yield hydraulic conductivity closely correlated with the effective hydraulic conductivity. Huysmans and Dassargues ^[23] found that complex depositional processes lead to heterogeneity of hydrogeological conditions. The research demonstrated the multipoint geostatistics to real-world applications to determine the impact of complex geological heterogeneity on groundwater flow and transport. The results showed that heterogeneity leads to spatial uncertainty of hydraulic conductivity distribution in an aquifer, resulting in significant uncertainty in groundwater-related calculations. Chuang et al. ^[24] divided the leaky aquifer system into several horizontal layers such as the heterogeneous aquitard and the underlying aquifer. A one-dimensional (1D) analytical model was then developed to describe hydraulic head fluctuations in such a heterogeneous leaky aquifer system. They found that the length and location of the discontinuous aquitard have a significant effect on the magnitude of hydraulic head fluctuations in the lower aquifer. Zlotnik et al. ^[25] considered the effect of anisotropy and proposed a general 2D solution for calculating terrain-driven groundwater flow, including both small- and large-scale hydraulic conductivity. The large-scale anisotropy is caused by changes in hydraulic conductivity with depth. This solution can be applied to a variety of systems with arbitrary hydraulic head distribution and irregular geometric boundaries. The interaction between small- and large-scale anisotropy controls overall groundwater flow and should be incorporated into the analysis of terrain-driven flow. Liang and Zhang ^[26] derived analytical solutions for the groundwater levels in a 1D heterogeneous unconfined aquifer with time-varying water sources and fluctuating river levels. The heterogeneity of the aquifer significantly increases the spatial variability of the aquifer and changes the groundwater level gradient, whereas the impact of river water level changes on groundwater flow are less significant. Das et al. ^[27] proposed a transient semi-analytical solution of the linearized Boussinesq equation describing the development of groundwater mounds in a 2D finite heterogeneous aquifer under vertical recharge. The aquifer consists of two rectangular basins that share a common impermeable or permeable boundary at the midplane. The regional development of groundwater mounds showed that the effect of heterogeneity becomes significant in the short term but not in the long term. Wang et al. ^[28] pointed out that in alluvial aquifers, hydraulic conductivity is usually inhomogeneous, claiming that finer sediments are deposited preferentially in the downstream direction. They found that aquifer heterogeneity may lead to groundwater flow instability, representing an inconsistent relationship between theoretically calculated head decay and hysteresis. Lately, Hsieh and Lee ^[29] proposed an analytical model considering groundwater flow variability in an unconfined aquifer system with multiple vertical soil zones. The model can simulate any vertically multi-layered soil stratification problem under temporally varying recharge.

The aquifer heterogeneity has a significant impact on hydrological balance ^[30]. Current large-scale hydrological models fail to adequately explain this subsurface heterogeneity. Research findings suggest that management strategies in these areas cannot rely on most current predictions of groundwater recharge, as spatial variability and concentration of recharge lead to estimation errors up to fourfold. Joshi et al. ^[31] studied detailed spatiotemporal variations in groundwater storage in northwestern India. By analyzing historical groundwater level data from the field, significant spatial heterogeneity in groundwater level and storage variations was observed. It was also found that the heterogeneity of the aquifer system is highly correlated with patterns of groundwater level variations. Wang et al. ^[32] proposed an approach to address groundwater flow problems in heterogeneous aquifers by using the finite element method. Using stochastic variables to parameterize the uncertainty of groundwater flow medium characteristics, extensive program training is required to solve the problem and obtain statistically meaningful results. The proposed method was rigorously performed in the study, and numerous numerical examples were provided to demonstrate the accuracy and efficiency of the local-global analysis approach.

Based on the above survey, most current analytical models are built based on 1D governing equation, which cannot be applicable extensively. Therefore, this study aims to develop a 2D groundwater model considering heterogeneity and anisotropy in a sloping aquifer with a semi-infinite domain under any type of surface recharge. Various factors, complex recharge distributions, and distinct medium combinations that affect groundwater flow variability in aquifers are demonstrated. The recharge patterns can be various in time and space. The recharge is represented by a series of unit step functions, capable of expressing any form of recharge in real-world scenarios. Multiple simulations are conducted to demonstrate the effects of heterogeneity and recharge patterns on groundwater flow variations. The present 2D analytical solutions are beneficial to realize the extent and shape of groundwater mounds in heterogeneous and anisotropic unconfined aquifers, providing insights into the impact for each hydrologic and hydrogeologic factor.

2.   Methodology

A conceptual semi-infinite rectangular domain illustrated in Figure 1 is presented for the simulation of groundwater in a heterogeneous semi-infinite unconfined aquifer with the thickness $d$ and anisotropic hydraulic conductivity ${K}_{x}$ and ${K}_{y}$ in the $x$ and $y$ directions, respectively. The interface is at $x = {L}_{s}$, and the initial water depth ${h}_{0}$ is constant in the whole aquifer. Presumably, the water level at the boundary $x\to \mathrm{\infty }$ is constant head ${h}_{0}$ and no-flow condition occurs at $y = {L}_{y}$. The surface recharge $r(x, t)$ varies with time and space. Under the Dupuit-Forchheimer assumptions, the streamlines closely align with the sloping impervious bed, and the discharge rate per unit width ${q}_{x}$ [L²/T] of the aquifer along the x-axis can be estimated using the following equations ^[33]:

where ${h}_{s}$ represents the variable height of the water head measured vertically from a horizontal reference point [L], $h$ represents the vertical measurement of water head height from the sloping impermeable bed [L]; the subscripts 1 and 2 denotes zones 1 and 2, respectively. ${\theta }_{x}$ is slope angle in the $x$ direction; $t$ is time.

According to the research hypothesis ${h}_{s} = h+x\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_{x}$, the seepage fluxes per unit width ${q}_{x}$ and ${q}_{y}$ [L²/T] in the 2D heterogeneous aquifer become

Considering the inflow and outflow through the aquifer with recharge, the continuity equation can be written as

where ${S}_{y}$ is specific yield. The recharge $r\left(x, t\right)$ can be expressed as follows:

where $u$(-) denotes the Heaviside function. $M$ denotes the total number of increments in space; $P$ denotes the total number of increments in time.

The 2D nonlinear Boussinesq equation for groundwater flow in heterogeneous aquifers can be obtained by substituting (3)–(6) into (7) and (8).

Bansal ^[34] refers to the work of Marino ^[35] to successively approximate the mean groundwater saturation depth $\stackrel{-}{h}$ of an aquifer using the iterative formula

where ${h}_{t}$ is the height of the groundwater table at the end of time $t$. Therefore, the linearized Boussinesq equation and associated boundary conditions ^[34] can be expressed as

Initial condition:

Boundary conditions:

Introducing the dimensionless variables

into (12)–(20) yields

Initial condition:

Boundary conditions:

Equations (21) and (22) associated with initial and boundary conditions (23)–(29) can be solved analytically by distinct techniques (e.g., Laplace transform, Fourier transform, etc.). In this study, the advective terms in (21) and (22) can be eliminated by the following substitution:

with

Equations (21)–(29) become

Initial condition:

Boundary conditions:

where

Expressing ${H}_{v}$ as the multiplication of the spatial terms $\phi \left(X\right)\psi \left(Y\right)$ and the time term $\mathit{\Gamma} \left(T\right)$ by the method of separation of variables results in

which satisfies the following eigenvalue problem

The general solution of (42) and (43) in the $x$ direction is

Substituting (50) and (51) into (44)–(46) yields

Setting ${c}_{1}\mathrm{sin}\left(\alpha \right)$ and ${c}_{1}\alpha \mathrm{cos}\left(\alpha \right)$ as $u$ and $v$, respectively, (53) as (54) can be rewritten as

Solving (55) and (56) for ${c}_{3}$ and ${c}_{4}$, we have

Substituting (57) and (58) into (50) and (51), the eigen functions become

where $1\le X\le \mathrm{\infty }$.

The eigenfunction in the $y$ direction can be easily found as

with eigenvalue ${\beta }_{n} = \frac{n\pi }{2}, \;\mathrm{a}\mathrm{n}\mathrm{d}n$ is a natural number.

After expanding $R\left(X, Y, T\right)$ and ${H}_{v}(X, Y)$ by the eigenfunctions $\phi \left(X\right)$, ${\psi }_{n}\left(Y\right)$, and ${{R}_{n}}^{\mathrm{*}}\left(T\right)$, we get

where

in which

Substituting (62)–(64) into (42) and (43) results in

where

can be replaced by $-{\alpha }^{2}\phi$ and $-{{\beta }_{n}}^{2}\psi$, respectively, so (65) and (66) can be rewritten as

Equations (67) and (68) were solved as follows:

Substituting (59)–(61) and (69) into (41) results in

Substituting (70) into (30) and (31), we can obtain the 2D groundwater level $H\left(X, Y, T\right)$.

The convergence of analytical solutions obtained by the method of separation of variables depends on the specific characteristics of the discussed differential equation and the imposed boundary conditions. The convergence of all analytical solutions was examined by reaching the accuracy ${10}^{-3}$ (dimensionless), and the required number of the eigenvalues varied slightly with the given parameters, typically around 50.

3.   Results and discussion

In the previous studies, some researchers employed the Laplace transformation method and different linearization techniques to derive the analytical solutions to the nonlinear Boussinesq equation (see ^{[35,36,37,38,39]}). After performing numerical experiments with varying space discretization $\mathrm{\Delta }x$, $\mathrm{\Delta }y$ and time discretization $\mathrm{\Delta }t$, the results show that the calculation scheme is stable if the following criteria are met:

The present solution was verified with the analytical solution presented by ^[34] and the nonlinear numerical solution proposed by ourselves using the finite difference method. The adopted parameters are ${h}_{0} = 3 \;\mathrm{m}$, ${K}_{x} = {K}_{y} = 1.5\;\mathrm{m}/\mathrm{d}$, ${S}_{y} = 0.18$, $d = 10\;\mathrm{m}$, $t = {t}_{d} = 1\;\mathrm{d}$, and $r = 30\;\mathrm{m}\mathrm{m}/\mathrm{d}$ for ${\theta }_{x} = 0°$ and $5°$ in Figure 2. In Figure 2, there exists discrepancy of the water level changes among the three solutions. The difference between the analytical solutions and the numerical solution was conjectured to be due to the linearization of the nonlinear term in the governing equations. For further discussion, an error analysis was performed as shown in Figure 3. Figure 3 reveals that both the analytical solutions are not consistent with the numerical solution, but the L2 norms of both analytical solutions are below 0.05, indicating an insignificant difference between the analytical solutions and the numerical solution. We attributed the discrepancy mainly to the distinct versions of the governing equations besides the linearization. Both the analytical solutions are for the linearized Boussinesq equation, whereas the numerical solution is for the original nonlinear Boussinesq equation. The convergence and stability of the numerical solution must adhere to the constraints outlined in Eqs (71) and (72). Moreover, these equations suggest that smaller $K$ and larger ${S}_{y}$ may enhance the simulation accuracy, though achieving this in real-world scenarios is challenging. Bansal et al. ^[34] analytically solved the Boussinesq equation and found the L2 norm of the solutions exceeding 0.174 in certain cases owing to the linearization technique used.

When the aquifer configuration is not uniformly distributed, the recharge rate infiltrating into the aquifer is distinct for different properties. Figure 4 shows the distribution of groundwater level in two aquifers with different compositions. The discrepancy arising from the nonlinear effect between the analytical and numerical solutions at the interface is more significant, but the overall water level distribution is similar. According to the analysis of L2 norm, the difference for the heterogeneous case (see Figure 4) is smaller than the homogeneous case (see Figure 2). Furthermore, the discrepancy between the analytical solution and the numerical solution is more obvious near the interface ($X = 1$) of the heterogeneous aquifer, with the overall L2 norm lower than 0.03. Figure 5 also shows that the gap is smaller for the zone of higher permeability, but it becomes larger for the zone of lower permeability.

In general, it is difficult to directly observe the variation of groundwater flow through the interface. Figure 6 illustrates the distribution of groundwater levels in the $X$ direction over time across the interface. The parameters for the simulation are ${h}_{0} = 3\;\mathrm{m}$, ${K}_{1x} = {K}_{1y} = 2.5\;\mathrm{m}/\mathrm{d}$, ${K}_{2x} = {K}_{2y} = 1.5\;\mathrm{m}/\mathrm{d}$, ${S}_{y1} = 0.28$, ${S}_{y2} = 0.18$, $d = 10\;\mathrm{m}$, $r = 30\;\mathrm{m}\mathrm{m}/\mathrm{d}$, $Y = 0.5$, and ${L}_{s} = 150$ m. The contours are relatively dense along the interface, indicating that the hydraulic gradient is large, and the water level changes greatly. Further, due to the influence of the aquifer slope, groundwater flows toward the boundary $X = 0$, and thus water level rises over there. The rate of water level rise increases with the slope in the whole region. In addition, T = 1 indicates the end of recharge, and when T > 1, groundwater level near the downstream area (X = 0) rises to the peak for a while and then gradually declines as expected.

In heterogeneous aquifers, permeability varies with distinct composition of soils, and the flow rates in the aquifers become non-uniform, resulting in distinct variations of water levels. Further, the fluctuation of groundwater level is largely influenced by surface recharge patterns. Figure 7 demonstrates the effects of two recharge patterns, unimodal recharge and multimodal recharge, with a total amount of 300 mm on the groundwater level. In Figure 7, the groundwater level fluctuates with different recharge modes over time, but the water level gradually rises in zone 2 because of the lower permeability. The average groundwater level in zone 1 is about 61% of that in zone 2, as shown in Figure 7a. Groundwater level changes under multi-peak recharge patterns appear to be greater than those under unimodal recharge patterns. The second peak recharge is 30% of the first peak one, and the second peak water level is about 50% of the first peak one, as shown in Figure 7b. The average groundwater level in Zone 1 is about 55% of that in zone 2.

Figure 8a indicates that the average groundwater level in zone 2 is about 91% of that in zone 1. This is due to the different textures of the soil within the aquifer, which causes the change of the flow discharge of groundwater in the sloping heterogeneous aquifer. In Figure 8a, because of the low permeability in zone 1, the groundwater in zone 2 will slow down and accumulate slightly at the interface when it enters zone 1. Changes in surface recharge patterns and soil configuration in heterogeneous aquifers can both affect groundwater level fluctuations. In Figure 8b, the average groundwater level in zone 2 is about 93% of that in zone 1. The fluctuations in groundwater levels over time are clearly caused by changes in surface recharge. The recharge rate also directly affects the peak of the groundwater level. When the recharge rate increases by 30%, the peak water level increases by 50%.

In Figure 8, only the time effect of recharge patterns on groundwater level change was considered. This study intends to consider the effect of recharge patterns varying in both space and time simultaneously. Figure 9 shows the effect on the groundwater level change under the spatiotemporally variable recharge patterns. When the recharge areas are relatively dispersed, the groundwater level changes are more consistent with the time-varying recharge. If the recharge in different areas is relatively concentrated in space, the groundwater level will become higher overall due to the superposition effect. The results indicate that the varying recharge over time and space has a significant impact on groundwater flow.

Except for the bottom slope of the aquifer, groundwater flow is greatly affected by aquifer anisotropy, as shown in Figure 10. Considering different hydraulic conductivities in $X$ and $Y$ directions, respectively, the effect of anisotropy of an aquifer on groundwater flow can be realized. The groundwater flow has no preferential direction in Figure 10a because of the aquifer isotropy. When ${K}_{x} > {K}_{y}$ in Figure 10b, the flow rate in the $X$ direction is larger as expected than that in the $Y$ direction, i.e., the flow in the $X$ direction suppresses that in the $Y$ direction. If the effect of slope is also considered, the flow moves much faster in the $X$ direction.

When ${\theta }_{x} = 5°$ as shown in Figure 11a, groundwater gradually flows towards $X = 0$. The flow discharge is defined as "negative", i.e., groundwater flows from zone 2 (sandy loam) to zone 1 (loam). Because the aquifer permeability in zone 1 is lower than that in zone 2, the discharge will be diminished before passing the interface, and the discharge becomes the lowest near $X = 1.2$. After passing the interface, the discharge gradually increases and the flow accelerates rapidly near $X = 0.65$. The phenomenon is especially obvious at the mid-section ($Y = 0.5$) because of no boundary effect. When ${\theta }_{x} = -5°$ as shown in Figure 11b, the flow discharge is defined as "positive". The groundwater steadily flows in zone 1 until it passes the interface into zone 2. The discharge in zone 2 slows down slightly near $X = 1.0-1.3$, and then enhances quickly because of higher permeability.

In the exploration and planning of various water resources, it is indispensable to realize the relationship between groundwater level and flow discharge. Once the water level is observed, the corresponding discharge can be found. Figure 12 can be used to estimate the flow change via the known water level change. The change of groundwater level decreases with the increasing slope, whereas the flow rate increases with the increasing slope. Such figures can be set up for different compositions of aquifers. The figure also indicates that the loops cannot become close for heterogeneous aquifers. Near the interface, the water table will increase more rapidly due to the groundwater flow from the high-permeability zone (zone 1) to low-permeability zone (zone 2). The flow rate of groundwater flow increases with the water level in zone 1. When the groundwater flow reaches the interface between different soil zones in the heterogeneous aquifer, the flow rate does not change, but the water level suddenly rises. After the groundwater flow passes the interface, the flow rate synchronously decreases with decreasing water level in zone 2, but the water level does not decrease as much as it rises in zone 1.

4.   Conclusions

This study presents an analytical model for the simulation of 2D groundwater flow through a conceptual rectangular semi-infinite domain. The variations of groundwater level and flow discharge are closely related to surface recharge, bed slope, and aquifer configuration. This study investigated the effect of variable lateral sources on groundwater flow in a heterogeneous semi-infinite domain. The 2D analytical mathematical model for groundwater was constructed, using the change-of-variable technique, the separation of variables method, and the integral transform.

The verification of the analytical solution with the numerical solution showed the agreement between them, and the discrepancy was conjectured to arise from the linearization of the nonlinear term in the Boussinesq equation. The number of eigenvalues required for the convergence of the analytical solution to reach an accuracy of 10^-3 is about 40–60, depending on the influencing factors, indicating a fast convergence speed.

The medium heterogeneity including the alignment greatly influences the distribution of groundwater level as well as discharge and flow direction, indicating the importance of consideration of heterogeneity. In heterogeneous aquifers, soil composition alters permeability, leading to non-uniform flow rates and substantive water level changes. As the slope increases, groundwater level decreases, but the flow rate increases. Surface recharge patterns strongly influence groundwater level fluctuations. Groundwater levels greatly respond to varying recharge over time and space. We also found that the loops delineating the relationship between discharge and groundwater level for different bottom slopes cannot become close for heterogeneous aquifers. Additionally, setting up the relationship between flow discharge and water level for different media can facilitate our exploration and planning of groundwater resources. The loops of discharge and groundwater level curves appear unclosed for heterogeneous media rather than closed for homogeneous media.

Use of AI tools declaration

The authors declare they have not used artificial intelligence (AI) tools in the creation of this article.

Acknowledgments

This study was financially supported by the National Science and Technology Council of Taiwan under Grant No. MOST 111-2313-B-005-037.

Conflict of interest

The authors declare no competing interests.