Qanat's hydraulic harim determination using the meshless numerical method

Today, due to Iran's annual rainfall of 250 mm yearly, equivalent to one-third of the world average, groundwater resources are becoming increasingly critical (Jafarzadeh et al. 2022). It has created many issues regarding the provision of water for agriculture, drinking and other uses. One of the oldest and most efficient methods of harvesting these resources is the use of qanat or kariz systems (Al-Ghafri et al. 2003), especially in semi-arid and arid areas (Samani et al. 2022). Although modern methods have surpassed them, the qanats still play a major role in the continuity of settlements, even today. This is why they have been used for more than 2,500 years, especially in Iran (Maghrebi et al. 2022).

These systems are naturally and socially sustainable (Khaneiki 2022). Sadly, they have been threatened for the past decades and faced unprecedented challenges for various reasons, including aggressive groundwater pumping which reduces qanat flow rates and its productivity, increased development and urbanization (Esmaeili et al. 2022). Another important reason is the excavation of a new qanat on the hydraulic harim of the other qanats which threatens the older qanat with drought (Salek 2019). Therefore, the qanat systems are under a significant threat of dying out. However, they are vital to water management in Iran (Khosravi et al. 2016).

A key legal concept relating to the qanat system is the law of harim, or borders, which protects the owner of a qanat over a territory surrounding it and prohibits the construction of new mother wells or any other structures within a defined distance of existing qanats (Albert et al. 1998; Hein 2020). On the other hand, harim, a protected area around each well that prevents new wells from affecting the yield of existing wells, is defined for tube wells as well as qanat mother wells, but the harim for qanat mother wells is much broader than that for tube wells (Lightfoot 1996).

The harim law does not specify this distance. In Al-Karaji's discussion, it is apparent that distance depends on the type of hydro-infrastructure, type of qanat, its length and depth, geological condition of the soil, location of the mother well, etc. (Albert et al. 1998). So, harim is more than just protecting ownership rights; it also defines the framework for protecting groundwater resources.

As a result of not paying attention to such critical issues in the last few decades, Iran has encountered a number of serious problems. According to Khosravi et al. (2016), the urbanization of the northern areas of Tehran after the 1960s disrupted the practice and application of harim to older qanats: the deep foundations of many buildings blocked underground water flow.

The importance of harim in qanats has been stressed in many studies over the years. Lightfoot (1996), in order to speak about the shaft and well's placement alongside each other, presented the concept of hydraulic harim and expressed that this area for the mother well is larger than other shafts (Lightfoot 1996). Balali et al. (2009) indicated the effects of modernization and urbanization problems on qanats. Therefore, they expressed the concept of harim and noted that it is usually within 100 and 300 m (Balali et al. 2009). Wessels (2011) brought the regulation about the hydraulic harim and stated the laws about it in order to avoid problems for qanats (Wessels 2011). Aroua (2014) presented the concept of ‘harim-al-ma’ and its importance in order to reduce the health problems of water exploited from qanats. The case study was M'Zab Valley, including four cities covering an area of 300 km² (Aroua 2014). Nabavi (2017) spoke about the harim law. The first benefit of the law is that it preserves established water rights across the country through a concept called ‘harim’ and, a new qanat or well cannot be excavated in that area (Nabavi 2017). Salek (2019), after an explanation of the harim concept, mentioned the necessity of harim determination and stated that this area usually started from wells and shafts and continued until 1.2 km from them (Salek 2019). Mirnezami et al. (2019) presented the definition of harim as the protection zone of qanats and wells according to ancient Islamic tradition (Mirnezami et al. 2019).

Hein (2020) illustrated the parameters affecting the harim of qanats and discussed the challenges that come from building something in the harim area (Hein 2020). Gharios (2020) wrote about harim, its definition and its aim for all water resources (Gharios 2020).

Simonen (2021) indicated that harim is another important aspect of qanats which directly affects inter-communal and intra-disputes (Simonen 2021). Mohan & Kuipally (2021) explained the harim rule for some water resources like wells and springs and stated that this rule identifies the minimum distance between two water point sources (Mohan & Kuipally 2021).

It was reported that the length of Dash-Bolagh qanat in Hamadan plain – an Iranian plain – is around 1,150 m. They also stated that Ali-Khan-Kahrizi qanat has the least length of 13 m. Mirzaei Yazdi et al. (2020) used FEFLOW software in order to simulate groundwater flow. Then with the help of ArcGIS, delineated the hydraulic harim for a standard aquifer. The Sichardt equation was used for the computation of hydraulic harim on that aquifer (Mirzaei Yazdi et al. 2020).

Previously, the qanat's harim area was mainly determined using these formulas, but there were never any standard, user-friendly and more accurate methods for determining this area. Currently, there is an urgent need to find a model that is compatible with all complex conditions of aquifers, has the flexibility to consider different hydrodynamic coefficients near the qanat system, has the ability to take into account the passage of time and is close to real conditions. For this reason, numerical models based on groundwater flow equations were investigated in order to be able to determine the harim of these resources in a more accurate and faster manner to achieve more realistic results.

In this regard, some numerical methods can be used for this aim. Finite difference, and finite element are the traditional numerical methods (Wang & Anderson 1995). Due to their dependency on the meshing stage, some errors exist in their results. However, a new group of numerical methods named the meshless method, presented in the last decade which removes this drawback and releases results with higher accuracy than the traditional methods (Liu & Gu 2005). Meshless local Petrov–Galerkin (MLPG) is one of the numerical methods which can be used for fluid dynamics for this aim (Atluri & Zhu 1998). Mohtashami et al. (2017, 2021) developed MLPG for groundwater simulation in both steady and unsteady states. They found their model more accurate than the finite difference method (Mohtashami et al. 2017, 2022a).

Delineation of hydraulic harim using the numerical method is not investigated in the previous studies. In this study, the hydraulic harim of a qanat system is depicted using the MLPG numerical method. This was done for the first time during this research; after introducing the qanat system to the developed model, the hydraulic harim will be delineated for the collection section. This developed model is applied to two standard and real field study aquifers.

Figure 2 illustrates the structure of a qanat system in an unconfined aquifer. The aquifer in the study area is saturated. The first well is the mother well, which has the maximum flow rate. Other wells are located in both the collection and conveyance sections.

During the generation of the cone of depression, the shaft in the collection section has a hydraulic harim. Our purpose in this study is to define this area.

Between 1955 and 2021, the average yearly rainfall was 145.9 mm. The average amount of rainfall in the plains was estimated at 166.3 and 173.4 mm/year in the mountainous area. Based on regional classifications, this area is categorized as arid, and its average yearly temperature is 16.6 °C. The amount of summation of evaporation in the heights of the area is calculated to be 2,234.26 and 4,607.19 mm in the plains.

MLPG is one of the weak-form meshless methods that enforce governing equations in nodes instead of elements. In 1998, Atluri and Zhu introduced it for the first time to solve the potential equation. Their emphasis was on the effectiveness and accuracy of this method in solving equations, especially those relating to fluid mechanics. MLPG enforces the governing equation on each node and then assembles the equations in one matrix in order to find the unknown coefficient (Atluri & Zhu 1998).

This method for groundwater flow simulation in the Birjand aquifer was first developed by Mohtashami et al. (2017). They calibrated and validated their model and found their model to be more accurate than traditional methods.

Here, H is groundwater head, k denotes hydraulic conductivity coefficient, Q shows the discharge (+) or recharge (–) rate, S is the storage coefficient and q indicates the precipitation and evaporation.

In Equation (3), W and represent weight and shape functions. Here, K, F and U represent the stiffness matrix, force body vector and the unknown variables vector.

Steps for delineating hydraulic harim are presented in the following figure.

Details like the location of a mother well, shafts and the gallery are entered into the model; also, with the usage of the Darcy equation, the flow rate of a mother well and shafts are calculated and assigned to the model. After all, the groundwater head is computed. With the initial groundwater head, the drawdown value presented by ‘s’ is computed.

At the final stage, the nodes placed in the area with the determined drawdown are illustrated. Based on the drawdown values of these nodes, contour lines will be drawn and hydraulic harim will be depicted in the aquifer.

The MLPG numerical method is going to be used to delineate hydraulic harim within the collection section of qanats. In this way, the drawdown value for the aquifer which is caused by the extraction of qanats is computed and based on that, an area for different drawdowns is plotted.

In order to find the drawdown values, the groundwater equation must be solved. To this end, MLPG is used. Solving an equation with MLPG has three steps: pre-processing, processing and post-processing. In the pre-processing stage, in the aquifer and its boundary, defined numbers of nodes are uniformly distributed. Each node has five layers of information containing: hydraulic conductivity coefficient, specific yield, initial groundwater head, boundary condition and recharge or discharge value (Mohtashami et al. 2017). After these layers are assigned, the processing stage begins, and the stiffness and force body matrices are calculated. During the last stage, groundwater head and after that drawdown values are computed and contour lines are drawn. According to these lines, this area can be achieved.

It is noteworthy to state that MLPG was a calibrated and verified model for groundwater flow simulation based on previous articles by the authors. This model in both unconfined and confined conditions, steady and unsteady forms, was calibrated and validated and applied to standard and real aquifers. Results are achieved by Mohtashami et al. (2017, 2021). In these studies, groundwater head which is presented by MLPG was compared to MODFLOW results (which is based on the finite difference method) and other meshless methods. MLPG had higher accuracy than other numerical methods. This was shown with three evaluation indices. Root mean square error (RMSE), mean absolute error (MAE) and mean error (ME) were these indices. Table 1 shows the performance of MLPG for standard aquifer and Table 2 is for Birjand unconfined aquifer as a real test case according to the mentioned evaluation indices.

In order to draw the hydraulic harim of a qanat system, we first move through the standard aquifer and consider that the standard aquifer has one qanat system and is unconfined. This system includes one mother well, three shafts in the collection section and three shafts in the conveyance part. The volume of water drained from the mother well and shafts placed in the collection section is computed with the Darcy formula. In this way, Table 3 details the location and flow rate of the mother well and collection shafts. There is no drained flow rate in shafts placed in the conveyance section because the groundwater table is at a lower depth than these shafts.

Simulations were conducted and a model was run for this standard aquifer. Hydraulic harim was depicted for three time sections, 30, 60 and 120 days after exploitation from the qanat.

Hydraulic harim for the collection section of qanat is depicted. The area was drawn for 30 days of extraction from qanat. In the collection section, the mother well and two shafts have remarkably large areas of harim, but at the third shaft, no harim is visible.

It should be noted that the drawdown values for the mother well and the two shafts are 4.5, 4.1 and 4 m, respectively. We delineated this area after 30 days extraction from qanat.

According to this figure, there is no hydraulic harim for the third shaft. This is because of the short time of extraction and the low level of water table on the node.

The following figure presents the hydraulic harim for 60 days extraction.

The area of the harim is 42,500 m², the width and length are 200 and 400 m, respectively. Both of these values have increased over the previous figure.

Figure 9 shows hydraulic harim after 120 days of exploitation.

You can see that in this case, harim is created around the third shaft, but it is too small in comparison to the others. This is the final harim, because after 109 days, the aquifer reaches its steady state and the cone of depression does not change. Therefore, this area remains constant after 109 days, even after 120 days.

This harim has the most area, length and width values in this situation, also the mother well and the three shafts have the greatest drawdown value. However as mentioned previously, the drawdown after 109 days remains constant.

As a result, Table 4 shows the area, width and length of hydraulic harim for standard aquifers in different time sections.

The length of this qanat is 2,500 m. Its flow rate range is [15–45] L/s, the depth of the mother well is 75 m and the number of shafts is 150. Water extracted from this qanat is used for agricultural irrigation. The path of this qanat is depicted in Figure 11.

Figure 11 shows the nodes scattered uniformly around the qanat outlet, as well as its path, tunnel and mother well. In order to draw the hydraulic harim of Haji-Abad qanat, first, the conceptual model for the Birjand aquifer was created. Then, the qanat package was entered into the MLPG simulation model. The location of the mother well, collection shafts and outlet were defined in the meshless method and is assigned to the nearest node. They are represented in blue, green and red symbols in Figure 11. The waterway of qanat is also shown in the same figure with a blue line. It is noticeable that the flow is heading toward the west. Figure 11 also shows the outlet of Haji-Abad qanat.

Simulation of groundwater flow was done for 11 years with a monthly time step (2011–2022). Then, the hydraulic harim for the chosen qanat in different time steps was computed and depicted. The simulation was run under the assumption that exploitation started in 2011 and would continue until 2022. Time steps were 3, 6 and 11 years.

Some points are derived from Figure 12: 1 – hydraulic harim area gets larger due to time passing. In other words, the number of yellow nodes is higher for harim of 11 years than 6 years, and for harim of 6 years than 3 years. 2 – hydraulic harim near the mother well is more extended than the other collection shafts. This is due to the higher groundwater table around the mother well which causes more flow rate than others.

3 – Extension of hydraulic harim is toward the regions with higher conductivity coefficient and specific yield. According to the values of hydraulic conductivity and specific yields, the role of hydraulic conductivity coefficient is more significant than the specific yield.

The hydraulic harim extended about 500 m from the mother well in 2040. The extension of this area is near the regions that have higher hydraulic conductivity coefficients such as the upper reaches of qanat.

In this way, the calculation of the harim area is important for preventing some probable problems faced by many of our qanats. The hydraulic harim principle prohibits drilling a water well within a specified distance of another water project, which is the main traditional law to prevent qanats from drying out. And this is achieved just with the help of simulation techniques, such as the accurate one. As such, the findings of this study can be applied to all types of aquifers and qanats, which are subject to varying conditions, as well as to arid aquifers that harvest water using this traditional method.

As a result, a comparison is carried out between the proposed model and other empirical equations like the Sichardt equation and drainage channels for hydraulic harim depiction of a qanat system. Table 6 presents the advantages and disadvantages of the proposed model in comparison to empirical equations.

The proposed model is applicable to both steady and unsteady states. It considers different values of hydraulic conductivity coefficient and specific yield near the qanat system, however, in empirical formulas, only one value can be entered into the equation. Also, as the engaged model, first, simulates groundwater flow and then depicts the hydraulic harim is more accurate than the empirical equations.

The only disadvantage of the proposed model is the high computational cost. Due to the heavy computation of this model, the runtime is much higher than empirical relations. The reason is solving a set of algebraic equations in MATLAB software.

In arid and semi-arid regions, qanat systems are used to collect groundwater. Qanats that drain water from aquifers have borders around them named hydraulic harim. Investigating this area in order to protect qanat from changes in flow rates and quality risks is critical and necessary for all qanat systems. Hence, in this study, we aimed to achieve the hydraulic harim of qanat by using a precise numerical method named meshless local Petrov–Galerkin. This is the first time it has been done. For this purpose, two aquifers, one standard and one real field aquifer, are chosen. Then, the qanat package has been added to the MLPG groundwater model. Finally, hydraulic harim is computed for the qanat. In the standard aquifer which has one mother well, three shafts in the collection section and three shafts in the conveyance section, harim is depicted. Hydraulic harim is drawn for 30, 60 and 120 days after exploitation. The hydraulic harim shapes for the three time periods are different. The length and width of harim for 60 days are greater than those for 30 days. Harim for 120 days is also larger than both 60 and 30 days. This is due to the high drawdown value for 120 days. Harim retains its constant shape after 120 days due to the equilibrium condition that governs the aquifer. Haji-Abad, one of the most important qanats in Birjand, is chosen for the real test case. Its mother well has a depth of 75 m. It has 10 shafts in the collection section. The hydraulic harim for this qanat is determined after 3, 6 and 11 years after extraction. Results show that the harim has expanded over time because of more extraction and drawdown. For instance, the length of hydraulic harim after 3, 6 and 11 years are 400, 700 and 900 m, respectively. The width of harim for these years varies from 200 to 300 m as well. This area is also larger for the mother well in comparison to other collection shafts. This is due to the high amount of water drained from the mother well. Results indicate the harim's development towards regions with a higher hydraulic conductivity coefficient. We also predict the hydraulic harim for this qanat 18 years down the track and find some buildings will be built in that area, and the qanat may face some problems. The length of harim this year is 1,700 m.

The authors want to thank the regional water company of South-Khorasan province according to providing the sufficient data.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.