Theoretical foundations of network methods: First, it is shown that the governing equation of the network model is the same as the governing equation of groundwater in the saturated state. Therefore, considering this point, instead of solving the partial differential equation, we can solve the network model. Then, the governing equations of the network method are extracted and how to use these relations to solve various problems will be shown [1].

After presenting a comprehensive method for using network models as a tool and a numerical method, with changes such as 1) the way of connecting the nodes to each other, 2) the way of modeling the boundary conditions and 3) the geometric shape of the channels connecting the pores, the accuracy of the presented method and model is improved [2].

 Governing equation of the network method

In this study, it is initially assumed that the porous medium can be simulated by a network model consisting of pipes that are connected to each other only horizontally and vertically. In this model, all pipes that are considered as channels of the porous medium are assumed to be cylinders that all have the same cross-sectional area. It should be noted that if the hydraulic conductivity of different points of the porous medium is different, the cross-sectional area of ​​these cylinders will also be different from each other. another assumption made in the structure of the network model is that no volume is considered for the pores of the porous medium, which are the intersections of the horizontal and vertical pipes. Figure 1 shows the structure of the assumed Rectangular Pore Network Model (RPNM).

The assumptions required to obtain the governing equation for the grid model shown are as follows:

1) The fluid in question is Newtonian and incompressible.

2) The flow is steady and steady.

3) Surface tension can be neglected.

Figure 2 shows the basic element of the RPNM used to obtain the governing equation [3].

The continuity equation for the element shown is written as follows:

Where Q_in and Q_out are equal to the inlet and outlet flow rates, respectively, u and v are equal to the horizontal and vertical velocities, respectively, and a_x and a_y are equal to the cross-sectional areas of the horizontal and vertical pipes, respectively. Given that the cross-sectional areas are the same in the x and y directions:

Since groundwater movement is usually slow, the movement of water inside the RPNM model pipes is also considered slow, and therefore the Hagen-Poiseuille law can be used:

Where v_ij is the fluid velocity in the conduit connecting nodes j and i, l_ij is the length of this conduit, h is the head of each node, γ is the specific gravity of the fluid, and μ is the viscosity of that fluid [4].

Assuming that the lengths of the pipes (y_Δ, x_Δ) are much smaller than the length (L) and width (H) of the porous medium, by substituting equation (4) into equation (3), the following equation is obtained:

and if the coefficient  is assumed to be equal to the soil permeability (K):

Assuming the environment to be homogeneous and uniform and considering d_x and d_y equal, the above equation will be written as follows:

The above equation, which is the governing equation of RPNM, is the same as Laplace's equation. Given that Laplace's equation is the governing equation for fluid flow in a saturated and steady state porous medium, it can be concluded that replacing the porous medium with the considered model will not introduce errors in the calculations.

The modeling and simulation of groundwater flow are crucial components of hydrological and environmental studies. Over the years, several numerical methods have been developed to solve the groundwater flow equation, including the Network Method (NM), Finite Difference Method (FDM), and Finite Element Method (FEM). This section presents a review of key international and national studies related to the application of the Network Method for solving the groundwater equation [5].

International Research Background

This seminal book introduces the theoretical foundations of groundwater flow and presents numerical methods for solving the governing equations, including the Network Method. Bear describes the discretization of the flow domain into a network of interconnected nodes and links, enabling efficient calculation of hydraulic heads and fluxes.

The authors provide a comprehensive overview of modeling techniques, emphasizing the importance of spatial discretization using network structures. They highlight how network-based methods are useful in handling complex boundary conditions and varying aquifer properties.

MODFLOW, developed by the U.S. Geological Survey (USGS), incorporates a discretized grid system that resembles the principles of the Network Method. The model’s modular structure allows for robust simulation of transient and steady-state groundwater flow using interconnected grid networks.

In their study on groundwater flow modeling in heterogeneous aquifers, the authors implemented adaptive network meshes to enhance computational accuracy. Their approach demonstrated the potential of network-based discretization in capturing localized hydraulic behavior.
Source: Water Resources Research

This USGS report evaluates analog and digital network methods for groundwater simulation. The authors emphasize the versatility of electrical network analogs in modeling complex flow domains, especially in early pre-computer-era studies.

 National (Iranian) Research Background

Their study applied the Network Method to simulate groundwater flow in the Mashhad plain aquifer. The results, validated against field observations, demonstrated high accuracy and computational efficiency. This master’s thesis at the University of Tehran explored groundwater flow modeling in karst aquifers using network grids. Their findings showed that irregular network structures better represented the heterogeneity of karstic domains.

Conducted at Isfahan University of Technology, this study evaluated the impact of over-extraction from aquifers using network simulations. The method accurately predicted long-term drawdowns under various stress scenarios [6].

Their comparative study between NM, FEM, and FDM revealed that the Network Method performs favorably in scenarios involving complex boundary conditions and anisotropic aquifer properties.
Journal: Iranian Journal of Soil and Water Research

 Analysis and Current Trends

1- Advantages of the Network Method:

ü  Simplicity in implementation.

ü  Flexibility in mesh design (structured and unstructured).

ü  Compatibility with widely used models such as MODFLOW.

ü  Useful in conceptual and analog modeling.

2-    Recent Advancements:

ü  Development of adaptive and irregular mesh networks.

ü  Integration with artificial intelligence for model calibration and uncertainty analysis.

ü  Enhanced visualization and post-processing tools for network-based models.

Algebraic equation governing the grid method in steady state

In general, the conventional numerical methods for solving groundwater problems are the two finite difference and finite element methods (Bailout, 2008). To solve this type of problem using the two methods mentioned, the problem in question is first divided into different nodes or elements, and then, using the specific techniques of each of these two methods, a linear algebraic relationship is written between each node and its surrounding nodes [7].

1- Internal nodes: To obtain the governing algebraic equation for RPNM, the continuity equation is applied to an arbitrary node (i) as shown in figure 3-3. Since the lengths of the pipes in the x and y directions are considered equal, we call this model the Square PNM (SPNM) network model.

By applying the continuity equation (Equation 9) and using the Hagen-Poiseuille law (Equation 3-10), the governing algebraic equation for the steady-state lattice method (Equation 11) is obtained:

Where:

If the desired network is assumed to be square and all pipes have the same diameter, equation (11) will be as follows:

The above equation is the governing algebraic equation for all internal nodes of the SPNM and is exactly similar to the differential equation of the central finite difference method. The computational representation of equation 3-13 is shown in figure 3-4.

While the utilization of groundwater for storage and energy extraction presents exciting opportunities, it is also fraught with challenges rooted in hydrogeology, engineering, and environmental science [8]. The literature reveals that technical feasibility, hydrological integrity, monitoring capacity, and policy support are essential for successful implementation. Future research should prioritize site-specific assessments, long-term system monitoring, and the integration of renewable energy technologies with groundwater-based infrastructure.

1- Boundary nodes: In order to solve the governing equation, the boundary conditions must be fully specified. Boundary conditions can be divided into three general categories. These three categories are specified in figure (5). It should be noted that to investigate the Newman and Robin conditions, since the flux value is needed, a node located outside the boundary and can be considered as a pipe that causes the flux to pass through that boundary should be considered.

A) Dirichlet condition: In these types of boundaries, the value of the solution on the boundary is known and therefore there is no need to write an algebraic equation on the boundaries. These types of boundaries will be influential in the equation written for the unknown points around them. The governing equation for node (i) shown in figure 5 will be written using equation 13 as follows:

B) Neumann conditions: In this type of boundary, the derivative of the solution (flux) is known . Therefore, algebraic equations must be written for the boundary points. These types of boundaries are first examined for the impermeable case and  then for the general case of a  certain flux.

 Zero flux: This type of condition in the aquifer arises from the presence of an impermeable layer such as slate or a layer with very low permeability such as a layer of compacted clay. In this type of boundary, there will be no flow to the depth of the impermeable layer and therefore the pipe inside the layer can be eliminated (Figure 6).

By using the continuity equation, the governing equation for the node shown in figure 3-6 will be equal to:

Non-zero flux: To model these types of nodes, we consider a hypothetical node outside the problem domain that transfers the desired flux to the node on the boundary [9].

According to equation 4, the flow rate into node (j) shown in figure 7, through virtual node 2 is equal to:

By writing the continuity relation for the node in question, equation 19 is obtained:

Which, after simplification, is written as follows:

In the above equation, by setting β equal to zero, the impervious condition, which is the same as equation 16, will be obtained. It should be noted that equation 19 is the same as the equation obtained from the previous finite difference method [10].

 C) Robin condition: In this type of condition, the linear combination of the solution and its derivative is specified at the boundary . This type of boundary can be used to determine the boundaries of a leaky aquifer.

Using equation 4, the governing algebraic equation for this type of condition for point k can be obtained as follows:

And by writing the continuity equation:

Equation 24 is exactly the governing equation for the Rabin condition through the finite difference method [11].

Reilly and Franke conducted early evaluations of analog and digital network analog methods for groundwater simulation, highlighting their robustness in representing complex flow domains en.wikipedia.org+13sciencedirect.com+13arxiv.org+13.

In “Benchmark for numerical solutions of flow in heterogeneous groundwater formations”, the authors compare several approaches including finite difference, finite element, discontinuous Galerkin, spectral and random walk but also discuss network-based discretization, providing insight into their accuracy and convergence in heterogeneous settings arxiv.org.

Although primarily a deep learning surrogate model, this study uses a network architecture (Attention U‑Net) to predict hydraulic head in heterogeneous aquifers, offering a mesh‑based analog to traditional numerical networks with significant speed advantages arxiv.org.

A comparative hydrogeologic study assessed semi distributed pipe-network and fully distributed finite-difference models (using MODFLOW–USG with Connected Linear Network, or CLN), showing how embedded linear networks (pipes) can represent karst conduits within unstructured grids mdpi.com+3pubmed.ncbi.nlm.nih.gov+3link.springer.com+3.

The GW‑PINN method applies a mesh-free neural network that directly solves the groundwater PDE akin to solving a continuous "network" of points informed by physics demonstrating the versatility of network-like numerical representations sciencedirect.com.

An MDPI study compared one-dimensional analytical solutions with numerical network-style discretization’s in the FEMME environment, highlighting numerical accuracy using cell-based network techniques mdpi.com [12].

Summary: Table of Key Foreign Studies

Emerging Themes

ü  Hybrid & multiscale modeling: Integration of pipe-like conduits or attention-based network layers within larger grid systems.

ü  Benchmarking rigor: Focus on accuracy, convergence, and performance in heterogeneous aquifers.

ü  Surrogate & mesh-free techniques: Modern analogs to network methods aiming to reduce computational costs. 

The effect of heterogeneity and asymmetry on the governing algebraic equations

The least aquifer can be found whose hydraulic conductivity coefficient is completely the same at different points and in different directions. Therefore, the relations obtained in the previous section cannot be applied exactly to this type of aquifer. By comparing relations 6 and 7, it can be concluded that the coefficient in the network method is equivalent to the permeability coefficient (K) in the aquifer. Therefore, different values ​​of K will cause a change in the pipe diameter. Given that in this case the coefficients C_ij are no longer the same, the continuity relation for a node such as i located in a non-homogeneous aquifer (Figure 9) can be written as follows:

For example, if the hydraulic conductivity of the right side of the node in question is half of the hydraulic conductivity of its left side                                        ( ), the algebraic equation governing node (i) is written as follows:

1. Adams & Parkin (2002):

ü  “Development of a coupled surface‑groundwater‑pipe network model for karstic groundwater”
They innovate by coupling a 3D saturated groundwater model (SHETRAN) with a 1D pipe network to simulate flow in karst conduits, surface-subsurface interactions, and preferential bypass flow in epikarst systems dzkjqb.cug.edu.cn+11ijci.avestia.com+11sciencedirect.com+11pmc.ncbi.nlm.nih.gov+2link.springer.com+2link.springer.com+2.

2. Kresic & Panday et al. (2020):

ü  “Evaluation of semi-distributed pipe-network vs. distributed finite-difference models in karst”
Using MODFLOW‑USG with Connected Linear Networks (CLN), the study embeds 1D pipe networks within unstructured grids, simulating turbulent laminar flows in conduits and their exchange with porous matrix cells link.springer.com [13].

3. He, Zhao et al. (2023):

ü  “Parallel groundwater flow simulation based on Discrete Fracture Network (DFN)”
Although focused on fractured rock, this method treats fractures as a network of discrete elements, solved in parallel for efficient simulation of heterogeneous anisotropic flow .

4. Smith et al. (2020 – ArXiv) & Shadab et al. (2021 – ArXiv):

ü  “GW‑PINN (2022)” and “PINNs for steady unconfined flow (2021)”
These Physics‑Informed Neural Networks (PINNs) model groundwater flow without explicit mesh by sampling in space-time. They act like a continuous network of collocation points, solving flow PDEs directly and outperforming traditional discretized solvers en.wikipedia.org+2sciencedirect.com+2arxiv.org+2.

5. Lykkegaard et al. (2020):

ü  “Deep neural network surrogate for uncertainty quantification”
Introduce a deep learning proxy (DNN) that mimics MODFLOW’s behavior. While not a physical network, it uses an internal node-based architecture to represent spatial flow relationships arxiv.org.

6. Vázquez‑Báez et al. (2018):

ü  “Numerical solution of groundwater flow for a 2D aquifer”
Although built using Finite Difference Method, their approach effectively builds a nodal grid system akin to a network of interconnected points solving Darcy’s law arxiv.org+1scirp.org+1.

7. Chen & Chen (2012):

ü  “Gradually varied functions”
They propose a mesh-free numerical technique based on graph theory, where nodes form networks not limited to rectangular grids, ideal for arbitrary geometries link.springer.com+1pmc.ncbi.nlm.nih.gov+1.

8. Dual-path CNN‑MLP (Omar et al., 2023):

“Hybrid DNN for regional-scale flow in Qatar”
Combines convolutional and fully connected layers to predict hydraulic heads across a coarse network of spatial locations, achieving ~99% speed gain over MODFLOW arxiv.org+11iwaponline.com+11ascelibrary.com+11.

Themes & Insights

ü  Pipe/Conduit Network Embedding: Integration of 1D flow paths (e.g., karst pipes, fractures) within 2D/3D network grids enhances model realism.

ü  Mesh-free PINNs & Graph-based Approaches: Nodes are leveraged directly to fulfill the flow equations, bypassing traditional meshing.

ü  Neural-network Surrogates: Use internal node connections to approximate spatial relationships, effectively acting as computational networks.

ü  Parallel & High-performance Models: Discrete fracture network (DFN) and pipe-network models run efficiently with parallel computing.

Injection and Withdrawal

Injection and withdrawal in aquifers are usually carried out using well rings. Due to the small radius of the wells compared to the size of the aquifer, they can be considered as point loads. The continuity equation can be written for the point where injection or withdrawal occurs (point i) as follows:

Where the positive sign indicates injection and the negative sign indicates withdrawal.

The governing algebraic equation of the network method in the unsteady state

By applying equation 13 for the internal nodes and equations 20 and 24 for the boundary nodes, the Poisson equation is transformed into a system of linear algebraic equations using the SPNM method. This means that the operator can be replaced by . To solve the temporal problems, spatial decomposition is performed using the obtained stiffness matrix, while the finite difference method is used for temporal decomposition. Temporal decomposition using the finite difference method can be performed in both forward and backward forms, as well as in the Crank-Nicolson form, which are shown below. The temporal equation governing a confined aquifer in the saturated state is as follows:

 Where K is the permeability coefficient of the porous medium.

 Leading method

Backward method

 Crank-Nicolson method

It should be noted that, like the finite difference and finite element numerical methods, one must be careful in choosing the time step in the forward method so that the solution does not become unstable. In other words, it should be noted that the numerical value of the time step in the forward method is always less than 0.5 in the one-dimensional case and less than 0.125 in the two-dimensional case. It is worth noting that the use of the two backward and Crank-Nicolson methods makes the solution obtained always stable and therefore the numerical value of the time step can be any number [14].

1.2. Nationwide Recharge Decline (2000–2018): Research using MODFLOW and MT3DMS across Iranian basins indicates an alarming drop in groundwater recharge, averaging 39.6 mm/year barely exceeding annual runoff (~32 mm/year) repository.sharif.edu+4civilica.com+4scientiairanica.sharif.edu+4pmc.ncbi.nlm.nih.gov+1link.springer.com+1. This strain is primarily anthropogenic: rising extraction rates, reduced surface-water inputs, land subsidence from sediment compaction, and altered land use. Approximately 56% of rivers have lost streamflow, further reducing recharge pmc.ncbi.nlm.nih.gov.

2.2. Subsidence and Sinkholes: Extensive extraction, especially for agriculture, has triggered land subsidence up to 35 cm/year in areas like Kerman province theguardian.com. Studies using satellite Interferometric Synthetic Aperture Radar (InSAR) confirm that ~3.5% of Iran's land is subsiding significantly, threatening structural integrity and limiting further aquifer utilization reddit.com+3theguardian.com+3pmc.ncbi.nlm.nih.gov+3.

 3. Aquifer Thermal Energy Storage (ATES):

Several Iranian studies evaluate ATES performance and challenges:

1.3. Economic and Environmental Modeling (Tehran, 2017): A key study analyzed ATES in a confined aquifer to meet thermal needs of a residential complex in Tehran. Four configurations were tested: cooling-only ATES; solar collector heating-only ATES; combined heating via solar collectors and ATES; and integrated ATES with heat pump for both needs. Results show rapid economic payback (≈2.41 years for cooling-only, lifecycle costs ~$16,000), and the combined ATES–heat-pump system yielded lowest emissions (~359 ton CO₂/year) civilica.com+6scientiairanica.sharif.edu+6reddit.com+6.

2.3. Parametric Technical Simulation (2013): Ghaebi et al. used finite difference temperature/flow simulations to assess ATES performance under varying aquifer parameters. They concluded the ATES–heat-pump–solar collector combo offers optimal energy recovery and Coefficient of Performance (COP). Multiple aquifer geometries were tested to refine system design.

3.3. Operational Efficiency & Thermal Storage (2017): A follow-up study by Ghaebi and Bahadori examined system operation under seasonal load, assessing thermal energy recovery (~300–120 kW over 60 days), and energy/exergy efficiencies (~35.6% & 27.4%) academia.edu+4researchgate.net+4civilica.com+4.

4.3. Heat Pump–ATES Coupling (2017 Persian): A separate article confirmed ATES’s effectiveness, highlighting a COP of 10 in cooling-only use, and COPs of 17 (cooling) and 5 (heating) when heat pump coupled, achieving annual energy yields of ~8.7 TJ cooling and ~1.9 TJ heating academia.edu+5civilica.com+5repository.sharif.edu+5.

 4. Groundwater Storage & Recovery (ASR):

Iranian researchers are exploring ASR in agricultural and urban water systems:

1.4. Qazvin Aquifer ASR Case (2023): Hosseini Jolfan et al. investigated ASR within Qazvin’s irrigation system using a coupled MODFLOW MATLAB approach.

 5. Geothermal Energy:

While not directly ATES, Iran’s geothermal initiatives inform subsurface energy use:

1.5. Meshkinshahr Geothermal Plant (2010 ongoing): Iran's first geothermal power station a 5 MW plant near Meshkinshahr extracts heat via >3,000 m wells (>250 °^C temperatures), reinjects cooled water, and combines electricity generation with district heating en.wikipedia.org. The model demonstrates direct-energy extraction feasibility but also highlights a need for extensive drilling and capital. Local authorities estimate ~70 billion IRR invested since 2006 [15].

2.5. National Geothermal Potential & Constraints: Iran’s physical geography supports geothermal, with ≈9% of land having potential. The country has prioritized 10 prospective zones, including near Mount Damavand. Despite promise, analyses show high upfront costs, deep drilling requirements, and environmental considerations (e.g., reinjection integrity) arxiv.org+7en.wikipedia.org+7scientiairanica.sharif.edu+7.

 6. Technical, Environmental & Hydrogeological Concerns:

1.6. Aquifer Integrity & Thermal Mixing: ATES requires low groundwater flow rates to avoid heat dispersion. Iran’s heterogeneous aquifers and seasonal recharge fluctuations complicate design and reduce storage efficiency civilica.com+3scientiairanica.sharif.edu+3academia.edu+3.

2.6. Subsidence & Storage Capacity: Rapid land subsidence (e.g., Darabad & Varamin >25 cm/year around Tehran) demonstrates aquifer structural compromise, reducing storage capacity and increasing seismic risk.

3.6. Economic Feasibility: While ATES payback can be rapid in specific contexts (~2.4 years), geothermal requires deep drilling and high investment (~70 billion IRR), making large-scale deployment a financial risk repository.sharif.edu+2scientiairanica.sharif.edu+2researchgate.net+2.

4.6. Regulatory & Institutional Barriers: Iranian studies note that ATES and geothermal are nascent fields needing clearer energy sector regulations, incentives, and institutional infrastructure to scale.

5.6. Monitoring Limitations:

Simulations using MODFLOW or coupled models highlight uncertainty in hydraulic input data, aquifer geometry, and recharge patterns hindering accurate long-term predictions link.springer.com [16].

 7. Synthesis & Research Gaps:

1.7. Integrating Energy Storage & Grid Analytics: While ATES studies have optimized energy flows, they often omit grid-level integration and seasonal climate variability.

2.7. Scale-up & Regional Planning: Beyond demonstration sites (Tehran, Meshkinshahr, Qazvin), there is limited research on region-wide deployment, multi-site interference, or cross-aquifer impacts.

3.7. Multi-Purpose Subsurface Use: Iranian literature seldom addresses combined use (heat storage + water supply + geothermal power) within a single aquifer system.

4.7. Institutional Framework Development: Despite technical viability, wide adoption is constrained by lack of policy, regulatory clarity, and collaboration among ministries, universities, and energy operators.

 8. Recommendations for Advancement:

ü  Site Selection: Prioritize low-flow, thermally stratified confined aquifers for ATES, with robust geological surveys.

ü  Hybrid Systems: Combine solar–heat pump–ATES systems in urban zones; replicate Qazvin-style ASR in agricultural regions to bolster recharge.

ü  Economic Analysis: Conduct full levelized cost of energy (LCOE) studies for geothermal vs. fossil/gas alternatives.

ü  Monitoring Innovations: Integrate sensors (temperature, pressure) linked to MODFLOW/FEFLOW for real-time adjustment.

ü  Policy & Incentives: Introduce tax credits, feed-in tariffs, and water-energy regulations to support ATES and geothermal ventures.

Iranian research demonstrates that groundwater-based energy storage through ATES, ASR, and geothermal is technically feasible and potentially cost-effective in select sites. However, challenges like aquifer subsidence, thermal dispersion, economic risk, and policy vacuums hinder national-level implementation. Addressing these requires multidisciplinary strategies, combining hydrogeology, energy systems engineering, environmental monitoring, and governance reforms.

There is substantial potential in expanding pilot projects, refining numerical models, strengthening institutional frameworks, and conducting large-scale feasibility analysis positioning Iran to harness its considerable subsurface energy resources sustainably [17].

 Confined and free aquifer

Since the equation of the free aquifer is in the form , it can be written as the equation by changing the variable , which is the same as the Laplace equation or the governing equation  of the grid method. The difference between the free aquifer solution method and the confined aquifer is that in the free aquifer, one of the boundaries is the free water surface, which itself is unknown. In other words, one of the boundaries of this type of aquifer is one of the outputs of the problem, not the input. This type of problem, known as free surface boundary, can be solved in several ways, two of which are explained below.

1- In this method, the free water surface is first drawn hypothetically. By doing this, practically all the boundaries will be known. Then, using this hypothetical boundary, the obtained range is divided and then solved with the desired numerical method. After solving the problem, if the assumption is correct, the pressure of the points on the water surface should all be zero, otherwise, the height of these points should be changed. This change is such that the new height of the points is considered equal to the head of those points. After several iterations, the problem will converge to a free water surface.

2- This method is very similar to the previous method, with the difference that first the boundary head that has a smaller value is replaced with a boundary that has a higher head. By doing this, the obtained range will be larger than the real range and boundaries are formed that will certainly not exist in the real problem. Then, assuming that the new boundaries are of the impenetrable type, the problem is solved. After obtaining the pressure at all nodes, all nodes with negative pressure are removed from the domain and the problem is solved again. This process continues until no node has a negative pressure [18].

 Network Method Modification

Improvement of the obtained equations is possible in two ways. The first modification can be done by changing the way a node is connected to its neighboring nodes. The second modification can be related to the way boundary nodes are modeled. In this section, these two modifications are explained in detail.

 Improvement by increasing the connection of nodes

In a study conducted by (Raoof and Hassanizadeh, 2009), it was pointed out that until then, in most of the research conducted, neighboring nodes were located on horizontal and vertical lines. This means that in a two-dimensional network, a node was connected to at most four nodes and in a three-dimensional network, to six nodes. With this type of connection, due to the absence of diagonal connections, flow in the diagonal direction is not possible at all. However, in this direction, it is possible to have a gradient head. Considering this point, neighboring nodes in the diagonal direction have also been considered in this study. The method of simulating a porous medium using the Square Diagonal PNM (SDPNM) method will be as shown in figure 10.

 To obtain the governing algebraic equation for this method, the continuity equation for an isolated node such as (i) (Figure 3-11) is written as:

Assuming that the desired porous medium is homogeneous and uniform, and considering that the length of the diagonal tubes   is equal to the length of the horizontal or vertical tubes, equation 32 will be written as 33:

 The computational representation of the above equation is shown in figure (12).

Since the number of node connections has increased compared to before, the algebraic equations obtained for the boundary nodes will also change. It should be noted that for modeling the flux passage used in the second and third type conditions, points outside the boundary are considered as before. The method of obtaining these equations according to figure 13 is as follows [19].

 1- Drishle conditions: In this type of boundary, since the value of the function on the boundary is known, there is no need to write an algebraic equation on the boundaries. The equation of a node like (i) (Figure 13) located in the vicinity of this boundary will be as follows:

2- Newman conditions: The governing equation for this boundary for a node like (j) shown in figure 3-13 is equal to:

By equating all three fractions above and writing the equation 3-33, the governing equation for this boundary is written as follows:

3- Rabin's condition: The governing equation for this boundary for a node like (k) shown in figure 3-13 is equal to:

By equating all three fractions above and writing the equation 3-33, the governing equation for this boundary is written as follows:

Improvement using the boundary node modeling method

As mentioned in equations 16, 20 and 24, the obtained algebraic equations were all the same as those obtained from the leading finite difference method. Considering that the central finite difference method has higher accuracy than the leading method, it can be expected that higher accuracy can be achieved with changes in the way the boundary node formulation is formulated. It is observed that the impermeability of the shown layer has been considered in such a way that the water flowing in the pipe (j4) after hitting the horizontal pipe 1-3 does not have the possibility of moving in the vertical direction, and therefore the impermeability condition is satisfied. To satisfy this condition, a virtual pipe that is in the direction of the vertical axis and has exactly the same characteristics as the pipe (j4) can be used. To establish the impermeability condition, it is assumed that the same flow that is present in pipe (j4) also flows in pipe (j2) in the opposite direction, and thus these two flows neutralize each other. This assumption has the advantage over the previous assumption that the momentum of the water flowing in pipe (i4) is neutralized in a gentler manner than before. This assumption can be written as follows:

In other words, the head of nodes 1 and 4 are equal, and no flow is passed between these two nodes. By substituting h₄ for h₂, equation 3-13 becomes:

The conclusion that can be drawn from the above is that for the formulation of boundary nodes when it is necessary to place the flux in a specific direction, the total flux of that node should be written or in other words, the flux through the node should be obtained by the difference of the heads of the two nodes surrounding the node in question [20]. Therefore, the formula for the Newman and Robin conditions will be written as follows according to figure 5.

1- Newman conditions:

And equation 13 will be transformed into equation 42:

2- Robin's conditions:

                                                              (43

By writing equation 13 and substituting equation 43 into it, the following equation is obtained:

The obtained equations (42 and 45) are all equal to the equation obtained through the central finite difference method and therefore have greater accuracy than equations 20 and 24.

Using the above information, the boundary conditions of SDPNM can also be modified by applying the mentioned method. The formulation of the Newman and Robin conditions according to figure 13 will be written as follows:

1- Newman conditions:

Assuming that the expressions of each of the above fractions are equal, the continuity equation (3-33) will be written as equation 3-47:

Which after simplification will become as follows:

 2- Robin's conditions: 

Assuming that the expressions of each of the above fractions are equal, the continuity equation (33) will be written as follows:

Which after simplification will become as follows:

Summary: Core Challenges in Groundwater Energy Storage