Evolution of low-karstified rock-blocks and their influence on reservoir leakage: a modelling perspective

Hydraulic structures such as dams and reservoirs pose significant construction challenges in karst areas due to severe and costly leakage issues. In this study, we apply a numerical model to test the hypothesis that karst aquifers in water divide areas may contain intrinsically low-karstified rock-blocks (LKB), which form due to the specific evolution of unconfined aquifer with recharge distributed to the water table. We develop, test, and apply a model of flow, transport, and dissolution in a 2D fracture network with a fluctuating water table. The model's structure and boundary conditions are based on the conceptualization of the Luojiaao (China) interfluve aquifer. First, we simulate the evolution of an unconfined network, representing the interfluve, up to a stage resembling the present conditions at Luojiaao. We then analyse leakage through the evolved aquifer from a reservoir at different water levels and simulate further aquifer evolution under reservoir conditions. Our results demonstrate the formation of the LKB and highlight its role in mitigating leakage.

1 Introduction

Karst areas occupy approximately 14 % of Earth's ice-free land surface. In these regions, the chemical erosion of rock by surface water and groundwater – known as karstification – produces distinctive surface and subsurface features. Karst aquifers are among the most productive freshwater resources, supplying drinking water to around 20 % of the world's population – and more than 50 % in some countries (Chen et al., 2017; Ford and Williams, 2007). These aquifers are characterized by networks of solution conduits formed within initially fractured and porous rock. The development of these conduits, a process known as speleogenesis, is complex and spans nearly the entire period during which the rock is exposed to groundwater circulation (Dreybrodt, 1996; Yuan et al., 1993).

Karst aquifers can be conceptualized as hierarchically organized conduit systems embedded within a fractured porous matrix. These conduits serve as high-transmissivity flow pathways, channelling water from points of recharge to karst springs. Karstified rocks are also widespread in young orogenic belts, where the high-relief terrain provides favourable elevation differences that facilitate the construction of hydraulic structures such as dams and reservoirs (Milanović, 2000; Yuan et al., 1993; Shen et al., 1997).

Figure 1 Simplified cross section of an unconfined aquifer in water divide region. (a) Low flow zone in the water divide area. (b) Karstification increases permeability and progressively lowers the water table over time, but leaves the water divide zone less karstified (LKB). (c) In case of reservoir construction, the LKB can effectively mitigate leakage across the aquifer (Modified from Zhang et al., 2018).

However, the construction hydraulic structures in karst areas pose significant engineering challenges. The position, size, and connectivity of existing networks of solution conduits are extremely difficult to predict or detect. As a result, the extent of initial leakage from a reservoir is also hard to estimate (Milanović, 2000, 2018). Furthermore, the high hydraulic gradients induced by such structures can accelerate the formation of new conductive pathways, potentially leading to a rapid increase in leakage – sometimes reaching intolerable levels within the operational lifespan of the structure (Dreybrodt, 1996; Gabrovšek and Dreybrodt, 2001, 2010; Romanov et al., 2003, 2007).

Basic flow solutions in unconfined porous aquifers with constant recharge as shown in Fig. 1a, suggest a relatively stagnant flow zone in water divide area (Rhoades and Sinacori, 1941; Tóth, 1962; Liang et al., 2010). Studies on porous aquifers have revealed a zone of low permeability in divide areas, characterized by slower flow velocities (Kang et al., 2003; Wang et al., 2011; Jiang et al., 2012; Nogues et al., 2013).This also applies to fractured aquifer and can have important implication for karstification, where low flow zone may also result in less karstified zone, as shown in Fig. 1b. Such low-karstified rock-blocks (LKB) have been recognised in water divide regions of real karst aquifers and proven to be effective in mitigation of leakage from reservoirs in karst areas (Yuan et al., 1993; Milanović, 2000; Xu and Yan, 2004). The LKB means that there is a relatively low permeable region formed in the middle of karst aquifer during the natural karst evolution. The concept is illustrated in Fig. 1, showing a low karstified block formed in the water divide region during natural karstification (Fig. 1a and b). As depicted in Fig. 1c, this zone can effectively mitigate leakage from man-made reservoirs.

Figure 2 The flowchart of the whole modelling process. At each time step, the new position of the water table is determined through an iterative process. Subsequently, the coupled flow, dissolution, and transport equations are solved at both the individual fracture scale and the fracture network scale. Based on these results, the fracture apertures are updated accordingly. The modified fracture network then serves as the basis for calculations in the following time step.

Due to the intricate nature of karst leakage, there have been limited reports on quantitative modelling studies of this phenomenon (Milanović, 2000, 2018; Smerdon et al., 2005). Modelling karst aquifer leakage is a challenging endeavour, often deemed impractical despite its importance for engineering planning. In engineering construction, predictions of reservoir leakage primarily rely on geological surveys and field experiments. For example, one study analysed extensive observational data and predicted that leakage would remain minimal (Yuan et al., 2002; Shen et al., 1997).

To model evolution of karst aquifers, it is essential to simulate flow, dissolution, and transport processes within soluble fractured media. When modelling an unconfined aquifer that includes a water divide, the problem involves a free surface, requiring the determination of the water table position at each time step of the model. Some studies have addressed this by treating discrete fractures as equivalent continuum media (Desai and Li, 1983). However, such models are unable to capture the processes of fracture dissolution and widening. Several studies have reported on modelling discrete fracture flow with free surfaces. For instance, Wang (1993) employed an initial flux method to locate the free surface within the fracture system, while Jing et al. (2001) used a discontinuous deformation analysis (DDA) method for modelling fracture flow. These intuitive methods often encounter numerical instability due to their lack of a solid mathematical foundation. The variational inequality method, which introduces the continuous Heaviside function in the finite element method, has proven effective in modelling the free surface in discrete fractures (Zheng et al., 2005; Jiang et al., 2013). But these models do not consider situations with recharge to the water table.

Several models have been used to simulate the evolution of unconfined karst aquifers with constant recharge at the water table (Clemens et al., 1999; Gabrovšek and Dreybrodt, 2001, 2010, 2021; Kaufmann, 2003; Dreybrodt et al., 2005). As demonstrated in these works, the evolution of unconfined karst aquifers is most pronounced near the water table. These models do not explicitly account for the details of flow and dissolution in the vadose zone – the region between the surface and the water table – but instead assume gravitational flow that reaches a specific level of chemical saturation with respect to calcite. Nevertheless, within the phreatic zone, water is typically least saturated with respect to calcite near the water table, and thus possesses the highest potential for dissolution (often referred as aggressivity) in this region. The numerical methods employed in these studies primarily focus on the evolution of relatively homogeneous fracture networks, rather than on the randomly distributed fractures typical of natural karst aquifers.

In this work, we introduce a numerical model to study the evolution of karst aquifer under distributed (rainfall) recharge. The algorithm is built upon an earlier study by Yu et al. (1999). We propose an intuitive method that uses iterative search for the water table and couples flow, transport and dissolution to simulate evolution of discrete fracture network. We first simulate the natural evolution of the aquifer under rainfall recharge up to selected stage of maturity. Then, the obtained aquifer is subjected to boundary conditions imposed by the construction of the reservoir. We analyse how the LKB has evolved and assess its effectiveness in mitigating leakage from the reservoir.

2 Modelling evolution of unconfined discrete fracture network

Models of fracture networks' karstification necessitate the integration of flow, dissolution, and transport processes within the system, while adhering to specified hydraulic and hydrochemical boundary conditions. In our research, we simulate the evolution of a stochastic fracture network under constant recharge and fixed head boundary conditions. The computational procedure is shown in flowchart of Fig. 2 and outlined as follows:

1. Generate random fracture network and boundary conditions;

2. Locate the water table and compute the flow within the network;

3. Calculate dissolution and solute transport;

4. Update fracture apertures: compute the changes in the width of each fracture within the specified time step;

5. Iterate or terminate: proceed to the next time step and repeat steps 2 to 4, or terminate the simulation as required.

Figure 3 Conceptual presentation of flow calculation in the fracture network. (a) Boundary conditions for confined flow calculation performed at every iterative step. (b) Node-by-node iteration: testing wetting of a dry node across a layer. (c) Layer-by-layer iteration.

2.1 Flow in fractures

2.1.1 Fracture network generation

The initial fracture network is created through Monte Carlo stochastic simulation, as described by Yu et al. (1999). Fracture parameters, including centre location, length, aperture, and direction, are randomly generated using Gaussian and lognormal distributions. Subsequently, the algorithm searches for fracture intersections (nodes) and excludes all dead end segments, which do not contribute to flow. In this work the fractures are uniform plan parallel planes, defined by length and aperture.

2.1.2 Flow in one fracture

When a fracture with length L_ij and aperture b_ij connects nodes i and j with hydraulic heads H_i and H_j, respectively, the flow rate q_ij within the fracture conforms to the cubic law for laminar flow, whereas for turbulent flow, the Lomize equation (Lomize, 1951) is used to relate the flow rate to the hydraulic gradient J_ij.

$$\begin{array}{}\text{(1)}& {\displaystyle }{q}_{ij}={\displaystyle \frac{\mathrm{g}{b}_{ij}^{\mathrm{3}}}{\mathrm{12}v}}J,\text{for}\mathit{Re}<\mathrm{2020}\text{(2)}& {\displaystyle }{q}_{ij}=\mathrm{4.7}\sqrt[\mathrm{7}]{{\displaystyle \frac{g^{\mathrm{4}}}{\mathit{\upsilon }}}{b}_{ij}^{\mathrm{12}}{J}^{\mathrm{4}}},\text{for}\mathit{Re}>\mathrm{2020}\text{(3)}& {\displaystyle }{J}_{ij}={\displaystyle \frac{H_i-H_j}{L_{ij}}}\end{array}$$

where g is the gravitational acceleration, v is the kinematic viscosity of water, and Re is Reynolds number. The transition from the laminar to the turbulent flow regime is here introduced at the point where turbulent frictional losses exceed those associated with laminar flow. The factor in Lomize equation (Eq. 2) is set so that transition occurs at a Reynolds number around 2020. This simplification omits the transitional flow zone, resulting in a smooth increase in flow rates as turbulence sets in Su et al. (1994) and Finenko and Konietzky (2024). However, in simulations presented in this work, flow always remains laminar.

2.2 Solution for water table in random fractures under distributed uniform recharge conditions

Each evolution time step requires the calculation of flow, dissolution, solute transport, and the corresponding changes in fracture aperture. Accurate flow calculation depends on determining the position of the water table at each time step. This is achieved through an iterative process, in which the water table position is updated until specific convergence criteria are met, as described in Sect. 2.2.2. The procedure is illustrated in Fig. 3. At each iteration, the current approximation of the water table defines a set of boundary conditions (Fig. 3a): prescribed recharge at water table nodes, seepage face (boundary nodes, where hydraulic head equals the elevation), and constant-head conditions elsewhere. The flow solver (confined flow solution) is then invoked to compute flow rates and hydraulic heads at all nodes. Water table nodes are subsequently searched node by node within each layer (Fig. 3b, node-by-node iteration) and layer-by-layer throughout the domain (Fig. 3c, layer-by-layer iteration).

2.2.1 Confined flow solution

During each node-by-node iteration, an approximate water table is used to classify boundary nodes as having either a specified head, a prescribed recharge, or as part of the seepage face. At every node, a flow balance must be satisfied:

$$\begin{array}{}\text{(4)}& \sum q_{ij}+\sum q_w=\mathrm{0}\end{array}$$

where the first part is the sum of flow from/to all inner (non-boundary) nodes (see Eqs. 1–2) and the second part is the sum of direct recharge q_w from the vadose zone, applied to the water table nodes.

Assuming that the equations are non-linear, we use the Newton–Raphson iteration to solve the system of equations, with convergence criteria as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }{H}_{\max }\le H_{\mathrm{tol}},\text{for}{H}_{\max }=max\left(\right|{\mathit{H}}_{\mathrm{error}}\left|\right)\text{(6)}& {\displaystyle }{Q}_{\mathrm{erro}}\le Q_{\mathrm{tol}},\text{for}{Q}_{\mathrm{erro}}=\sum q\end{array}$$

where H_error is the water head error vector between the n and n+1 iterations. H is the water head vector. H_max is the maximum water head difference between two iterations and H_tol is the allowable error. Q_erro is the sum of water flow budget of all inner nodes, and Q_tol is the allowable error.

The difference in the water heads between two iterations and the water budget is checked to validate the flow results. When H_max is less than the allowable error H_tol and the water budget Q_erro is less than the allowable error Q_tol, the iteration converges. The C$++$ open source linear algebra package Armadillo (Sanderson and Curtin, 2016, 2019) was employed to perform matrix computation.

Figure 4 Location and setting of the Geheyan Reservoir, situated on the upstream side of a bend in the Qingjiang River, Changyang County, Hubei Province, China. (a) The location of Geheyan Reservoir. (b) Simplified map of the reservoir area and Luojiaao interfluve aquifer (marked by a dashed rectangle in a). (c) Cross-section through the Luojiaao interfluve aquifer along the line AB in (b). The dashed line shows position of the water table before the reservoir construction.

2.2.2 The node-by-node and layer-by-layer iterations: search for the water table

One of the challenges in modelling unconfined aquifers with distributed recharge lies in accurately identifying the water table. In the presence of dissolution and increasing permeability of fractured media, the water table changes in time also when recharge is constant. We use a robust “intuitive” algorithm for computing the position of a water table as shown in Fig. 3. The algorithm employs a two-level iterative strategy: (1) layer-by-layer iterations and (2) node-by-node iterations within each layer. When the algorithm solving for the confined water head is included, the method effectively operates as a three-level iterative strategy, ensuring robustness and accuracy in resolving wetting and free surface dynamics. The algorithm is applied for all time-steps of a simulation.

Step 0:

Initially, the water table is positioned at the lowest constant head boundary in the domain (Fig. 3a). The nodes of the water table receive direct recharge from the vadose zone. The heads at all nodes in the phreatic part are calculated by above described confined flow solution.

Step 1:

Dry nodes, directly connected to water table nodes are identified. For each pair, the difference between the hydraulic head of the nearest water table node $H_i^{\mathrm{WT}}$ and the elevation of the dry node $z_j^D$ is calculated as ${\mathrm{\Delta }}_{ij}=H_i^{\mathrm{WT}}-z_j^D$. These differences are then sorted into an ordered sequence.

Step 2:

Following the sequence each dry node is changed to a water table node and the heads are recalculated by confined flow solver, assuming new boundary conditions with constant input to the new WT node (Fig. 3b). If the head at the node is higher than its elevation, and no other wet node becomes dry, the transition of the node from dry to wet is kept. Otherwise, the node remains dry and the procedure continues with the next node in the sorted sequence.

Step 3:

Repeat Steps 1 and 2 layer by layer until no dry node becomes wet (Fig. 3c).

The numerical model was verified successfully against a MODFLOW model (Harbaugh et al., 2000) and an analytical model using the Dupuit assumption. Please see the Supplements for details.

Table 1 Statistical parameters of the original fractures.

Figure 5 The initial aquifer (a), and aquifer after 10 ka (thousand years) (b) and 20 ka (c) of evolution. The widths of the lines present the aperture, colours the head and arrows the flow direction. To facilitate the discussion the aquifer is divided into 3 primary vertical sections and each of them into 5 horizontal subsections, as shown in panels (a) and (d).

Figure 6 Evolution of fracture apertures under natural karstification at 0 ka (thousand years) (a), at 10 ka (b) and at 20 ka (c).

2.3 Dissolution widening of fractures

The dissolution rates of limestone under both laminar and turbulent flow conditions have been the focus of extensive research (Dreybrodt, 1990, 1996; Liu and Dreybrodt, 1997; Eisenlohr et al., 1999; Dreybrodt and Kaufmann, 2007). These rates are governed by a complex interplay of surface reactions, transport mechanisms, and CO₂ conversion processes (Dijk et al., 2002; Detwiler and Rajaram, 2007). In this study, we adopt the rate equations developed by Dreybrodt (1996):

$$\begin{array}{}\text{(7)}& {\displaystyle }F\left(c\right)=k_{\mathrm{1}}\left(\mathrm{1}-{\displaystyle \frac{C}{C_{\mathrm{eq}}}}\right),\left(C<C_{\mathrm{s}}\right)\text{(8)}& {\displaystyle }F\left(c\right)=k_n{\left(\mathrm{1}-{\displaystyle \frac{C}{C_{\mathrm{eq}}}}\right)}^n,\left(C>C_{\mathrm{s}},n>\mathrm{1}\right)\end{array}$$

where k₁ and k_n are rate constants (in mol cm⁻² s⁻¹), C represents the concentration of Ca²⁺ ions (in mol L⁻¹), and C_eq denotes the equilibrium concentration in the H₂O-CO₂-CaCO₃ system. The C_eq is determined by the concentrations of calcium and dissolved CO₂ in the solution at the water table. These concentrations result from complex flow and dissolution processes occurring from the surface down to the water table – including flow and reactions within the soil, epikarst, and vadose zone – which are beyond the scope of this study. The reaction follows a linear rate law up to the switch concentration C_s, transitioning to a nonlinear rate law between C_s and C_eq. Here we assume that rates are entirely controlled by surface reaction, and ignore the concentration gradient perpendicular to flow. This approximation is valid for situations where the solutions are close to equilibrium (Gabrovšek, 2000), which is mostly the case in the presented scenarios (Yuan et al., 2002; Wan et al., 1999). The reaction ordern and the switch concentration C_s are experimentally determined and influenced by the impurity content in the limestone (Eisenlohr et al., 1999). Experimental studies have reported reaction orders ranging from 3 to 11. For this study, we use n=4, C_s=0.9 C_eq, $k_{\mathrm{1}}=\mathrm{4}\times {\mathrm{10}}^{-\mathrm{11}}$ mol cm⁻² s⁻¹, and $k_{\mathrm{4}}=\mathrm{4}\times {\mathrm{10}}^{-\mathrm{8}}$ mol cm⁻² s⁻¹ (Gabrovšek and Dreybrodt, 2000; Dreybrodt et al., 2005). Studies on the dissolution kinetics of calcium carbonate minerals under turbulent flow (Buhmann and Dreybrodt, 1985a, b; Liu and Dreybrodt, 1997) indicate that the rate constant in turbulent conditions is an order of magnitude greater than that under laminar flow.

The change in concentration ΔC along an individual fracture can be determined by applying the principle of mass conservation within a water parcel passing along the fracture. This can be analytically calculated for a plane parallel fracture, as shown in the Supplements. The change in concentration along the fracture is converted to mass removed from the fracture walls and the change of aperture during time Δt:

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }b={\displaystyle \frac{\mathrm{\Delta }Cq{M}_{{\mathrm{CaCO}}_{\mathrm{3}}}}{\mathit{\rho }L}}\mathrm{\Delta }t\end{array}$$

where L is the fracture length, ρ is the density of calcite (2.5 g cm⁻³), $M_{{\mathrm{CaCO}}_{\mathrm{3}}}$ is the molar mass of calcite (100 g mol⁻¹), and q is the flow rate.

To calculate the change in concentration along a fracture, the calcium concentration at the input node must be known. This input concentration is determined as the flow-weighted average of all inflows to that node. To ensure that all required concentrations are available when needed, the computation proceeds in a cascading sequence, starting from nodes with the highest hydraulic heads and moving downstream (see the Supplement). A uniform widening rate is assumed along each fracture. Once the fracture apertures in the entire phreatic zone are updated to reflect dissolution widening within a time step, the evolution time step is completed.

Figure 7 The box-whisker plot of fracture aperture under natural evolution for 20 thousand years in three primary sections (see Fig. 5).

3 Karstification of an unconfined interfluve aquifer with uniform recharge

As an example of an aquifer, where the concept of LKB formation can be applied is Luojiaao interfluve aquifer, located between the two reaches of the Qingjiang river meander. Hubei Province, China (Fig. 4a and b). On the western side of the interfluve, the Geheyan Reservoir was constructed, which is one of three karst reservoirs along the river (Xu and Yan, 2004). The water level in the reservoir is above the original water divide, which was determined by boreholes (Fig. 4c). However, the leakage and its trend could not be determined with field methods.

Previous studies (Yuan et al., 2002; Shen et al., 1997; Wan et al., 1999) revealed cave development near both riverbeds, with an absence of caves in the water divide region. Although we lack knowledge on the possible previous karstification stages at this particular site, we expect that the water divide in the interfluve is a result of karstification under rainfall recharge, and that the LKB may play a crucial role in reservoir leakage.

3.1 Construction of conceptual and numerical model of an initial aquifer

The conceptual model, translated into a numerical one, is based on a simplified interpretation of field mapping data. The modelled domain represents a two-dimensional cross-section of a fractured limestone aquifer, 2000 m in length and 500 m in height. The initial conductive structures are gently dipping bedding planes (8–10°) and a set of near-vertical fractures with an average dip of 85°, both generated randomly using the statistical parameters listed in Table 1. Their initial apertures are uniformly distributed, with a mean value of 0.01 cm. After identifying fracture intersections (nodes) and removing all dead-end segments, the resulting fracture network consists of 4073 intersection nodes and 7401 fracture segments, with segment lengths ranging from 0.01 to 143 m. The average segment length is 15.74 m.

Boundary conditions include a constant hydraulic head of 100 m applied to the left and right sides of the domain (from elevations −100 to 100 m), and seepage boundaries between 100 and 400 m elevation. A uniform infiltration rate of 200 mm yr⁻¹ is applied to the water table nodes.

Reconstructing hydrogeochemical boundary conditions at the water table is inherently challenging. Based on measurements of soil CO₂ and Ca²⁺ concentrations from three boreholes and two springs, we assume a closed-system equilibrium concentration (C_eq) of 2 mmol L⁻¹ (Yuan et al., 2002; Wan et al., 1999). The saturation ratio of infiltrating water is assumed to be $C_{\mathrm{in}}/C_{\mathrm{eq}}=\mathrm{0.92}$. The total simulation time is set to 100 000 years, with model convergence evaluated at each time step.

Figure 8 The relative change in the average aperture $b/b_{\mathrm{0}}$ and the change in discharge rate q (flow out) in different subsections.

3.2 Model of karstification under distributed uniform recharge

The modelling results are shown in Fig. 5. We observe the expected drop of water table due to increasing transmissivity of the aquifer. The groundwater divide remained close to the middle of the domain, whilst the elevation of the highest wet node, which represents the divide ridge, descended from 350 to 230 m at 10 ka and to 185 m at 20 ka.

As expected, there is a fringe of high dissolution rates close to the water table. Below the water table, the fracture widening is relatively slow. As the water table gradually decreases, the fringe of high dissolution rates migrates vertically through the aquifer and increases the hydraulic conductivity.

Figure 9 Changes and increasing trends of the equivalent conductivity K in the different subsections during 20 ka natural evolution process.

3.3 Changes in fracture aperture and discharge during karstification

The fracture aperture distributions at three evolution time steps are shown in Fig. 6. To facilitate the discussion, we have divided the aquifer into 15 subsections along three vertical primary sections, as shown in Fig. 5a. The variation of aperture statistics in three vertical sections is shown in Fig. 7. The water table with dissolution fringe mainly descends through the upper section, which experience the evident change in aperture (Gabrovšek and Dreybrodt, 2001; Hubinger and Birk, 2011).

The changes of apertures and discharge rates in several representative subsections are shown in Fig. 8. Generally, the change of apertures drops along vertical directions and is higher near the left and right boundary. In the lower section, the change of aperture is minimal (A13).

The discharge rate in each subsection was computed by checking the fractures that cross the boundary and flow into the subsection. The lateral aquifer width was assumed to be 100 m. The A13 subsection has the lowest discharge, which changes the least over time. The discharge in the A8 subsection changed more slowly than that in other parts of the middle primary section. The A3 subsection, in the middle part of the upper primary section, exhibited a much greater leakage trend than did the lower parts. From the aperture and discharge analysis, we concluded that a low-karstified area exists mainly in the middle and lower primary sections.

3.4 Changes in conductivity K and formation of the LKB

To further elucidate the formation of an LKB, we calculated vertical and horizontal hydraulic conductivity in each subsection. To do this, we “cut out” a section at certain stage of evolution, and calculated confined flow rate for a given head difference along horizontal and vertical direction, respectively. Thus horizontal and vertical equivalent K values, and the changes and trends in the representative subsection were computed (Fig. 9). Notably, the fractures below 50 m in the lower subsections (A11 to A15), remained unevolved. The horizontal K varied between 0.005 and 0.008 m d⁻¹ in the different subsections, and the vertical K varied between 0.001 and 0.002 m d⁻¹.

The middle subsections A7, A8, and A9 experienced low level of karstification, in particular, subsection A8 (less than 0.01 m d⁻¹), which combined with the lower section forms an LKB, which could effectively prevent horizontal flow. In the A6 and A10 subsections, the K increased by more than one order of magnitude.

The upper primary section was also tested to determine whether it included a low permeable part. Note that here only phreatic karst process is considered, therefore most fractures in subsection A1 did not experience widening. The water table descended through subsection A3, where the conductivity trends increased in both the horizontal and vertical directions compared with those in the other subsections. Although the upper primary area does not exhibit a low conductivity, its underlying low permeable rock-blocks have already been recognized as an LKB.

Figure 10 Flow fields within the aquifer, which has evolved for 20 000 years under natural conditions, after the application of reservoir boundary conditions on the left side at: (a) 200, (b) 250, and (c) 300 m.

Figure 11 Leakage as a function of reservoir level for aquifer after 20, 30 and 40 ka of evolution. The lateral width is 100 m.

4 The influence of LKB on reservoir leakage

To verify the relevance of LKB formation on reservoir leakage, as in the case of Luojiaao interfluve, we now apply reservoir boundary conditions to an evolved aquifer. The borehole data revealed that the natural water divide level is 50 to 80 m above the river, which roughly corresponds to the situation in our synthetic aquifer after 20 ka of evolution.

4.1 Reservoir leakage with various reservoir levels

First, we test the leakage through the aquifer at different reservoir water by applying constant head of 200, 250 and 300 m, on the left boundary (presenting reservoir) of the aquifer which has evolved for 20 ka. The flow fields are shown in Fig. 10. In case of 200 m, the groundwater divide still exists, although it is shifted towards the left boundary. However, the water still flows from the aquifer to the reservoir. When the water is raised to 250 m, as shown in Fig. 10b, the groundwater divide moves almost to the left boundary, so that the flow to the reservoir is only local (upper left), however, deeper the flow is directed from the reservoir to the aquifer. At 300 m, as shown in Fig. 10c, the flow direction is almost entirely from the reservoir to the aquifer.

Figure 11 shows the reservoir leakage at different water levels through the aquifer at different stage of evolution, assuming the lateral width of the aquifer is 100 m. For the aquifer after 20 ka of evolution under rainfall recharge, no leakage is observed below 225 m, as the water divide is still present and high enough to drive water from the aquifer to the reservoir. For 30 and 40 ka “old aquifers”, leakage is present at all levels applied and increases linearly with the reservoir level. The actual water level of the Geheyan reservoir is 200 m, which is below the level, where significant leakage may be expected from the model.

Figure 12 Variation of hydraulic conductivity (K) in the karst aquifer before and after reservoir filling: horizontal (a) and vertical (c) distribution of hydraulic conductivity after 20 ka of natural evolution; horizontal (b) and vertical (d) distribution of K after one thousand years of reservoir conditions. The saturation ratio ($C/C_{\mathrm{eq}}$) of reservoir water is 0.8.

4.2 LKB changes before and after reservoir filling

If water in the river is aggressive to calcite, evolution of aquifer continues after the reservoir filling. We assume that the saturation ratio ($C/C_{\mathrm{eq}}$) of reservoir water is 0.8. Figure 12 shows the distribution of hydraulic conductivities in an aquifer, which has evolved naturally for 20 ka and then under the influence of reservoir for one thousand years. The subsections A8 and A13 maintained the lowest conductivity also after reservoir filling. Based on the conductivity analysis of the middle and lower primary sections, we conclude that the developed LKB remained a continuous, low-permeability barrier throughout 1000 years of reservoir filling conditions.

Figure 13 Reservoir leakage through karst aquifers over one thousand years (assuming that the aquifer lateral width is 100 m).

Figure 13 illustrates the increasing trend of reservoir leakage over a 1000-year period of reservoir operation. Leakage rates were estimated assuming a lateral aquifer width of 100 m. At the 200 m reservoir level, no leakage occurs initially, resulting in no further evolution over time. For the 250 and 300 m levels, leakage remains below 2 L s⁻¹ over the entire 1000-year span.

While the severity of leakage under these idealized conditions may appear acceptable during long-term reservoir operation, the final determination of a permissible water level must consider a broader range of factors. These include geological uncertainties, site-specific topographic features, and social safety requirements, all of which play critical roles in the decision-making process.

The fracture flow velocity distribution during the natural condition evolution process and under different water levels of the reservoir is shown in Fig. 14. Areas with low flow velocity in the central part indicate low permeability and weak karstification. The high-velocity flow region forms an arch-shaped region along the permeable fringe, which has developed along the descending water table.

In the water divide region, the LKB forms, mitigating flow beneath the fringe. Under reservoir conditions, leakage becomes substantial only when the reservoir level rises well above the LKB, causing the water divide to shift towards the reservoir and initiating flow along the permeable arch.

Figure 14 Velocity distribution in the aquifer at different stages and conditions: (a) at the onset of karstification, (b) after 20 ka of natural evolution prior to reservoir impoundment, (c) after 20 ka of natural evolution and then with a reservoir level of 200 m, and (d) after 20 ka of natural evolution and with a reservoir level of 300 m.

Figure 15 Conceptualisation of the low-karstified rock-blocks (LKB) in the water divide region (shaded area between dotted lines b and c). The zone below line c exhibits minimal karstification. Line a indicates the naturally evolving water table prior to reservoir impoundment, with conduits and caves forming on both sides of this water level.

5 Discussion and conclusions

The concept of LKB, also demonstrated by our model, is summarized in Fig. 15. The dotted line is the evolving water table under natural conditions. Below line c, the fractures have experienced almost no karstification. Between lines a and b, well-karstified zones have developed on both sides of the descending water table area. In contrast, the region between lines b and c exhibits only moderate fracture widening. This section – designated as the LKB – acts as a naturally low-permeability zone that helps mitigate potential leakage from the reservoir.

In a broader sense, karstification presents a form of nucleation, where flow-induced dissolution and changes in porosity are coupled through feedback mechanisms (Edery et al., 2021; Molins et al., 2014). In unconfined aquifers under constant recharge conditions, dissolution is most intense near the water table. This process creates a highly permeable fringe that effectively channels inflow toward both sides of the water divide. As this fringe migrates downward across the aquifer cross-section, it leaves behind a distinctive porosity imprint. Simultaneously, it inhibits deeper penetration of the inflowing solution, favouring the preferential development of horizontal fractures. Moreover, flow along the water table increases progressively from the water divide toward the discharge points. As a result the water divide zone is less karstified than the regions close to the output. Similar anisotropic, directional changes, including fingers or preferential flows, has also been observed through experimental studies and other numerical simulations (Rege and Fogler, 1989; Shavelzon and Edery, 2022; Singurindy and Berkowitz, 2003).

Reservoir leakage is highly sensitive to the relative elevation of the reservoir water level and the elevation of the low-permeability-block (LKB). Identifying the water table position within the aquifer prior to reservoir impoundment can help locate the LKB and potentially reduce the need for costly engineering measures to prevent leakage (Zhou and Li, 1996).

The presented model simplifies the system by neglecting many of the complexities observed in natural settings. Karst aquifers evolve in different climatological, hydrological and geological settings, and they often undergo multiple phases of karstification over geological timescales, driven by varying hydrological and geochemical boundary conditions. In active tectonic areas, the structure and boundary conditions may change during the karstification, which may continuously adapt to the terrain uplift and/or valley entrenchment (Gabrovšek et al., 2014). Therefore, any modelling related to a specific site must include all information on the environment in which an aquifer has evolved (Class et al., 2021; Hartmann et al., 2013). The initial structure of an aquifer is typically far more intricate than assumed in the model, potentially including pre-existing preferential flow paths, geochemical heterogeneity (e.g., interbedded layers of less soluble rock), and temporal variability in recharge. Additionally, the model does not account for the development of the vadose zone, which plays a crucial role in modifying the permeability of layers above the water table and influences the geochemical characteristics of the water at the water table. These features could be systematically incorporated into the model to investigate their influence on aquifer evolution.

Despite its limitations, the model effectively serves as a proof of concept, offering valuable insights for hydraulic engineering in karst regions.