GIS-Integrated Groundwater Flow Modeling for Heterogeneous Media: Application to the Calera Aquifer

1. Introduction

Water resources management is a critical issue, particularly in regions where populations are growing. This directly affects the water demand for industrial, agricultural, commercial, and domestic use [1]. This issue is generally addressed without considering environmental implications or presenting sustainable water management, which directly affects regional development.

Recent studies on water availability [2,3] estimate that only $2-3$% of all the Earth’s water is freshwater, whereas most of it is in the oceans. From the freshwater, nearly 30% is groundwater, 69% is in glaciers and ice caps, and $\approx 1$% is in surface water bodies (lakes, rivers, dams, etc.). In arid or semi-arid regions in particular, the greatest availability is in the subsurface. For this reason, groundwater exploitation has increased in Mexico, often with limited knowledge of its features [1]. There are many Mexican cities in semi-arid and arid regions where groundwater constitutes the only drinking water supply [4]. Overall, in Mexico, groundwater comprises almost 40% of the supply, according to the Mexican National Commission on Water [CONAGUA, 5].

Scientific tools have been continually developed within the hydrogeological community [4] to achieve optimal characterization of this vital resource. Such methods involve evaluating groundwater flow directions using chemical and piezometric data, sometimes to determine recharge zones or to analyze the quality of water extracted from the subsoil.

Despite significant contributions in the field of groundwater, there are still regions where water management agencies estimate water budgets using outdated techniques, with limited data, and often without considering a conceptual model that incorporates all the region’s hydrogeological variables [6]. This leads to the misclassification of an administrative aquifer as over-exploited.

Among the elements necessary for the optimal construction of a conceptual hydrogeological model are delineation of the basement, determination of water chemistry, identification of recharge and discharge zones, definition of saturated and unsaturated zones, and determination of subsurface heterogeneities in terms of hydraulic conductivity. For the determination of hydraulic parameters, such as conductivity, analytical techniques have been developed [7] based on ideal solutions that are generally not found in real aquifer scenarios; for example, a homogeneous, infinite-extent aquifer. Numerical techniques for determining this parameter can be statistical or involve a data inversion, which can lead to inadequate models due to the non-uniqueness of inverse problems.

In this work, we developed an integrated scheme that uses geographic information system (GIS) tools to obtain a heterogeneous conductivity model. This scheme was subsequently applied to groundwater modeling to determine the implications for flow systems in a heterogeneous region. Furthermore, a simulation code for confined aquifers was developed that enables direct transition between GIS and groundwater modeling without requiring a graphical interface. This aims to facilitate hydrogeological data inversion with model constraints, based on prior information about surface geology or geophysical data.

This paper is structured into four sections: a description of the study area; an overview of groundwater modeling and GIS integration; groundwater flow results and discussions for homogeneous and heterogeneous media in the Calera Aquifer; and our conclusions.

2. Study Area

The Calera Aquifer is one of the administrative aquifers in the state of Zacatecas established by CONAGUA, which is the Mexican government-dependent organization that manages water in the country. Calera’s boundaries were mainly delimited considering geopolitical factors [5]. It is located on the Mexican Plateau, in the state of Zacatecas, Mexico (Figure 1). Geographically, the aquifer includes 3 of the region’s most populous cities: Zacatecas City (the state’s capital), Fresnillo, and Calera. Within its limits, there is an extensive presence of industrial, mining, commercial, and agricultural activities that, alongside the dense population, make it one of the most studied aquifers in the State. This aquifer’s region is also where we find most of the agricultural activities, predominating irrigation methods by gravity and trickle. Besides, this region is home to one of the most important breweries and one of the most important silver mines in the world.

Most of the research in this aquifer is based on chemical water analyses, from water quality for human consumption [8], water quality for agricultural uses [9], identifying groundwater flow directions [5,10], to determining the presence of recharge zones [11]. Agricultural studies are also available [12], given the extensive use of the soil in the region. However, groundwater modeling studies are scarce. The single study in the region, with 21 piezometric data points, dates back to 1994 [13]; therefore, the data is outdated and unreliable for comparison with modern-day data. Furthermore, geophysical prospecting studies in the region are scarce in the literature. There is a single work that employs Transient-Electromagnetic (TEM) measurements [14], which could lead to a 3D conceptual model in the future.

The main goal of this work is to contribute to the study of water resources in this region by providing a GIS-integrated tool and conceptual models to better understand the aquifer, with the hope of changing the way water is managed in this region and increasing awareness of the water availability problem.

2.1. Hydrogeological background

The Calera aquifer is located in the Sierra Madre Oriental, which is a mountain range in the Mexican Plateau characterized by submarine volcanic and sedimentary deposits. In terms of tectono-stratigraphy, the aquifer lies near the triple junction of the Guerrero Composite, Central, and Oaxaquia Terranes [15], but is mostly within the Central Terrane, where the nature of the basement remains unknown. The geological map in Figure 2 was obtained from shapefiles downloaded from the Mexican Geological Services (SGM).

Regionally, the aquifer presents a horst-graben system [14], with alluvial and lacustrine deposits in the central part of the valley, where almost all well data are located (Figure 2). The sedimentary basin is surrounded at the south by two important geological formations composed of igneous rocks such as Tuff, Rhyolite, and Andesite; the Chilitos Formation [16] at the southwest part of the aquifer and the Zacatecas formation [17] at the southeast, where recharge zones have been identified by isotope analysis [11].

There are no piezometers installed in the region for time-dependent observations, or even for stationary. Usually, all piezometric data are collected with level loggers or conventional water-level meters; hence, there is no information about the well properties in the region. Despite this, water budgets in the aquifer are often performed assuming homogeneous media, without information on the vertical component of flow directions, using simplified groundwater frameworks [7] and, sometimes, a single pumping test experiment with no observation wells.

Static Water Level (SWL) data is available in the CONAGUA databases, which provide a geographical viewer to download piezometric data for a single state and/or administrative aquifer. We downloaded the data within the limits of Figure 2, that is, not only the inner data in the polygon, but also the well data around the aquifer. The database contains sparse information for the years 1996 to 2025, with many missing data points (82.9% are zero). This is due to poorly organized data acquisition campaigns in the region, which occur at irregular intervals and at different well locations.

To estimate the hydraulic head, the well’s surface elevation is needed. For this, we consider the Digital Elevation Model (DEM) of Figure 3 for each data point, then we subtract the elevation from the SWL to obtain the Static Water Elevation (SWE) in meters, as presented in Figure 2, where SWE is shown in the well’s colors. This indicated that groundwater flows from south to north towards the central part of the basin, a hypothesis that will be explored through numerical modeling.

2.2. Digital Elevation Map and Gravity of Calera Aquifer

There are elevation points up to 2700 m.a.s.l. in the aquifer, and lower elevation points on the central part of the basin (2010 m.a.s.l.), as presented on the DEM of Figure 3, obtained from NASA Earthdata free repository. The topography suggests that runoff water flows towards the center of the aquifer, and the recharge zone can be located on the two principal higher parts of the area. The topography is consistent with satellite-gravity data downloaded from the Satellite Geodesy research group at Scripps Institution of Oceanography, at the University of California, San Diego. The higher altitudes in the region, with the higher-density (igneous rocks), produce the highest gravity anomalies on the aquifer. This is an expected feature from a geophysical perspective. There is a negative gravity anomaly in the central part of the medium, where the density is constant (alluvial deposits) and the topography is relatively flat. This could be related to a mass deficit in the medium associated with the presence of extraction wells in the zone (Figure 2); however, at this time, we lack the data to support this hypothesis. Incorporation of additional geophysical information and well logging data is needed.

3. Methodology

3.1. Groundwater Flow Modeling

We consider the partial differential equation for a confined, heterogeneous, and anisotropic aquifer, given by

$$\nabla ·(K\nabla h)+W=S_S{\displaystyle \frac{\partial h}{\partial t}},$$

(1)

where $K\left(\mathbf{x}\right)$ is the hydraulic conductivity tensor, $h(\mathbf{x},t)$ is the hydraulic head (m), $W(\mathbf{x},t)$ is the source and/or sink term (s⁻¹), $S_S\left(\mathbf{x}\right)$ is the specific storage of the porous material (m⁻¹) and t is time (s). For non-orthogonal anisotropy, only the diagonal components of K are considered, therefore, equation (1) can be written as

$${\displaystyle \frac{\partial }{\partial x}}\left[K_{xx}{\displaystyle \frac{\partial h}{\partial x}}\right]+{\displaystyle \frac{\partial }{\partial y}}\left[K_{yy}{\displaystyle \frac{\partial h}{\partial y}}\right]+{\displaystyle \frac{\partial }{\partial z}}\left[K_{zz}{\displaystyle \frac{\partial h}{\partial z}}\right]+W=S_S{\displaystyle \frac{\partial h}{\partial t}},$$

(2)

where $K_{xx}\left(\mathbf{x}\right)$, $K_{yy}\left(\mathbf{x}\right)$ and $K_{zz}\left(\mathbf{x}\right)$ are the hydraulic conductivities along the x, y and z directions (m/s). Equation 2 is the main equation for groundwater modeling, and can be solved using Finite Difference Methods [FDM, 18] or Finite Element Methods [FEM, 19]; in this paper, we only consider FDM. The discretization of Equation 2 yields the following ODE system

$$\mathbf{M}{\displaystyle \frac{d\mathbf{h}}{dt}}+\mathbf{Ah}=\mathbf{q},$$

(3)

where $\mathbf{h}$ is hydraulic head (vectorized), $\mathbf{q}$ is the vector containing all the source-sink terms, and $\mathbf{M}$ and $\mathbf{A}$ are matrices obtained from the application of the stencil on the grid (see Appendix A), and include the hydraulic conductivity, specific storage, and boundary conditions [20]. We solve the matrix system of Equation 3 using Biconjugate Gradient Stabilized Method (BiCGSTAB) [21] for the non-symmetric problem.

In practice, the transient-state model (Equation 3) requires time stepping, which can be performed using implicit or explicit methods [22]. To realistically model a transient-state problem, data must be recorded over several time periods, sometimes in logarithmic increments from seconds to days, or from days to years, in well-known piezometric networks, sometimes requiring in situ pumping test data acquisition [7]. Nevertheless, in the Calera aquifer, there are regions where a well datum, i.e., Static Water Level (SWL), is used instead of piezometric monitoring, with data collected every 5 years. Additionally, the z-component information (well screen and casing depths) is never provided. With respect to pumping test data, in our area of study, there is sometimes no observation well for the experiment, leading to measurement errors, and most of the time, only a single pumping test is performed across the entire aquifer region for water budget evaluations. Given the limited data availability, we limit our study to 2D steady-state modeling; hence, equation 2 can be simplified as

$${\displaystyle \frac{\partial }{\partial x}}\left[K_{xx}{\displaystyle \frac{\partial h}{\partial x}}\right]+{\displaystyle \frac{\partial }{\partial y}}\left[K_{yy}{\displaystyle \frac{\partial h}{\partial y}}\right]+W=0.$$

(4)

Although our implementation is designed for 3D, anisotropic, and steady-state modeling, the 2D simulation is straightforward (see Appendix A).

3.2. Conductivity Model Creation

For the conductivity model, we used the geology of the area of interest (Figure 2), where alluvial deposits are most prevalent in the zone, and volcanic structures surround the southern part of the region. To obtain a conductivity model, Geographical Information System software is needed, such as ArcGIS [23] or QGIS [24]. An alternative for GUI software is Geopandas [25], which has recently become popular among code developers to avoid graphical interfaces [26]. Geopandas is an open-source Python 3 project that manipulates dataframes similarly to the Pandas package, with spatial operations on geometry types (points and polygons). This feature allows us to create a dataframe column with hydraulic conductivity values depending on the rock type in the zone. For this, we take the mean of the range of values of Bear [27] to obtain the hydraulic conductivities presented in Table 1.

We consider mainly 3 hydrogeological environments in our region of study: 1) a high hydraulic conductivity zone for sedimentary deposits (Alluvial and Lacustrine deposits), 2) an intermediate hydraulic conductivity zone for conglomerate, sandstone, and limestone, and 3) a low conductivity region for igneous rocks(mostly rhyolite and andesite) and shale. All zones are based on the mean values of the ranges presented by Bear [27] (see Figure 4a).

The alluvial region consists predominantly of coarse- to medium-grained sands, hence the high-conductivity values. For lacustrine, we set high conductivity values because this is the region where dry lakes and rivers occur in the Calera aquifer. For the intermediate zone, we chose a value that accounts for poor cementation and high porosity. Finally, igneous and non-fractured rocks yielded the lowest conductivity. However, fractured igneous rocks could have conductivities comparable to those of gravel deposits. This would require in situ cartography of faults or fractures in the region, or geophysical data acquisition and processing, such as aeromagnetic or seismic data, to identify fractured rocks. We only consider in this work non-fractured rocks, but the model creation could be easily updated in our workflow if more information becomes available.

To determine the grid resolution, we examine piezometric data in the zone. To change to metric units, we transform the coordinates to WGS 84 Zone UTM 13N. Then, we compute the minimal distance among all well data points in the aquifer, which is 273.31 meters. Using a resolution of 250 meters yields a grid size of $95\times 306$, which we considered low-resolution compared with the geology map in Figure 2. Besides, the minimal distance would ignore effects between adjacent wells. For this reason, a resolution of 100 meters was chosen, i.e., at least 1 cell will be computed between piezometric data points. This resolution leads to a grid size of $489\times 767$, which is negligible for our computational resources.

We wrote Python code using the Geopandas package to transform the GIS file into a grid matrix ($K_{489\times 767}$), using the hydraulic conductivity values from Table 1 and an external CSV file. The model obtained using this method is presented in Figure 4b, and resembles the geology in the zone. This implementation provides a fast, accessible way to update the hydraulic conductivity value if additional information becomes available, such as well logging, geophysical inversion, in situ cartography, or even groundwater inverse modeling for K estimation [28]. The same procedure, with the same grid size, was performed to compute the head-constant matrix where the wells are located, to feed the FDM code.

4. Results and discussions

Here, we present the results from 3 experiments: first, we test our implementation on a simple 2D homogeneous model for isotropic and anisotropic media, comparing the results with those of MODFLOW. Second, we simulate groundwater flow in the Calera Aquifer using a homogeneous conductivity model and compare the results with those from MODFLOW. Finally, we test our GIS-integrated groundwater model using FDM, accounting for heterogeneous conductivity.

4.1. Validation of our scheme

The main purpose of this experiment is to validate our FDM implementation by comparing it with simulations using MODFLOW 6 [29,30] for isotropic and anisotropic media with different grid sizes (Figure 5). The 2D domains of the following tests are $-1000\le x\le 1000$ m and $-1000\le y\le 1000$ m. For boundary conditions, we set a head-constant value $h\left(\mathbf{x}\right)=50$ m for $\partial \mathrm{\Omega }$, and a discharge rate $Q=-1500$ m³/day in the center of the domain, hence the source term is

$$\begin{array}{c}\hfill W=\{\begin{array}{cc}\hfill -1500{\mathrm{m}}^3/\mathrm{day},& \mathrm{if}\mathbf{x}=[0,0]\hfill \\ \hfill 0,& \mathrm{otherwise}.\hfill \end{array}\end{array}$$

(5)

In MODFLOW, the boundary conditions are implemented with the CHD package and the sink/source implementation with the WEL package [18]. Since we use a steady-state problem, the time units do not play an important role in the simulations; hence, we set the hydraulic conductivities to $K_x=K_y=10$ m/day. For a grid size of $21\times 21$, the results for a homogeneous medium are almost the same using MODFLOW (Figure 5a) and our implementation (Figure 5b), where the contours present hydraulic conductivity circles centered at $x=0$ and $y=0$. The differences between simulations are negligible, as seen in the percentage error in Figure 5c.

Then, we simulate for an anisotropic medium with $K_x=10$ m/day and $K_y=20$ m/day, with the same grid size and boundary conditions. Again, the results are similar using MODFLOW (Figure 5d) and our implementation (Figure 5e). In this example, the head contours clearly exhibit the hydraulic conductivity ellipsoid for an anisotropic medium, as presented by Cherry and Freeze [31], i.e., there is faster convergence to the drawdown point in the center along the y-direction than along the x-direction. The percentage error, Figure 5f, is still negligible (up to $\pm 1\times {10}^{-4}$%).

Finally, we calibrate our code with a larger grid, $201\times 201$, for the same anisotropic media. Because of the higher resolution, the head contours are smoother than in the lower-resolution examples (see Figure 5g and Figure 5h). Although the percentage errors are larger than in the low-resolution examples, they are still negligible (Figure 5i).

We present horizontal and vertical profiles of the head for the anisotropic case, using both lower- and higher-resolution simulations for a more detailed comparison. As shown in Figure 6, the results are indistinguishable. These profiles are not to be compared with pumping-test data, because they are for a steady-state regime and do not account for infinite boundaries [7].

4.2. Calera Aquifer groundwater modeling for homogeneous medium

We performed steady-state groundwater modeling of the administrative aquifer in Calera, Zacatecas, Mexico (Figure 1), using the piezometric information from Figure 2 (Static Water Height) and a homogeneous hydraulic conductivity model with $K=1$ m/day. This consideration allows us to compare the experiment with MODFLOW, where homogeneous media is straightforward to implement. Mathematically and numerically, this setup is a Poisson problem; see Equation 2. We set the known piezometric data as the head-fixed value and applied no-flow boundary conditions on the boundaries in our code and in MODFLOW (using the CHD package).

We present the groundwater results in Figure 7, using a grid of $489\times 767$. As in the previous experiment, the results are very similar for both methodologies, showing flow from the center region of the aquifer towards the northeast and southwest regions.

The flow towards NE (inner contour level of 2010 meters) is due to the presence of the deepest water static levels in the region (around 2000 meters). It is consistent with the topography and the gravity anomaly of Figure 3. The flow towards the southwestern part of the model is due to the interaction with the wells in that region (approximately at 22°45’ N and 103°W) for the contour levels 2130 m and outer 2100 m. The simulation results present the highest head value at 23°N and 102.48°W (approximately) on the aquifer, which will be interpreted in the following subsection.

4.3. Calera Aquifer Modeling Heterogeneous Media

With the final calibration of our code, we proceeded to perform groundwater modeling for a heterogeneous medium, taking the hydraulic conductivity model $K_{489\times 767}$ of Figure 4b, based on Table 1. Changing the table values would not require changes to our implementation; hence, it can be easily modified for future work.

Before performing the numerical simulation, we note a head anomaly in the central part of the region in the homogeneous result. This could be interpreted as a point source in the media, such as a recharge zone. However, this finding is not consistent with the topography or previous studies [5,11] in the region. For this reason, we perform an exploratory data analysis of the piezometric information using a box plot, as shown in Figure 8. The box plot presents a single outlier in the CONAGUA database. Since we can not assess the veracity of this outlier (it could be a measurement error or a spring well), we decided to eliminate this datum from the piezometric database.

Then, we perform the modeling using our GIS-integrated code, as shown in the map in Figure 9a, for 2D heterogeneous media. To clarify, we do not set any cells to be inactive in our model; the simulation naturally shows inactive zones in the gray areas, where an umbral $h\left(\mathbf{x}\right)=$ NaN is imposed for $h\left(\mathbf{x}\right)<0$. This inactive area corresponds to regions with the lowest hydraulic conductivity values (Figure 4), which are associated with non-fractured igneous rocks and shale, especially in the southern part of the model.

We consider this result the most important in our simulation, since geology provides natural limits for an aquifer, which should be necessary to define its administrative delimitation, but this is not always the case.

Hence, geology, through appropriate hydraulic-conductivity models, provides non-interaction zones where wells are separated by low h values, such as the region in the southwest part of the model, which is separated by tuff rocks. In this case, the numerical limit of the aquifer agrees with the administrative one, in combination with the low conductivity of the andesite rocks (southeast). These findings highlight the importance of using heterogeneous rather than homogeneous media in aquifer characterization, as the homogeneous simulation (Figure 7) shows groundwater interactions in zones where, in reality, there are none. On the other hand, in the northern part of the model, there is no numerical barrier in the simulation, which we interpret as a continuity of the geological aquifer outside the administrative aquifer of Calera.

With respect to the lowest contour values of head in the simulation, there is a drawdown zone in the north-center part of the domain, as mentioned in literature on the region and consistent with the mass deficit presented in the gravity anomaly of Figure 3b. This persistent result could support the hypothesis of Bredehoeft [6] that the aquifer is in an equilibrium state; however, z-component information is needed to ascertain this.

5. Conclusions

We developed a FDM groundwater modeling scheme to simulate flow in realistic aquifers. The FDM code was validated by comparing results with those of MODFLOW 6 for simple, homogeneous, and heterogeneous examples with different grid sizes. The scheme is integrated with GIS and provides a quick tool for creating hydraulic conductivity models based on surface outcropping geology, generating grid-sized, heterogeneous models tailored to simulation needs. The main goal of this work was to evaluate groundwater flow patterns in the Calera Aquifer. Using a homogeneous medium, numerical flow simulations can lead to an erroneous interpretation of groundwater movement in the subsoil, which is frequently the case in water-budget evaluations in this area by water-management organizations; namely, homogeneous, 1D, and full-horizontal flow are assumed. When heterogeneities are considered, the flows are more consistent with the geology and piezometry, yielding results similar to those of hydrogeochemical and hydraulic works in the area. In low hydraulic zones, the simulations show numerical flow barriers corresponding to igneous rocks, consistent with the aquifer’s administrative boundary. On the other hand, in high hydraulic-conductivity zones, groundwater interactions are evident in piezometric data outside the aquifer’s administrative limits. We conclude that, geologically speaking, the aquifer does not match the boundaries set by the water organism, which considers only geopolitical factors. We recommend implementing our methodology to define these boundaries to better understand groundwater in the region and to use this vital resource more sustainably. Current limitations arise in the quality and availability of the piezometric data. There are no piezometers in the region, where hydrogeological vertical information is often missing. Even though the code is written for 3D transient, anisotropic, and heterogeneous media, the lack of data constrains us from running more realistic simulations. An important issue is incorporating faults and fractures into the model, which we leave for future work. Also, it is imperative to incorporate additional geophysical data to develop a 3D conceptual model and to perform groundwater modeling that accounts for the vertical component.

Appendix A.1.

In this appendix, we describe our Finite Difference implementation for a general 3D groundwater flow problem in a heterogeneous, anisotropic medium, in steady state. First, we expand the derivatives in Equation 2 as follows

$${\displaystyle \frac{\partial K_{xx}}{\partial x}}{\displaystyle \frac{\partial h}{\partial x}}+K_{xx}{\displaystyle \frac{{\partial }^{2}h}{\partial x^2}}+{\displaystyle \frac{\partial K_{yy}}{\partial y}}{\displaystyle \frac{\partial h}{\partial y}}+K_{yy}{\displaystyle \frac{{\partial }^{2}h}{\partial y^2}}+{\displaystyle \frac{\partial K_{zz}}{\partial z}}{\displaystyle \frac{\partial h}{\partial z}}+K_{zz}{\displaystyle \frac{{\partial }^{2}h}{\partial z^2}}=S_S{\displaystyle \frac{\partial h}{\partial t}}-W,$$

(A1)

renaming ${\displaystyle \frac{\partial K_{xx}}{\partial x}}\to q_x$ and $K_{xx}\to K_x$, and similarly for x, y, and z, we get

$$q_x{\displaystyle \frac{\partial h}{\partial x}}+K_x{\displaystyle \frac{{\partial }^{2}h}{\partial x^2}}+q_y{\displaystyle \frac{\partial h}{\partial y}}+K_y{\displaystyle \frac{{\partial }^{2}h}{\partial y^2}}+q_z{\displaystyle \frac{\partial h}{\partial z}}+K_z{\displaystyle \frac{{\partial }^{2}h}{\partial z^2}}=S_S{\displaystyle \frac{\partial h}{\partial t}}-W.$$

(A2)

Now, we applied the centered discretization following Spitzer [32], using

$$\begin{array}{ccc}\hfill {\left.{\displaystyle \frac{\partial h}{\partial x}}\right|}_{i,j,k}& =& {\displaystyle \frac{a_{i-1}^2{h}_{i+1,j,k}+\left(a_i^2-a_{i-1}^2\right)h_{i,j,k}-a_i^2{h}_{i-1,j,k}}{a_{i-1}{a}_i(a_{i-1}+a_i)}}\hfill \end{array}$$

(A3)

$$\begin{array}{ccc}\hfill {\left.{\displaystyle \frac{{\partial }^{2}h}{\partial x^2}}\right|}_{i,j,k}& =& 2{\displaystyle \frac{a_{i-1}{h}_{i+1,j,k}-\left(a_{i-1}+a_i\right)h_{i,j,k}+a_i{h}_{i-1,j,k}}{a_{i-1}{a}_i(a_{i-1}+a_i)}},\hfill \end{array}$$

(A4)

where $a_i=\Delta x_i$ is the grid spacing for the $x-$direction. The FD expression for y and z can be easily applied using $b_j=\Delta y_j$ or $c_k=\Delta z_k$ for j- or k-indexes, respectively. Equation A3 is also applied to the derivatives $q_x$, $q_y$, and $q_z$ from above. Thus, equation (A1) is rewritten as

$$\begin{array}{c}\hfill q_x{\displaystyle \frac{a_{i-1}^2{h}_{i+1,j,k}+\left(a_i^2-a_{i-1}^2\right)h_{i,j,k}-a_i^2{h}_{i-1,j,k}}{a_{i-1}{a}_i(a_{i-1}+a_i)}}+2K_x{\displaystyle \frac{a_{i-1}{h}_{i+1,j,k}-\left(a_{i-1}+a_i\right)h_{i,j,k}+a_i{h}_{i-1,j,k}}{a_{i-1}{a}_i(a_{i-1}+a_i)}}\\ \hfill q_y{\displaystyle \frac{b_{j-1}^2{h}_{i,j+1,k}+\left(b_j^2-b_{j-1}^2\right)h_{i,j,k}-b_j^2{h}_{i,j-1,k}}{b_{j-1}{b}_j(b_{j-1}+b_j)}}+2K_y{\displaystyle \frac{b_{j-1}{h}_{i,j+1,k}-\left(b_{j-1}+b_j\right)h_{i,j,k}+b_j{h}_{i,j-1,k}}{b_{j-1}{b}_j(b_{j-1}+b_j)}}\\ \hfill q_z{\displaystyle \frac{c_{k-1}^2{h}_{i,j,k+1}+\left(c_k^2-c_{k-1}^2\right)h_{i,j,k}-c_k^2{h}_{i,j,k-1}}{c_{k-1}{c}_k(c_{k-1}+c_k)}}+2K_z{\displaystyle \frac{c_{k-1}{h}_{i,j,k+1}-\left(c_{k-1}+c_k\right)h_{i,j,k}+c_k{h}_{i,j,k-1}}{c_{k-1}{c}_k(c_{k-1}+c_k)}}\\ \hfill =S_S{\displaystyle \frac{\partial h}{\partial t}}-W.\end{array}$$

(A5)

Arranging the seven $h-$terms and applying the FDM to the time derivative, i.e., an implicit formulation, this leads to

$$\begin{array}{c}\hfill h_{i-1,j,k}\left[{\displaystyle \frac{2K_x-q_x{a}_i}{a_{i-1}(a_{i-1}+a_i)}}\right]+h_{i+1,j,k}\left[{\displaystyle \frac{2K_x+q_x{a}_{i-1}}{a_i(a_{i-1}+a_i)}}\right]+\\ \hfill h_{i,j-1,k}\left[{\displaystyle \frac{2K_y-q_y{b}_j}{b_{j-1}(b_{j-1}+b_j)}}\right]+h_{i,j+1,k}\left[{\displaystyle \frac{2K_y+q_y{b}_{j-1}}{b_j(b_{j-1}+b_j)}}\right]+\\ \hfill h_{i,j,k-1}\left[{\displaystyle \frac{2K_z-q_z{c}_k}{c_{k-1}(c_{k-1}+c_k)}}\right]+h_{i,j,k+1}\left[{\displaystyle \frac{2K_z+q_z{c}_{k-1}}{c_k(c_{k-1}+c_k)}}\right]-\\ \hfill h_{i,j,k}\left[{\displaystyle \frac{2K_x-q_x{a}_i}{a_{i-1}(a_{i-1}+a_i)}}+{\displaystyle \frac{2K_x+q_x{a}_{i-1}}{a_i(a_{i-1}+a_i)}}+{\displaystyle \frac{2K_y-q_y{b}_j}{b_{j-1}(b_{j-1}+b_j)}}+{\displaystyle \frac{2K_y+q_y{b}_{j-1}}{b_j(b_{j-1}+b_j)}}+{\displaystyle \frac{2K_z-q_z{c}_k}{c_{k-1}(c_{k-1}+c_k)}}+{\displaystyle \frac{2K_z+q_z{c}_{k-1}}{c_k(c_{k-1}+c_k)}}\right]\\ \hfill =S_S{\displaystyle \frac{h_{i,j,k}-h_{i,j,k}^0}{\Delta t}}-W_{i,j,k}^t,\end{array}$$

(A6)

where $h^0$ is the previous hydraulic head field, $\Delta t$ is the time step and $W^t$ is the source term for each t-time. Renaming the coefficients which multiply the $h-$terms,

$$\begin{array}{c}\hfill C_L{h}_{i-1,j,k}+C_R{h}_{i+1,j,k}+C_B{h}_{i,j-1,k}+C_F{h}_{i,j+1,k}+C_D{h}_{i,j,k-1}+C_U{h}_{i,j,k+1}+C_C{h}_{i,j,k}=W_{i,j,k}^t-{\displaystyle \frac{S_s}{\Delta t}}{h}_{i,j,k}^0\end{array}$$

(A7)

where the letters in the stencil stand for Left, Right, Back, Front, Down, Up, and Center, respectively. thus

$$\begin{array}{c}\hfill C_L=\left[{\displaystyle \frac{2K_x-q_x{a}_i}{a_{i-1}(a_{i-1}+a_i)}}\right],C_R=\left[{\displaystyle \frac{2K_x+q_x{a}_{i-1}}{a_i(a_{i-1}+a_i)}}\right],\\ \hfill C_B=\left[{\displaystyle \frac{2K_y-q_y{b}_j}{b_{j-1}(b_{j-1}+b_j)}}\right],C_F=\left[{\displaystyle \frac{2K_y+q_y{b}_{j-1}}{b_j(b_{j-1}+b_j)}}\right],\\ \hfill C_D=\left[{\displaystyle \frac{2K_z-q_z{c}_k}{c_{k-1}(c_{k-1}+c_k)}}\right],C_U=\left[{\displaystyle \frac{2K_z+q_z{c}_{k-1}}{c_k(c_{k-1}+c_k)}}\right],\end{array}$$

and

$$\begin{array}{c}\hfill C_C=-\left[C_L+C_R+C_B+C_F+C_D+C_U+{\displaystyle \frac{S_s}{\Delta t}}\right].\end{array}$$

This discretization results in a system of equations to be solved in the form,

$$\mathbf{Ah}=\mathbf{q},$$

(A8)

we use Biconjugate Gradient Stabilized Method (BiCGSTAB) [21] for the non symmetrical solution of Equation A8, combinend with large-sparse system storage.