# Improving Fault Zones Hydrodynamic Characterization and Simulation in Karstified Carbonate Environments with GLM and IES Invers Methods

Aurélie Boura, Yohann Cousquer, Victor Clauzon, Rémi Valois, Véronique Léonardi

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Aurélie Boura, Yohann Cousquer, Victor Clauzon, Rémi Valois, Véronique Léonardi. Improving Fault Zones Hydrodynamic Characterization and Simulation in Karstified Carbonate Environments with GLM and IES Invers Methods. Journal of Hydrology, 2024, 645, pp.132010. 10.1016/j.jhydrol.2024.132010. hal-04803492

HAL Id: hal-04803492

https://hal.science/hal-04803492v1

Submitted on 28 Nov 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Research papers

# Improving Fault Zones Hydrodynamic Characterization and Simulation in Karstified Carbonate Environments with GLM and IES Invers Methods

Aurélie Boura<sup>a, *,</sup> , Yohann Cousquer<sup>a</sup> , Victor Clauzon<sup>a</sup> , Rémi Valois<sup>b</sup> , Véronique Leonardi<sup>a</sup>

ARTICLE INFO

Keywords:

Fault zone

Carbonate

Parameter estimation

Karst

Parametric uncertainty

Predictive uncertainty

IES

GLMA

A B S T R A C T

A quantitative and qualitative estimation of groundwater resources is a challenge for water management worldwide. In the specific case of the carbonate karstified aquifer, the high sensitivity of this resource is partially due to its level of permeability. Moreover, the extreme heterogeneity of this type of aquifer makes it complex to define vulnerable areas, especially when fault zones exert a significant influence on groundwater hydrodynamics and transport. Improving the management of this type of aquifer involves better assessment of its hydrodynamics properties. Sustainability and management solutions can be estimated through numerical modeling. In this study hydrodynamic parameter of the fault zone in karstified carbonate environments are estimated from cross-hole pumping test data. These analyses allow the simulation and forecasting of fault zone flow and travel time. For this purpose, a flow and travel time model has been set up using MODFLOW6 and MODPATH7 codes. Inverse modeling of hydrodynamic parameters is performed with the PEST++ code suite, with both the Gauss-Levenberg-Marquardt Algorithm (GLM) and the Iterative Ensemble Smoother (IES) which is a recent approach with few applications so far. Simulation, parameter estimation, and forecasts of drawdown and pollutant travel time through the fault zone have been first evaluated on a synthetic fault zone and are eventually applied to a real-world fault zone of interest in karst carbonate environments. A comprehensive comparison and analysis of GLM and IES results are performed for parameter estimation, parametric and predictive uncertainties, ensuring a thorough assessment of results. The numerical inversion approach provided a spatial parametric estimation of the hydraulic conductivity of the fault zone, along with an uncertainty analysis.

This approach demonstrates its ability to estimate the hydrodynamic parameter field with a high level of accuracy within the internal structures of the fault zone. This study constitutes valuable and reproducible in- sights into the simulation and forecasts of fault zones in karst carbonate aquifers. The findings are a step forward in groundwater flow and travel time simulation of fault zones and provide relevant practical help for the management and protection of water resources in these critical areas of interest.

# 1. Introduction

Fault zones (FZ) play a significant role in groundwater circulation, which can have significant consequences on water resources availabil- ity, heat flux, and solute transport (Bense et al., 2013; Jones and Detwiler, 2016; Taillefer et al., 2021). They can either facilitate fast flow leading to preferential infiltration of pollutants and rapid transfers or act as barriers to water flow (Cherry et al., 2009; Gassiat et al., 2013). Fault zones can promote exchanges between aquifer units like deep and shallow aquifers, posing a risk to the quality of affected drinking water (Maslia and Prowell, 1990; Iskandar and Koike, 2011; Williams et al.,

2013). Fault zones can also act as barriers to groundwater circulation due to the physical damage of the rock, leading to a reduction in its permeability (Evans et al., 1997; Rawling et al., 2001; Micarelli et al., 2006; Anderson and Fairley, 2008). These differences in groundwater circulation behaviors, whether they reduce or facilitate flows, can be attributed to fault structure.

Fault zones are generally structurally divided into three entities: the core zone (CZ), the damage zone (DZ), and the healthy undamaged protolith (Shipton and Cowie, 2001; Billi et al., 2003; Mitchell and Faulkner, 2009; Savage and Brodsky, 2011). Each of these structures have distinct hydrodynamic properties (Caine et al., 1996; Evans et al.,

<footnote><sup>a</sup> HSM, Univ. Montpellier, CNRS, IMT, IRD, Montpellier, France</footnote>

<footnote><sup>b</sup> Avignon Université, EMMAH UMR1114, F-84000, Avignon, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: aurelie.boura1@umontpellier.fr (A. Boura).</footnote>
1997; Faulkner et al., 2010), which depend on the protolith rock, the type of the fault displacement (normal, reverse, strike-slip), and the type of the rock composing the CZ (gouge, breccia, cataclasite) (Sibson, 1977; Caine et al., 1996; Agosta et al., 2007; Bauer et al., 2015).

Fault zones in cohesive media, such as crystalline, clastic and carbonate rocks, typically exhibit a well-developed network of interconnected fractures within the DZ (Agosta et al., 2007). This fracture network results in increased permeability compared to the protolith effectively acting as conduits in groundwater flow (Faulkner, 2004; Micarelli et al., 2006; Balsamo et al., 2010; Faulkner et al., 2010). The CZ is characterized by low permeability and reduced porosity, acting as a hydraulic barrier. This reduced permeability is attributed to the mechanical and/or chemical deformation of the rocks (Agosta and Kirschner, 2003; Micarelli et al., 2006; Agosta et al., 2007). As a result, the fault zone can exhibit a dual behavior of both conduit and barrier, with the DZ acting as the conduit zone and the CZ acting as the barrier zone. In carbonate rocks, karst dissolution can occur and overlay on this specific fault zone structure, inducing either conduit or barrier behavior (Drogue, 1972; Ford and Ewers, 1978; Andre and Rajaram, 2005; Ford and Williams, 2007; Sauro, 2014; Littva et al., 2015; Li et al., 2019; Clauzon et al., 2020). This trend introduces a great deal of uncertainty regarding underground circulation. Karst conduits preferentially develop along fractures, significantly increasing the hydrodynamic properties of the reservoir by several orders of magnitude (PALOC, 1964; Larsen et al., 2010; Nooruddin et al., 2014; Robineau et al., 2018; Clauzon, 2019). Consequently, in fault zones, the karst network plays a crucial role in groundwater circulation at the aquifer scale. Considering the importance of the fault zones on groundwater hydrodynamics and for the protection and management of the groundwater resource, it is crucial to determine both the behavior of faults acting as conduit and/or barrier through the study of their hydrodynamic properties.

In this perspective, spatial estimation of fault zone hydrodynamic properties is generally based on geological, hydrogeological, or numerical approaches (Bense et al., 2013). Currently, the estimation of hydrodynamic properties of the fault zone is generally conducted through structural and geological analysis only, complemented sometimes by measuring permeability along the fault zone using a permeameter (Antonellini and Aydin, 1994; Rawling et al., 2001; Balsamo and Storti, 2010; Mitchell and Faulkner, 2012; Matsumoto and Shigematsu, 2018; Clarkson et al., 2019). Those methods are dependent on the representativity and scale of the measured samples, and they are often only representative of surface observations and not groundwater circulation (Evans et al., 1997; Agosta et al., 2007). The various hydrogeological methodologies used to characterize the hydrodynamics of the fault zone are primarily based on data related to the hydraulic properties of the reservoir (Singhal and Gupta, 2010). Consequently, these approaches tend to represent groundwater dynamics at the global reservoir scale (Bense et al., 2013) and therefore cannot fully reflect the actual impact of the fault zone in detail, due to the scale inherent to these methods. In highly heterogeneous environments such as karstified carbonate settings, these methods are not capable to capture the structural complexity and therefore the spatial distribution of their hydrodynamic properties.

The management of vulnerability zones in fractured aquifers, especially in karstic environments, is of crucial importance. Numerical applications provide prediction methods in these fractured environments. However, few consider the effects of conduits and barriers in fractures and fault zones (Day-Lewis et al., 2000; Yang et al., 2016). Models that simulate flows through these fault zones rely on highly detailed geological and structural descriptions, making their methodology difficult to reproduce (Wilson et al., 2003; Fredman et al., 2007; Leray et al., 2013; Romano et al., 2017). Without using hydrological observations to validate their results or propose a distribution of hydrodynamic properties. This underscores the significant absence of study for both hydrodynamically interpreting fault zones and accurately simulating and predicting flows, particularly in managing the vulnerability of fractured

aquifers.

Most hydrodynamic studies in fault zones are conducted in granitic, siliciclastic rocks or clastic sediments (Rawling et al., 2001; Jourde et al., 2002; Wibberley and Shimamoto, 2003; Childs et al., 2007; Rotevatn et al., 2008; Faoro et al., 2009; Balsamo and Storti, 2010; Mitchell and Faulkner, 2012; Bigi et al., 2013; Fisher et al., 2018; Lapperre et al., 2022) and a few in carbonate environments (Rotevatn and Bastesen, 2014; Korneva et al., 2014; Haines et al., 2016). Therefore, there is currently no description of the hydrodynamic properties and their magnitude within a fault zone in karstified environments, nor any work on the spatialized study of the hydrodynamic properties of fault zones in these environments. With around 20 % of the world's water resources coming from karst aquifers (Ford and Williams, 2007), there is a need to improve the understanding of hydrodynamic properties estimation of the fault zone for better and more appropriate resource management decision-making.

Hydrogeological modeling is a powerful decision-making aid for water resources management in general (Aster et al., 2013) and for karst aquifers (Scanlon et al., 2003; van Beynen et al., 2012; Jeannin et al., 2021; Sivelle et al., 2021; Alarcon et al., 2022; Cousquer and Jourde, 2022). Most hydrogeological models deal with fault zones by using expert knowledge based on field observations at the outcrop (Tveranger et al., 2005; Romano et al., 2017; de Souza et al., 2022), and hydrodynamics properties of the fault zone are rarely constrained by hydrodynamics observations.

The inverse methods are widely used in the estimation of the groundwater parameters (Aster et al., 2018) as in carbonated and karstified environments (Dörfliger et al., 2009; Lu et al., 2012; Wang et al., 2016; Giese et al., 2018; Kavousi et al., 2020; Fischer et al., 2017), as well as for estimated karst systems structure (Jaquet et al., 2004; Illman, 2014; Zha et al., 2015; Wang et al., 2016; Collon et al., 2017; Fischer et al., 2018).

However, inverse problems in hydrogeology are often, if not always, ill-posed, i.e. suffer from the non-uniqueness, non-steadiness, or non-existence of the solutions (Hunt et al., 2007; Doherty, 2011; Zhou et al., 2014). Algorithmic methods for history-matching associated with numerical solutions make it possible to tend towards a well-posed problem with for example mathematical regularization with singular value decomposition (SVD) or Tikhonov regularization (Carrera et al., 2005; Zhou et al., 2014; Doherty, 2015; White, 2018). These solutions are implemented in the codes PEST and PEST++ (Doherty, 2015; White et al., 2020a) and are widely employed for parameter estimation and groundwater model analysis (Zhuang et al., 2017; Cousquer et al., 2018; Essouayed et al., 2021, Pryet et al., 2022; Delottier et al., 2022). These codes aim to solve the inverse problem based on the deterministic method Gauss-Levenberg-Marquardt Algorithm (GLM) (Doherty, 2015) and more recently on the ensemble method the Iterative Ensemble Smoother (IES) (Chen and Oliver, 2013). The GLM method is widely used for parameter estimation in hydrogeological modeling (Commer et al., 2014; Burrows and Doherty, 2016; Aster et al., 2018) and the IES method is much more recent with fewer applications in hydrogeology (Coulon et al., 2022; Delottier et al., 2022) and any application at a local scale (i.e. 100 m – 1 km). Both methods are implemented in the PEST++ suite (White et al., 2020a). PEST++ GLM uses a generalized linear regression approach to estimate the relationship between observed and simulated values, employing the Jacobian matrix to solve optimization problems. On the other hand, PEST++ IES is based on the same linear regression approach, but it uses an ensemble-based method to estimate the posterior distribution of model parameters. PEST++ GLM and IES are designed for parameter estimation and for parametric and predictive uncertainty quantification, with the latter directly obtained from PEST++ IES ensemble and implemented in PEST++ GLM using the First-Order, Second-Moment (PEST++ GLM-FOSM package) method (Doherty, 2015).

The objective of this work is to propose a spatial distribution of hydrodynamic parameters within a fault zone in a carbonate and karstified
aquifer. This will be achieved using data derived from hydrogeological approaches, particularly the evolution of groundwater heads during pumping tests and the corresponding discharge rates.

In this aim, we propose to utilize groundwater level variations during aquifer tests to obtain a 2D spatial estimation of the hydraulic conductivity field of the fault zone. This approach facilitates our understanding of groundwater flow dynamics. We aim to determine if these fluctuations can effectively characterize hydrodynamic behavior in fault zones through the matrix and conduit system typical of karst environments. Subsequently, based on this preliminary estimation, our approach enables the simulation and prediction of flow through the fault zone, considering parametric and predictive uncertainties obtained from GLM and IES methods. The direct integration of GLM and IES methods into PEST++ provides access to a comparison of both approaches within the same methodological configuration. To assess the relevance of our approach, we conducted a comparison using a synthetic case study where the “truth” hydrodynamic field is known. The synthetic fault zone case was constructed based on existing knowledge, relying on known thickness relationships relative to fault scale and spatially described hydraulic conductivity data variations. This comparison was performed with a limited number of wells to evaluate the method’s performance, considering the constraints imposed by limited spatial information. This initial step allows us to test the method’s ability to capture the heterogeneity of the hydrodynamic properties of the fault zone, particularly concerning the magnitudes of hydraulic conductivity encountered in carbonated environments. The same validated approach is then tested on a real case study, supported by previously acquired geological information and geophysical data. This study aims to evaluate the methodology developed in a carbonate karst environment, determining its effectiveness despite the presence of a real karst network. It also provides an understanding of the interaction between the fault zone and the karst environment. The benefits and drawbacks of each of the proposed methods are subsequently examined.

## 2. Workflow

The numeric methodology used to address this issue is described in the following workflow. The developed workflow follows the steps presented in Fig. 1. This workflow is based solely on open-access and flexible codes, communicating with Python wrappers, and is designed to

be transposable to other studies. A transient MODFLOW6 flow model that simulates multi cross-hole pumping test is built using the Python wrapper FloPy (Bakker et al., 2016; Leaf and Fienen, 2022) where geographical entities are pre-processed through the geographic information system QGIS interface (Fig. 1 step 1). Model meshes are generated from the Python package GridGen (Fig. 1 steps 2) (Lien et al., 2015). Then, Flopy enables pre- and post-processes MODFLOW6 input and output files (Fig. 1 steps 2 and 3). The pre-and post-processing for parameter estimation code PEST++ input and output files is performed using the PyEMU Python package (Fig. 1 steps 2 and 3) (White et al., 2016) (Fig. 1 step 4). Parameter estimations are then performed with both GLM and IES methods using PEST++ suite (White et al., 2020). The steady-state MODFLOW6 and MODPATH7 models with estimated parameters for flow and particle tracking for advective transport simulation forecasts are pre and post process using FloPy Python package (Fig. 1 step 5). Both forecasts of steady-state flow and particle tracking with predictive uncertainty are eventually performed using both the PEST++ GLM-FOSM and the ensemble of PEST++IES models (Fig. 1 step 6). This workflow is identical to the synthetic and the real case study. Each step is detailed in the following sections.

### 2.1. Flow and transport model

In flow modeling, two models will be distinguished: i) the transient flow model which is employed for parameter estimation purpose, and ii) the steady-state flow and advective transport model for prediction based on the previous estimated parameter field (Fig. 1). The transient groundwater flow model is based on an equivalent porous media approach accounting for the hydrodynamic behavior of both fault and karst features at a medium scale (100 m – 1 km) (Abusaada and Sauter, 2013; Jeannin et al., 2021). Flows are simulated using the MODFLOW 6 computational code (Langevin et al., 2022), this model is used to simulate cross-hole pumping tests for the parameter estimation purpose. An identical but steady-state model is also used to simulate flow and advective transport for forecast purposes using MODFLOW6 and MODPATH7 codes (Pollock, 2016). The advective transport, with the particle tracking method, is used here as a surrogate model for advective-dispersive transport. Dispersion can be neglected here since it simulates transport between a point source and a pumping well (Cousquer et al., 2018; Pryet et al., 2022). This allows forecast

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Workflow description: 1) spatial information contained in QGIS 2) pre-processing of geometric information of step 1 and inputs for MODFLOW 3) cross-hole pumping tests simulations 4) Solving the inverse problem with both GLM and IES from PEST++ code 5) flow and travel time simulation with estimated parameters and parametric 6) predictive uncertainty quantification with both FOSM and PEST++ IES.</figcaption></figure>
uncertainty quantification.

The focus is on characterizing the influence of the fault zone on groundwater flow, rather than on the specific travel time estimation within the fractured environment such as karst. We employ an equiva- lent porous medium model, a decision motivated by the fact that fault rocks in general and in karst environment are predominantly matrix entities (Balsamo et al., 2010) and because the objective is for forecast purpose taking into account both fault zone entities and karstic objects in an equivalent porous media as it has already been made for karstic environment (Wang et al., 2016; Abusaada and Sauter 2013).

## 2.2. PEST++ GLM and PEST++ IES for parameter estimation

In this section, a brief overview is provided on how parametric estimation functions using the PEST++ software, along with the methods GLM and IES implemented in the developed workflow.

In the stage of parametric estimation, the hydraulic conductivity results from the inversion of the transient flow model. This parametric estimation is spatially adjusted during the inversion. It is parameterized using the widely employed pilot points method (de Marsily et al., 1984; Doherty, 2003; White and Lavenue, 2023). The method is extended to the calculation grid by kriging. An exponential variogram is defined with a threshold of one, a nugget of 0.1, and a range that is twice the distance between pilot points (Kitanidis, 1997). Storage is also estimated and considered homogeneous.

Two methods of parameter estimation by historical matching are examined i) one deterministic, based on the gradient method with the widely used Gauss-Levenberg-Marquardt (GLM) algorithm, and (ii) the ensemble method, specifically the more recent Levenberg Marquart form of the Iterative Ensemble Smoother (IES) (Chen and Oliver, 2013; White, 2018).

Estimation using PEST++ GLM and IES is based on the GLM method, which works by reducing an objective function using the gradient descent method. The objective function $\Phi$, used by the GLM method, corresponds to the weighted sum of squares of the residuals between the simulated and observed values:

$$ \Phi = \sum_{i=1}^{n} w_i (x_{i,\text{obs}} - x_{i,\text{sim}})^2 \quad (1) $$

where $\Phi$ is the objective function to be minimized, $n$ is the total number of observation, $w_i$ is the weight associated to the objective function, $x_{i,\text{obs}}$ is the observed value and $x_{i,\text{sim}}$ is the simulated value.

The differences between GLM and IES are, firstly, that for each iteration PEST++GLM linearly solves the inverse problem by calculating the Jacobin matrix for each iteration, whereas PEST++ IES uses the probabilistic technique of Monte Carlo analysis to obtain an approximation of the Jacobin matrix. The approximation of the Jacobian matrix used is IES algorithm is based on the empirical cross-covariance between the parameters and the observations, and systematically provides the deviations around the mean value of the parameters (Chen and Oliver, 2013; White, 2018).

For each iteration, IES probabilistically generates a set of parameters, which is then linearly calibrated using the same principal as GLM method. This approach thus yields a set of parameters that allows solving the inverse problem. Thus, the IES results are interpreted as the average of the parameter estimates at each pilot point, whereas the PEST++ GLM estimation results attest to a single parameter set.

In our case, where drawdown data is used to constrain the flow model, a large number of parameters, through pilot points, are needed to cover the spatial extent of the fault zone. In our application, the major advantage of IES is that it avoids the high costs associated with estimating a large number of parameters due to its low computational requirements, as demonstrated in large-scale cases (Hunt et al., 2021). This method's computational efficiency which remains nearly unaffected by parameter count, makes it suitable for high-dimensional problems like ours (White, 2018).

Mathematical regularization and preferential parameter Tikhonov regularization methods are used to overcome the non-uniqueness solution problems in such hydrogeological models (Doherty and Hunt, 2010; Doherty, 2015): Tikhonov regularization (Tihonov, 1963; Doherty and Hunt, 2010) and singular value decomposition (SVD) (Tonkin and Doherty, 2005).

Parametric and predictive uncertainty on flow and transport simulation is eventually quantified with both the first-order, second-moment method (FOSM) embedded in the GLM algorithm and with IES.

## 2.3. Parametric and predictive uncertainty

Parametric and predictive uncertainty analysis with GLM is based on FOSM (White et al., 2020), and conducted with pyEMU (White et al., 2016). This method is used for propagating uncertainty by considering parameter distribution and measurement error. The uncertainty analysis on the posterior parameter set and then model predictions are propagated from prior parameters and measurement error by model linearization with the application of Bayes' theorem. This kind of linear-based analysis is known to be approximative but less time-consuming than nonlinear methods (Hill et al., 2016; Anderson, 2015).

Uncertainty is also estimated with a non-linear method from IES inversion results. In this way, IES algorithm provides an improved analysis of uncertainty thanks to the analysis between the prior and the posterior parametric and predictive distribution of the model ensemble.

# 3. Synthetic case

The described method embedded in the workflow (Fig. 1) is first evaluated on a synthetic case where a fault zone hydraulic conductivity field in a carbonate environment has been generated. This hydraulic conductivity field is generated to represent a conduit-barrier fault zone characterized by an impermeable fault rock constituting the CZ of the fault, and a highly permeable DZ.

## 3.1. Synthetic hydraulic conductivity field generation

The hydraulic conductivity field is defined as a North-South trending fault zone with a 10 m wide CZ which corresponds to a fault displacement of 10 to 100 m (Wibberley et al., 2008; Childs et al., 2009). For these orders of magnitude, the associated DZ thicknesses vary between 50 and 100 m (Balsamo et al., 2019). The chosen width of the DZ is 80 m.

The CZ is known to have a reduction of hydraulic conductivity of about three orders of magnitude compared to the healthy protolith (Balsamo and Storti, 2010; Rath et al., 2011). In the DZ, the hydraulic conductivity can increase by up to two orders of magnitude compared to the protolith (Antonellini and Aydin, 1994; Balsamo et al., 2019). Within the DZ, hydraulic conductivity decreases exponentially away from the CZ (Bastesen and Braathen, 2010; Faulkner et al., 2010; Mitchell and Faulkner, 2012). In the context of a carbonate aquifer, studies such as those by Micarelli et al. (2006), Hausegger et al. (2010), Bauer et al. (2016), and Geng et al. (2021) have identified similar characteristics from fault zones to barrier conduits. However, fracturing can lead to the development of karst networks, where hydraulic conductivity varies based on factors such as connectivity, openness, and frequency (Korneva et al., 2014; Haines et al., 2016; Ferraro et al., 2020).

Different studies have shown significant variations in hydraulic conductivity between extreme values because of karstification (Wibberley et al., 2008; Leonardi et al., 2016; Clauzon et al., 2020; and Wu et al., 2020).

To represent these variations, we generated a hydraulic conductivity field using an exponential model based on the well-characterized patterns of hydrodynamic property changes perpendicular to the fault zone (Balsamo et al., 2010; Bastesen and Braathen, 2010; Faulkner et al., 2010; Mitchell and Faulkner, 2012). This model incorporates hydraulic
conductivities on the order of magnitude typical in karst environments. The highest and lowest permeability values were used to describe the DZ and the CZ, respectively (Fig. 2B). The hydraulic conductivity field along the axis perpendicular to the fault zone is characterized by an hydraulic conductivity of the protolith ($K_{\text{protolith}}$) of $10^{-5} \text{ m s}^{-1}$, the CZ by a $K_{\text{CZ}}$ of $10^{-7} \text{ m s}^{-1}$ and the DZ with a $K_{\text{DZ}}$ varying exponentially from $10^{-3} \text{ m s}^{-1}$ to $K_{\text{protolith}}$ over the width of the DZ (Fig. 3) (Clauzon, 2019).

## 3.2. Flow and transport model

The domain is discretized on a 300 x 300 m grid with a regular 1.25 m rectangular mesh. Boundary conditions are defined by head-dependent flux (Cauchy-type) with the General Head Boundary (GHB) MODFLOW6 packages on the North and South boundaries with a head of 50 m and a conductance of $10^{-3} \text{ m}^2 \text{ s}^{-1}$. East and West boundary conditions are no flow. Transient flow conditions are considered, and the transmissivity is assumed to be independent of water level fluctuations,

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Schematic representations of fault zone structures and associated variations in hydraulic conductivity for both karstified and non-karstified media. A: fault zone features in a non-karstified carbonate. B: average evolution of hydraulic conductivity perpendicular to the central segment of the fault (A). (After Faulkner et al., 2010).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Generated hydraulic conductivity of the synthetic case. A: Construct variations of the hydraulic conductivity perpendicular to the central segment of the synthetic fault. B: Map of the hydraulic conductivity of the synthetic fault zone and wells obtained for the direct simulation on MODFLOW 6 of the synthetic pumping test through the Fig. 3.A distribution of the hydraulic conductivity.</figcaption></figure>

which corresponds to the Boussinesq assumption. The five wells are considered fully penetrating the entire thickness of the domain, which is 30 m. The cross-hole pumping tests are simulated at a 10-minute time step for a total duration of 60 min, with flow rates imposed in wells W2, W4, and W5 set to $Q = 1.33 \times 10^{-3} \text{ m}^3/\text{s}$ and in wells W1 and W3 set to $Q = 1.7 \times 10^{-7} \text{ m}^3/\text{s}$ (Fig. 3). The imposed flow rates for the pumping tests were determined to avoid dewatering the wells and thus maintain realistic aquifer test conditions. Therefore, in the case of wells located in the impermeable core zone with a hydraulic conductivity of $10^{-7} \text{ m s}^{-1}$, the flow rates are very low. Direct simulation of pumping tests provides the observations used in the inverse model. Specifically, we simulate the pumping tests using MODFLOW, utilizing the synthetic hydraulic conductivity field and various discharge rates to generate our observations of hydraulic heads. This method ensures that the observations accurately reflect the synthetic case pumping scenarios, enabling us to construct a dataset for inverse modeling. The steady-state transport simulation with particle tracking is performed with a forward particle placed in W2 well. A flow rate of $Q = 1.4 \times 10^{-3} \text{ m}^3/\text{s}$ is extracted from well 1 and travel time between W2 and W4, which are located on both sides of the fault, is calculated.

## 3.3. Parameter estimation, parametric and predictive uncertainties

The hydraulic conductivity field was parameterized using a regular network of 195 adjustable pilot points spaced at 10-meter intervals, distributed directly across the mesh throughout the domain. The parameter estimations from the pilot points are interpolated to the remaining cells at the grid cell centers, depending on their locations. The 195 pilot points are free within the zone of influence of the pumping
tests, so they are considered independent. Beyond this zone of influence, the rest of the pilot points are tied in areas where information is lacking, thus creating a buffer zone between the boundary conditions and the zone of interest.

For parametric estimation using the GLM method, the initial parameters are defined based on exploring the objective function within the parameter space, yielding an initial value for hydraulic conductivity estimation of $10^{-5}$ m s$^{-1}$. For both parametric estimation methods, the prior parameter distribution is set between $10^{-3}$ and $10^{3}$, enabling the use of a high-dimensional parameter space to achieve robust parameter estimates (Doherty, 2011).

The choice of the number of iterations for the two inversion methods initially involved testing to find the right balance between sufficient reduction of the objective function and moderating the overfitting (White et al., 2020a). For parameter estimations based on the GLM method, we observed the reduction of the objective function alongside the evolution of the regularization objective function to select the iteration that provided a good compromise between the two (Fig. 4A). For IES, the ensemble approach allows for rapid reduction of the objective function in a small number of iterations (White et al., 2020b; Abbas et al., 2024). Additionally, increasing the number of iterations for this approach tends to make the ensemble of models converge into a single unique model, nullifying the method's benefit. Thus, in our work, the number of iterations imposed on parametric estimation with IES was defined as 3 iterations, which is most of the time the chosen iteration (White et al., 2020a).

In the interest of balancing the representativeness of the observations

in the objective function for the estimation of the hydraulic conductivity (Doherty and Welter, 2010), we defined weights for the set of observations by:

$$ \omega_{ij} = \frac{1}{\sigma_{ij}} \qquad (2) $$

where $\omega_{ij}$ the weight of the observation variable corresponding to the well $i$ and the period $j$, $\sigma_{ij}$ represents the standard deviation of the observations is calculated as the square root of a value derived from comparing observation points to the mean of head observations from each pumping test conducted at every well. Every observation of each pumping test within each well is appropriately weighted in the pre-calibration objective function.

## 3.4. Results

### 3.4.1. Fault zone parameter estimation with GLM and IES

The parameter estimation was parallelized on 90-core CPU for both GLM estimation and IES, with the estimation based on 2755 head observations. Quantitative results are summarized in Table 1. Fig. 4 shows the rapid reduction of the objective function in the first 10 iterations of the inverse model via PEST++ GLM and for the third with PEST++ IES according to White et al. (2020a) recommendation to avoid overestimation that leads to a unique model.

The estimation results demonstrate a successful convergence of the models, as evidenced by the correlations obtained between simulated

<figure><img src="image_7.png" /><figcaption>Fig. 4. Parameter estimation results for the synthetic case. A: evolution of the measurement objective function and the regularized objective function for GLM iteration. B: correlation between observed and simulated head from GLM estimation. C: evolution of the objective functions of the 300 IES models. D: correlation between observed and simulated head with prior and posterior IES models.</figcaption></figure>
<table><caption>**Table 1**
Summary of optimization results.</caption><thead><tr><th>Optimization results</th><th>PEST++GLM</th><th>PEST++IES</th></tr></thead><tbody><tr><td>Number of PEST++ iteration</td><td>30</td><td>3</td></tr><tr><td>Number of model execution</td><td>9554</td><td>900</td></tr><tr><td>Computing time on 90-core CPU (hours)</td><td>26.69</td><td>1.28</td></tr><tr><td>RMSE (m)</td><td>0.29</td><td>Range: 0.61 – 2.33<br>Ensemble: 1.48</td></tr></tbody></table>

and observed hydraulic heads in both the IES and GLM inverse models (Fig. 4B, 4D).

Fig. 4A shows the measurement objective function reduction according to the increase of the regularization objective function. This result determines the best trade-off that strikes a balance between measurement and regularization objective function. In this case, the 12th iteration is a good trade-off. The correlations achieved between the estimated and observed head values by IES indicate an average root mean square error (RMSE) between 0.61 and 2.33 m with an average value of 1.48 m. In comparison, GLM demonstrates an RMSE of 0.29 m.

Generally, the head estimates obtained from GLM and IES bracket the behavior observed in the pumping test simulations, as illustrated in Fig. 5 with pumping conducted in well W1 and observation in well W2. The head simulations obtained from the parameter estimation based on GLM show better calibration of the heads in the pumped well W1 than in W2 (Fig. 5). The posterior ensemble obtained from IES exhibits a wide range around the synthetic observed heads, but brackets them relatively well (Fig. 5).

The estimated hydraulic conductivity field (Fig. 6A) shows an alignment of low hydraulic conductivity values corresponding to the synthetic fault pattern (Fig. 3B, the synthetic alignment). The observation of the residual field is in Fig. 6B, i.e. the differences between the real synthetic and estimated conductivity fields show that the residuals are low in the area where the boreholes are close to each other on both sides

of the fault zone. The observation of the residual map (Fig. 6B) shows that the model is less efficient in estimating the correct location of the fault barrier zone in the southern part of the test area where a well is located only on one side of the fault zone.

A comparison of the hydraulic conductivity fields estimated by IES and GLM reveals certain differences, particularly in the lower part of the fault zone (Fig. 6A, 6C). In this zone, the average of the hydraulic conductivity fields obtained by IES shows a poorer definition of the synthetic fault segment, which, as in the residual field map (Fig. 6D), indicates the lack of spatialized data in this zone. The maps representing the residuals between observed and simulated values (Fig. 6B, 6D) indicate that the IES generally provides a better estimation of the conductivity field, especially for the DZ and the protolith. IES is significantly more successful in distinguishing the transition phase of the DZ with high hydraulic conductivity values. In contrast, GLM does not clearly show the DZ but only highlights the transition from the impermeable CZ.

### 3.4.2. Simulation and forecasts uncertainty analysis with FOSM and IES

Drawdown and travel time forecasts under uncertainty have been estimated for well 4 with pumping in well 2 with parameters previously estimated (Fig. 6). Predictive uncertainty from GLM estimated parameters is performed with FOSM and with IES for IES results. The posterior head and travel time obtained with FOSM and PEST++ IES is greatly improved over the prior and is close to the “true” travel time.

The posterior mean travel time estimated by FOSM is 1,214 days, with a standard deviation of 446 days, compared to a “true” value of 987 days (Fig. 7A, Fig. 7B). Using IES, the posterior mean travel time is 961 days with a standard deviation of 447 days, which is very close to the true value and similar to the standard deviation obtained with FOSM (Fig. 6C, Fig. 6D).

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Observed vs. simulated plots of synthetic case hydrographs for pumping in well W1 and observation in well W2. The actual synthetic case heads are in red, the head estimates obtained by GLM are in black, and all the heads obtained with IES are in gray. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
## PEST++ GLM estimation

<figure><img src="image_9.png" /></figure>

## PEST++ IES estimation

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Fault zone parameter estimation with GLM and IES. A and C: Estimated Hydraulic conductivity field for the synthetic fault zone with GLM (A) and IES (C). B: Residual field of the GLM (B) and IES (D) estimated hydraulic conductivity.</figcaption></figure>

# 4. Real world case study

The method described above has been applied on a case-study site located in the southeast of France in a karst aquifer, which supplies drinking water to the Montpellier urban area, i.e. 350 000 ca (Fig. 8). The experimental site is situated in a NW-verging monocline with stratigraphic dips ranging from 25° to 30°. It is characterized by surface geology composed of Berriasian limestone and marl-limestone alternations from the Valanginian. The section of the aquifer studied at this site is located along the eastern border of the catchment area, near the outfall (2.3 km). The aquifer layer represented corresponds to the Berriasian limestone.

An N140°E-oriented fault (Rieu Coullon) crosses the site, extending over several hundred meters. The observed 50 m dextral offset juxtaposes the Berriasian limestone with the Valanginian marl-limestone. A heterogeneous and altered zone, with a width ranging from 50 m to 60 m, is observed, potentially associated with the damage zone. The lack of outcrops and clear zones along this fault makes it challenging to directly observe its composition or find structural indicators. As a result, geological field approaches do not provide an initial insight into the hydrodynamics of the fault zone. Five boreholes (F1 to F5), penetrating vertically into the Valanginian and upper Berriasian for approximately 60 m, were strategically placed to study the fault and are positioned on either side of it (Fig. 8).

Clauzon (2019) observed the influence of the Rieu Coullon fault

fracturing through the realization of pumping tests and measurements with flowmeters at all boreholes on the site. This was evidenced by the NNW-SSE direction of the observed hydraulic gradient. The karst network is notably sparse around the fault, as revealed by optical borehole imaging.

## 4.1. Geophysical information

Geophysical investigations have been conducted on the site. The resistivity panel crosses the Rieu Coullon fault perpendicularly, over a profile length of 240 m, with 120 electrodes spaced 2 m apart. The electrical resistivity tomography (ERT) profile is made from the Syscal Pro datalogger (IRIS instruments) using the Wenner-Schlumberger array acquisition configuration. Measurements with a standard deviation higher than 5 % were removed. Data was inverted using Res2dinv software (Loke, 1999).

The electromagnetic (EM) resistivity map was made from the acquisition of the average variations of the electrical conductivity of the soil, at depths between 0 and 5.5 m. Measurements were carried out with an EM31 (Geonics) in a vertical dipole configuration connected to a GPS. Measurements were recorded each second along lines oriented perpendicularly to the faults and with lines spaced approximately 5 to 10 m from each other.

Kriging using a fitted pentashperical model in Surfer (Golden Software) was used on a 1 m x 1 m 2D grid.
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Travel time prediction results and head estimation with FOSM (A and B) and IES (C and D).</figcaption></figure>

Geometric information from the ERT profile shows a sub-vertical structure in the center of the profile at the level of the central segment of the Rieu Coullon Fault (Fig. 9). Variations in apparent conductivity on the electrical conductivity map (Fig. 9) highlight the dextral strike-slip

of the Rieu Coullon fault, as interpreted by the surface geology (Clauzon, 2019).

The geological knowledge and the geophysical information acquired on the experimental site of Rieu Coullon (Fig. 8, Fig. 9), is used as a qualitative criterion for the posterior definition of the fault's shape for the estimation of the hydraulic conductivity field.

## 4.2. Model description

The Rieu Coullon experimental site is spatially discretized using a rectangle of 1200x800 meters. We are imposing an upper layer thickness of 30 m based on observations of variations in the site's water table. The transient flow for parameter estimation and the steady-state flow and transport model for forecasting purpose have the same boundary condition, temporal and spatial discretization. Parameter estimation methodology of this case study was carried out using the same workflow as the synthetic case one.

The head data set represents four well tests carried out in the different wells with observations in all wells for a total of 20 observation series (Table 2). Pressure data is recorded from TD-Diver probes at frequencies of 10 to 30 s. The drawdown data is corrected for the natural recession of the water table.

The pumping tests selected for our study were conducted during comparable hydrological periods, reflecting similar uninfluenced groundwater levels. In the presence of variations in the pumping regimes, each well reached a permanent state at different times. The imposed flow rates were adjusted to avoid well dewatering.

The bounds for the prior values in parameter estimation were defined based on the hydrodynamic knowledge of the site, with the minimum and maximum values of hydraulic conductivity ranging from $10^{-3}$ to $10^{-10} \text{ m s}^{-1}$ (Clauzon, 2019). Exploration of the objective function in parameter space for a homogeneous hydraulic conductivity and storativity identified initial prior values of $10^{-5} \text{ m s}^{-1}$ and $10^{-5}$ respectively.

## 4.3. Case study parameter estimation

The parametric estimations were based on a set of 7140 head

<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Rieu Coullon experimental site. Aerial map with geological interpretation (modified from Clauzon, 2019), and location of the five boreholes (F1 to F5).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Results of geophysical applications at the Rieu Coullon site. A: The electrical conductivity map of the ground, in red, is the axis of the fault interpreted on the geological map, and in the dotted line is the ERT panel line (realized with EM31). B: ERT panels intersect the Rieu Coullon fault perpendicularly with the interpretation of the fault zone (FZ) location. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>**Table 2**
Setting conditions of cross-hole pumping test realized in four boreholes (F2, F3, F4, F5) (location present in Fig. 8).</caption><thead><tr><td>Test Date</td><td>Pumped well</td><td>Flow rates (m³/s)</td><td>Test Duration</td></tr></thead><tbody><tr><td>12-July-2018</td><td>F5</td><td>7.30 10<sup>-3</sup></td><td>5 h 00</td></tr><tr><td>07-Sep-2018</td><td>F4</td><td>5.33 10<sup>-3</sup></td><td>51 h 14</td></tr><tr><td>18-Sep-2028</td><td>F2</td><td>4.25 10<sup>-3</sup></td><td>58 h 50</td></tr><tr><td>22-Oct-2018</td><td>F3</td><td>1.46 10<sup>-2</sup></td><td>120 h 57</td></tr></tbody></table>

observations. The correlations between the observed values and the simulated values are satisfactory in view of the results of the average RMSE from the PEST++GLM and PEST++IES estimates (Table 3).

Fig. 10A, and Fig. 10C show the hydraulic conductivity fields estimated with PEST++GLM and PEST++IES, both demonstrate the heterogeneity of the hydraulic conductivity field resulting from observations in cross-hole pumping tests. The hydraulic conductivity fields obtained (Fig. 10) are qualitatively compared with the available geological and geophysical information (Figs. 8 and 9). This highlights
<table><caption>**Table 3**
Summary of optimization results.</caption><thead><tr><th>Optimization results</th><th>PEST++GLM</th><th>PEST++IES</th></tr></thead><tbody><tr><td>Number of PEST iterations</td><td>30</td><td>3</td></tr><tr><td>Number of model executions</td><td>9954</td><td>200</td></tr><tr><td>Computing time (hours)</td><td>43.8</td><td>1.27</td></tr><tr><td>RMSE (m)</td><td>0.75</td><td>Range: 1.12 – 3.26<br>Mean: 2.31</td></tr></tbody></table>

the passage of a zone of low hydraulic conductivity where the fault has been described and observed in geophysical investigation. The results obtained for the hydraulic conductivity field suggest that the Rieu Coullon fault zone may have a less rectilinear direction than that indicated on the geological map. The observed hydraulic conductivity contrasts highlight the parts of the fault zone that are impermeable. These parts can be interpreted as CZ lenses (Géraud et al., 2006; Davis et al., 2011).

Fig. 10B and 10D illustrate the uncertainties of the model

permeabilities where a better estimation occurs within the area between and around wells. IES estimation results show lower uncertainty values than those obtained with FOSM estimation from GLM.

### 4.4. Hydrodynamic interpretation of the Rieu Coullon fault

The application of inverse estimation of hydraulic conductivity on the fault zone of the Rieu Coullon site yields a range of values for this hydrodynamic property, ranging from $1.10^{-3}$ m s$^{-1}$ to $1.10^{-9}$ m s$^{-1}$. The hydrodynamic estimation has been made on parameters estimated by IES due to the assumed best estimate of uncertainties and knowing that the results obtained by GLM and IES are concordant.

Three cross-sectional profiles across the Rieu Coullon Fault are depicted in Fig. 11. The location of these profiles was deliberately selected based on the hydraulic conductivity field estimated by IES and the areas with the lowest estimation uncertainties (Fig. 10D). The hydraulic conductivity variations obtained for the Rieu Coullon Fault and those derived from hydrodynamic knowledge in the fault zone, as

<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Comparison of hydraulic conductivity estimates in the Rieu Coullon fault zone using GLM and IES. (A) and (C) display the estimated hydraulic conductivity fields obtained with GLM and IES, respectively. It should be noted that the scales of the hydraulic conductivity fields differ between GLM and IES. (B) and (D) depict the maps of uncertainties associated with the hydraulic conductivity field estimates using GLM and IES, respectively.</figcaption></figure>
illustrated in Fig. 11A, 11B, 11C and in Fig. 2, are consistent. The hy-
draulic conductivity value profiles in Fig. 11 demonstrate significant
variations along the Rieu Coullon fault. In the peripheral zones of the CZ
shape, higher hydraulic conductivity values around 10⁻³ m s⁻¹ are
observed, which can be associated with a DZ where thickness varies
along the fault.

The zone ranges from low hydraulic conductivity (10⁻⁷ – 10⁻⁸ m s⁻¹) to impermeable values (10⁻⁸ – 10⁻⁹ m s⁻¹) (Fig. 11). Therefore, the core zone could be composed of impermeable fault rock such as cataclasite or gouge to exhibit such values (Sibson, 1977). The estimated orders of magnitude of hydraulic conductivity indicate the presence of a well-developed fault zone in the Rieu Coullon area, characterized by a continuous impermeable CZ and a highly permeable DZ. This observation suggests the successful development of a fracture network, forming a conduit-barrier fault zone (Micarelli et al., 2003; Micarelli et al., 2006).

On the survey area, the thickness of the impermeable zone varies significantly, ranging from 1.5 m to 20 m. The interpreted thickness of the fault zone based on the hydraulic conductivity field reveals a width of approximately 65 m. These dimensions of fault zone thickness correspond to on-site observations indicating altered and fractured zones within this interval (Clauzon, 2019). This width demonstrates a coherent scaling relationship with the apparent fault offset (50 m), which aligns with the linear scaling relationship observed in previous research by Mayolle et al. (2019) in a carbonate environment.

The approach carried out at the Rieu Coullon site has enabled the
expansion of point observations and their utilization for a spatial
interpretation of the fault zone.

The application of this developed approach to the fault zone of Rieu
Coullon highlights the significant heterogeneity of the fault zone, as well
as variations in the thickness of the CZ and DZ structures along the fault.

4.5. Case study forecasts and uncertainty

As this aquifer is highly sensitive, we are seeking to assess the hy-
drodynamic impact of this type of fault zone on the water resource in
terms of point-source pollution. For this, travel time forecasts have been
performed between wells 3 and 1 with both FOSM from GLMA estima-
tion and IES. Results show a mean travel time of about 300 days for both
methods (Fig. 12). Both differences between prior and posterior travel
time confirm the influence of an impermeable CZ of the Rieu Coullon
fault, as evidenced by long-term transport times.

In contrast to the synthetic case, the range of possible results esti-
mated by IES extends over a broader range of times compared to the
estimation by GLMA with FOSM, highlighting the greater capability of
IES to explore the wide null space of the inverse problem.

The comparison of prior and posterior distributions demonstrates the
effectiveness of IES and GLMA with FOSM in reducing predictive
uncertainty.

<figure><img src="image_15.png" /><figcaption>Fig. 11. Hydrodynamic interpretation: hydraulic conductivity perpendicular to the Rieu Coullon fault zone. A, B, and C: Trends in hydraulic conductivity along the northern (A), central (B), and southern (C) segments of the fault. In red is the average variation in hydraulic conductivity, and in grey lines is the set of hydraulic conductivity estimation models obtained with IES. D: Map of hydraulic conductivity estimated by IES and interpretation of CZ location. DZ and CZ represent the damage zone and the multicore zone. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Results of travel time prediction estimates on both sides of the fault zone on the Rieu Coullon Fault. Forecast uncertainty reduction with GLM (A) and IES (B). Time data before estimation are shown in grey, and time data after the reduction of uncertainties are shown in blue.</figcaption></figure>

## 4.6. Contributions of case study for karstic fault zones

The hydrodynamic interpretations derived from the estimation results suggest that the Rieu Coullon fault exhibits conduit/barrier behavior, supported by flow simulations. Groundwater flow appears to be preferentially directed through the damage zone along a parallel axis. There is a potential for flow channelization parallel to the core zone but through the damage zone. This dynamic may facilitate karst

development. The hydrodynamic interpretation of the Rieu Coullon fault reveals the presence of a highly permeable damage zone, making it highly vulnerable to pollution infiltration. This serves as an additional example of a fault zone in fractured carbonate environments (Micarelli et al., 2006; Hausegger et al., 2010; Bauer et al., 2016; Geng et al., 2021), contributing to our understanding of hydrodynamic properties in such complex settings.

This study provides an original proposal for a distribution of

<figure><img src="image_17.png" /><figcaption>**Fig. 13.** A: Conceptual diagram of the internal structures of a karstified carbonate fault zone. B: Evolution of hydraulic conductivity perpendicular to the central segment of the karstified fault zone. The red line represents the average variations in hydraulic conductivity resulting from structural heterogeneity and karstification. The dotted black line indicates the actual hydraulic conductivity, considering the extreme variations in hydraulic conductivity imposed by karstification and multiple core zones.</figcaption></figure>
hydraulic conductivities within the CZ and DZ, representative of the dynamics of groundwater flows and consistent with hydrological observations (Fig. 13). This application provides an initial insight into the magnitude variations encountered in these karst environments (Fig. 11A, 11B, 11C). The hydraulic conductivity variations illustrated in Fig. 13 in this case study range between $10^{-3}$ and $10^{-9}$ m s$^{-1}$, with values of $10^{-3}$ m s$^{-1}$ in the DZ, $10^{-7}$ to $10^{-9}$ m s$^{-1}$ in the CZ, and K values of $10^{-5}$ m s$^{-1}$ for the healthy protolith outside the fault zone. The extremely low values encountered in the CZ indicate that even in the context of rock dissolution in a fractured network, the fault rock in the CZ can remain intact and form an impermeable barrier. On the other hand, the K values encountered in the DZ indicate high permeability. In our choice of equivalent porous model, we average the hydrodynamic effects of the karst network and the matrix. In this case study, it is possible to hypothesize that the high K values may be associated with the karst network. This is further supported by the well-known fact that fractures and areas with high fracture density, like those found in the DZ, contribute to the development of the karst network (Lopes et al., 2022; Furtado et al., 2022).

## 5. Discussion and conclusion

The proposed approach to estimate the hydraulic conductivity field from cross-hole pumping tests for fault zone hydrodynamic simulation has proven its efficiency for a conduit-barrier fault zone of a carbonated and karstified environment with a sparse spatial dataset. The synthetic test case that validated the approach has shown that the use of solely hydraulic head measurements from cross-hole pumping tests estimates both the direction and the hydraulic conductivity field of the fault zone. Head and pollutant travel time forecasts are close to the “true” value for the synthetic case, and the associated posterior predictive uncertainty is greatly reduced from the prior. It is important to highlight that even the travel time predictions are close to the real value and the posterior is significantly improved, even though no transport information was incorporated for parameter estimation. This finding needs to be qualified, as mass transfer has not been observed and may be largely influenced by the karstic nature of the aquifer.

The observation of groundwater levels on either side of a fault zone can lead to a binary interpretation, either conduit or barrier, ultimately providing no insight into the permeability contrasts between the damage zone and the core zone (Berkowitz, 2002). However, by processing this information using the methodology presented, the possibility to achieve a interesting level of detail on the spatial heterogeneity of the fault zone has been shown, with precision to distinguish between the behavior of the damage zone and that of the core zone.

Our results, providing a spatial interpretation of fault zone hydrodynamics, demonstrate that the head variation data are effective for understanding fault zone hydrodynamic structures. The hydrodynamic properties of the fault zone are estimated directly from observation influenced by them, i.e. groundwater flow within the fault. We aim to encourage future hydrodynamic interpretations of fault zones to incorporate hydrological data more extensively in addition to structural observations. In fractured carbonate environments, the development of preferential karstic networks along fractures is a complex element to anticipate based solely on surface geological observations (Bonetto et al., 2008; De Waele et al., 2011; Parise et al., 2015). Depending on numerous factors like the topological characteristics of the fracture networks (Zimmerman and Main, 2004) or the type of fracture contact (Berkowitz et al., 1994), the development of the karstic network introduces significant uncertainties when its only morphologically estimate (Bakalowicz, 2005). Simply using groundwater levels as the main indicator of the hydrodynamic behavior of a fault zone allows to encompass flow through matrix, within karst conduits, across fractured media, and through fault zone structures. The hydrodynamic interpretation validation conducted on a real fault zone has been made possible by the findings, shedding light on the dimensions and internal structural

entities of the fault zone. This interpretation is further supported by results consistent with the geological and structural characteristics of faults in carbonate formations (Micarelli et al., 2003, 2006; Balsamo et al., 2016; Maqbool et al., 2016; Clauzon et al., 2020; Mayolle et al., 2019). By averaging the effects of the karst network and the matrix using the equivalent porous model to estimate the conductivity field, we provide a representation of the hydrodynamic behaviors of these two karst entities. Ultimately, this approach allows the obtained hydraulic conductivity field to link high levels of permeability to a potential influence of the karst network in these areas.

Our results offer a hydrodynamic estimation of the fault zone, overcoming the challenge of limited structural and geological information. This methodology addresses the issue of accessing geological data in fault zones and underscores the significance of hydrological information in enhancing our understanding of fault zones.

This approach improves decision-making on quantitative and qualitative issues for better management of water resources as it has been applied with success to a real-world case study where the results of the inversion match the geophysical and geological observations.

The proposed approach is designed to be easily replicable even with limited information of a fault zone, and with a limited number of on both sides of a fault zone. It is recommended to compare the results with complementary qualitative and/or quantitative observations such as geological and hydrogeological studies of the fault zone (Micarelli et al., 2006; Lockner et al., 2009; Ball et al., 2010; Bense et al., 2013). Presented parametric uncertainty quantification can help to find areas where observations are missing, and a data worth methodology could be used for better spatial implementation of complementary observations.

There are a few recommendations for the proper use of this method.

The estimation of hydraulic conductivity faces limitations when observations are confined to one side of the fault or when spatial information around the fault zone is scarce. As a result, this methodology requires access to hydrodynamic information from wells located on both sides of the fault zone. Additionally, this approach does not consider vertical variations in hydraulic conductivity. In cases where lithologies with differing permeabilities are intersected by wells, our method will be constrained by this vertical heterogeneity. The aim of this method is not to make a fine horizontal and vertical cartography of the hydraulic conductivity along the fault zone, but to simulate the flow and transport within it. The equivalent porous model works well in this respect, but a larger number of wells in the different zones of the fault would enable a finer hydraulic tomography. Although the equivalent porous medium model simplifies the approach and improves estimations and simulations by averaging the effects of the karst network and fault zone structures on estimated permeability, transport simulations may underestimate travel times due to the absence of a fracture network model.

Overcoming ill-posed problems necessitates the use of regularization, such as the preferred parameter Tikhonov regularization. However, for fault zone estimation purposes we aim to characterize a highly contrasted hydraulic field where the regularization has the effect of smoothing the field. The trade-off between measurement objective function and regularization objective function is justified by the desire to achieve a significant reduction of the objective measurement smooth enough while avoiding excessive adjustment that can lead to an overfitting with an unrealistic hydraulic conductivity field. One important aspect of the presented method is its ability to incorporate additional site-specific hard or soft knowledge to constrain the inverse model. If localized hydrodynamic property data is available, it could be directly informed through the pilot points, and the Tikhonov regularization can be used to enhance the model by incorporating the information.

A few points could be considered in the future. Diverse data types could be incorporated as observation through the inversion process. For example, observations of artificial or natural tracer tests or electrical observations with coupled hydro-geophysical modeling can improve parameter estimation and simulation of fault zones (Irving and Singha 2010; Cousquer et al., 2018). The use of vertically discretized head and/
or transport observation coupled with packer tests could extend this type of approach to 3 dimensions. The integration of information on transport times between boreholes in our application could significantly enhance particle tracking estimates for advective transport. This could better constrain posterior time estimates in the graphs of Fig. 12 and reduce the spread of models around a more probable value.

To further enhance the understanding of hydrodynamics within fault zones, we could investigate the scaling relationship between static observations such as structural geological information, and hydrodynamic observations on a large set of fault zones to upscale the hydrodynamic function of fault zones in a given context. From this kind of future work, knowledge along the fault zone area at the basin scale can be implemented by linking pilot points along the fault zone and assigning prior knowledge of hydraulic conductivity values.

This study also enabled us to compare the GLM and IES methods in a certain context. GLM needs a much higher computational time demand than IES due to the calculation of the Jacobian matrix at each algorithm iteration. Compared to GLM, which is a deterministic method, IES is an ensemblistic method based on the parameter space sampling that leads to a set of realizations that is particularly suitable for better uncertainty analysis and forecast purposes. IES is well suited for investigating the diverse range of possible solutions for highly parametrized models and capturing the variability inherent to the estimation process with a particularly short computation time. GLM has demonstrated its well-known capability for both estimation and uncertainty quantification close to the IES results. The two inverse problem-solving methods used demonstrate significantly different needs in terms of the number of iterations required. GLM requires a number of iterations primarily dependent on the reduction of the objective function. In contrast, IES significantly reduces the objective function within the first few iterations, necessitating a limit on the number of iterations to prevent the ensemble of models from converging into a single unique model. Furthermore, the computational capabilities of IES enable the extension of this strategy to larger scales, such as the aquifer scale, by incorporating a comprehensive dataset of geological and hydrogeological information and considering a wide range of multidisciplinary data set.

A reproducible workflow for flow and transport simulation of fault zones with parameter estimation and predictions under uncertainty has been proposed. We successfully demonstrated the efficiency of the methodology for a conduit-barrier fault zone in a moderately karstified carbonate environment, first on a synthetic case and then in a real-world case study application. This study is innovative in that this 2D distribution of hydrodynamic properties is proposed on a multi-decametric scale of the fault, rather than a multi-metric scale (geological approach) or a hecto-metric scale (reservoir approach). In terms of the dynamics of underground flow in a fault zone, it can be used as a basis for the reproducibility of hydrogeological models in faulted aquifers. Application of this approach to a real fault zone in a carbonate environment provides comprehensive insights into karstic environments within hydrodynamic studies of fault zones, thereby improving qualitative and quantitative groundwater management.

## CRediT authorship contribution statement

**Aurélie Boura:** Conceptualization, Methodology, Investigation, Data curation, Software, Writing – original draft, Visualization. **Yohann Cousquer:** Conceptualization, Methodology, Investigation, Validation, Software, Writing – review & editing, Visualization. **Victor Clauzon:** Data curation. **Rémi Valois:** Data curation, Visualization. **Véronique Leonardi:** Funding acquisition, Resources, Investigation, Supervision, Writing – review & editing, Visualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence

the work reported in this paper.

## Acknowledgments

This work has been funded by a French government grant managed by the Agence Nationale de la Recherche under the France 2030 program, reference ANR-22-PEXO-0012, as part of the OneWater - Eau Bien Commun national research program. We are grateful to the KARST Observatory Network, which aims to enhance knowledge-sharing and promote cross-disciplinary research on karst systems within the framework of the OZCAR Critical Zone network Research Infrastructure. This network includes the MEDYDYSS observatory and enables the instrumentation of the Lez aquifer. Special thanks to Clotilde Thompson for her assistance with the English. We thank John Doherty for constructive discussions and valuable insights. We also thank the two anonymous reviewers and Editor C. Corradini for their suggestions, which greatly helped improve the manuscript.

## References

Abusaada, M., Sauter, M., 2013. Studying the flow dynamics of a karst aquifer system with an equivalent porous medium model. Groundwater 51, 641–650. https://doi.org/10.1111/j.1745-6584.2012.01003.x.

Agosta, F., Kirschner, D.L., 2003. Fluid conduits in carbonate-hosted seismogenic normal faults of central Italy. J. Geophys. Res. Solid Earth 108. https://doi.org/10.1029/2002JB002013.

Agosta, F., Prasad, M., Aydin, A., 2007. Physical properties of carbonate fault rocks, fucino basin (Central Italy): implications for fault seal in platform carbonates. Geofluids 7, 19–32. https://doi.org/10.1111/j.1468-8123.2006.00158.x.

Alarcon, V.J., Linhoss, A.C., Mickle, P.F., Kelble, C.R., Fine, A., 2022. Estimation of Groundwater Salinity for the Central Biscayne Bay Coast, Florida, USA: in: Gervasi, O., Murgante, B., Misra, S., Rocha, A.M.A.C., Garau, C. (Eds.), Computational Science and Its Applications – ICCSA 2022 Workshops. Springer International Publishing, Cham, pp. 594–606. https://doi.org/10.1007/978-3-031-10545-6_40.

Anderson, T.R., Fairley, J.P., 2008. Relating permeability to the structural setting of a fault-controlled hydrothermal system in southeast Oregon, USA. J. Geophys. Res. 113, B05402. https://doi.org/10.1029/2007JB004962.

Anderson, M.P., Woessner, W.W., Hunt, R.J., 2015. Chapter 10 - Forecasting and Uncertainty Analysis. In: Anderson, M.P., Woessner, W.W., Hunt, R.J. (Eds.), Applied Groundwater Modeling (second Edition). Academic Press, San Diego, pp. 443–491. https://doi.org/10.1016/B978-0-08-091638-5.00010-9.

Andre, B.J., Rajaram, H., 2005. Dissolution of limestone fractures by cooling waters: early development of hypogene karst systems. Water Resour. Res. 41. https://doi.org/10.1029/2004WR003331.

Antonellini, M., Aydin, A., 1994. Effect of faulting on fluid flow in porous sandstones: petrophysical properties1. AAPG Bull. 78, 355–377. https://doi.org/10.1306/BDFF90AA-1718-11D7-8645000102C1865D.

Aster, R.C., Borchers, B., Thurber, C.H., 2013. Chapter Nine - Nonlinear Regression. In: Aster, R.C., Borchers, B., Thurber, C.H. (Eds.), Parameter Estimation and Inverse Problems (second Edition). Academic Press, Boston, pp. 217–238. https://doi.org/10.1016/B978-0-12-385048-5.00009-4.

Aster, R.C., Borchers, B., Thurber, C.H., 2018. Parameter Estimation and Inverse Problems. Elsevier.

Bakalowicz, M., 2005. Karst groundwater: a challenge for new resources. Hydrogeol J 13, 148–160. https://doi.org/10.1007/s10040-004-0402-9.

Bakker, M., Post, V., Langevin, C.D., Hughes, J.D., White, J.T., Starn, J.J., Fienen, M.N., 2016. Scripting MODFLOW model development using python and FloPy. Groundwater 54, 733–739. https://doi.org/10.1111/gwat.12413.

Ball, L.B., Ge, S., Caine, J.S., Revil, A., Jardani, A., 2010. Constraining fault-zone hydrogeology through integrated hydrological and geoelectrical analysis. Hydrogeol J 18, 1057–1067. https://doi.org/10.1007/s10040-010-0587-z.

Balsamo, F., Clemenzi, L., Storti, F., Mozafari, M., Solum, J., Swennen, R., Taberner, C., Tueckmantel, C., 2016. Anatomy and paleofluid evolution of laterally restricted extensional fault zones in the Jabal Qusaybah anticline, Salakh arch, Oman. Geol. Soc. Am. Bull. 128, 957–972. https://doi.org/10.1130/B31317.1.

Balsamo, F., Clemenzi, L., Storti, F., Solum, J., Taberner, C., 2019. Tectonic control on vein attributes and deformation intensity in fault damage zones affecting Natih platform carbonates, Jabal Qusaybah, North Oman. J. Struct. Geol. 122, 38–57. https://doi.org/10.1016/j.jsg.2019.02.009.

Balsamo, F., Storti, F., 2010. Grain size and permeability evolution of soft-sediment extensional sub-seismic and seismic fault zones in high-porosity sediments from the Crotone basin, southern Apennines, Italy. Mar. Pet. Geol. 27, 822–837. https://doi.org/10.1016/j.marpetgeo.2009.10.016.

Balsamo, F., Storti, F., Salvini, F., Silva, A.T., Lima, C.C., 2010. Structural and petrophysical evolution of extensional fault zones in low-porosity, poorly lithified sandstones of the Barreiras Formation, NE Brazil. Journal of Structural Geology, Fault Zones 32, 1806–1826. https://doi.org/10.1016/j.jsg.2009.10.010.
Bastesen, E., Braathen, A., 2010. Extensional faults in fine grained carbonates - analysis of fault core lithology and thickness-displacement relationships. *J. Struct. Geol.* 32, 1609-1628. https://doi.org/10.1016/j.jsg.2010.09.008.

Bauer, J.F., Meier, S., Philipp, S.L., 2015. Architecture, fracture system, mechanical properties and permeability structure of a fault zone in Lower Triassic sandstone, Upper Rhine Graben. *Tectonophysics* 647-648, 132-145. https://doi.org/10.1016/j.tecto.2015.02.014.

Bauer, H., Schröckenfuchs, T.C., Decker, K., 2016. Hydrogeological properties of fault zones in a karstified carbonate aquifer (Northern Calcareous Alps, Austria). *Hydrogeol. J.* 24, 1147-1170. https://doi.org/10.1007/s10040-016-1388-9.

Bense, V.F., Gleeson, T., Loveless, S.E., Bour, O., Scibek, J., 2013. Fault zone hydrogeology. *Earth Sci. Rev.* 127, 171-192. https://doi.org/10.1016/j.earscirev.2013.09.008.

Berkowitz, B., 2002. Characterizing flow and transport in fractured geological media: a review. *Adv. Water Resour.* 25, 861-884. https://doi.org/10.1016/S0309-1708(02)00042-8.

Berkowitz, B., Naumann, C., Smith, L., 1994. Mass transfer at fracture intersections: an evaluation of mixing models. *Water Resour. Res.* 30, 1765-1773. https://doi.org/10.1029/94WR00432.

Bigi, S., Battaglia, M., Alemanni, A., Lombardi, S., Campana, A., Borisova, E., Loizzo, M., 2013. CO2 flow through a fractured rock volume: insights from field data, 3D fractures representation and fluid flow modeling. *Int. J. Greenhouse Gas Control* 18, 183-199. https://doi.org/10.1016/j.iggc.2013.07.011.

Billi, A., Salvini, F., Storti, F., 2003. The damage zone-fault core transition in carbonate rocks: implications for fault growth, structure and permeability. *J. Struct. Geol.* 25, 1779-1794. https://doi.org/10.1016/S0191-8141(03)00037-3.

Bonetto, S., Fuorucci, A., Fornaro, M., Vigna, B., 2008. Subsidence hazards connected to quarrying activities in a karst area: the case of the Moncalvo sinkhole event (Piedmont, NW Italy). *Estonian Journal of Earth Sciences* 57. https://doi.org/10.3176/earth.2008.3.01.

Burrows, W., Doherty, J., 2016. Gradient-based model calibration with proxy-model assistance. *J. Hydrol.* 533, 114-127. https://doi.org/10.1016/j.jhydrol.2015.11.033.

Caine, J.S., Evans, J.P., Forster, C.B., 1996. Fault zone architecture and permeability structure. *Geol. 24*, 1025. https://doi.org/10.1130/0091-7613(1996)024<1025: FZAAPS-2.3.C02.

Carrera, J., Alcolea, A., Medina, A., Hidalgo, J., Slooten, L.J., 2005. Inverse problem in hydrogeology. *Hydrogeol J* 13, 206-222. https://doi.org/10.1007/s10040-004-0404-7.

Chen, Y., Oliver, D.S., 2013. Levenberg-Marquardt forms of the iterative ensemble smoother for efficient history matching and uncertainty quantification. *Comput Geosci* 17, 689-703. https://doi.org/10.1007/s10596-013-9351-5

Cherry, J.A., McWorther, D.B., Parker, B.L., 2009. Site conceptual model for the migration and fate of contaminants in groundwater at the Santa Susana Field Laboratory, Simi, California (draft), vols 1-4. Association with the University of Guelph, Toronto, ON; MWH, Walnut Creek, CA.

Childs, C., Walsh, J.J., Manzocchi, T., Strand, J., Nicol, A., Tomasso, M., Schöpfer, M.P. J., Aplin, A.C., 2007. Definition of a fault permeability predictor from outcrop studies of a faulted turbidite sequence, Taranaki, New Zealand. *Geol. Soc. Lond. Spec. Publ.* 292, 235-258. https://doi.org/10.1144/SP292.14.

Childs, C., Manzocchi, T., Walsh, J.J., Bonson, C.G., Nicol, A., Schöpfer, M.P.J., 2009. A geometric model of fault zone and fault rock thickness variations. *J. Struct. Geol.* 31, 117-127. https://doi.org/10.1016/j.jsg.2008.08.009.

Clarkson, C.R., Vahedian, A., Ghanizadeh, A., Song, C., 2019. A new low-permeability reservoir core analysis method based on rate-transient analysis theory. *Fuel* 235, 1530-1543. https://doi.org/10.1016/j.fuel.2018.07.115.

Clauzon, V., Mayolle, S., Leonardi, V., Brunet, P., Soliva, R., Marchand, P., Massonnat, G., Rolando, J.-P., Pistre, S., 2020. Fault zones in limestones: impact on karstogenesis and groundwater flow (Lez aquifer, southern France). *Hydrogeol. J.* 28, 2387-2408. https://doi.org/10.1007/s10040-020-02189-9.

Clauzon, V., 2019. Caractérisation en situ multi-échelles des transferts de fluide en zone de faible en milieu carbonaté.

Collon, P., Bernasconi, D., Vuilleumier, C., Renard, P., 2017. Statistical metrics for the characterization of karst network geometry and topology. *Geomorphology* 283, 122-142. https://doi.org/10.1016/j geomorph.2017.01.034.

Commer, M., Kowalsky, M.B., Doetsch, J., Newman, G.A., Finsterle, S., 2014. MPITOUGH2: A parallel parameter estimation framework for hydrological and hydrogeophysical applications. *Computers & Geosciences*, *TOUGH Symposium* 2012 (615), 127-135. https://doi.org/10.1016/j.cageo.2013.06.011.

Coulon, C., Lemieux, J.-M., Pryet, A., Bayer, P., Young, N.L., Molson, J., 2022. Pumping optimization under uncertainty in an island freshwater lens using a sharp-interface seawater intrusion model. *Water Resources Research* 58, e2021WR031793. https://doi.org/10.1029/2021WR031793.

Cousquer, Y., Jourde, H., 2022. Reducing uncertainty of karst aquifer modeling with complementary hydrological observations for the sustainable management of groundwater resources. *J. Hydrol.* 612, 128130. https://doi.org/10.1016/j.jhydrol.2022.128130.

Cousquer, Y., Pryet, A., Atteia, O., Ferré, T.P.A., Delbart, C., Valois, R., Dupuy, A., 2018. Developing a particle tracking surrogate model to improve inversion of ground water - surface water models. *J. Hydrol.* 558, 356-365. https://doi.org/10.1016/j.jhydrol.2018.01.043.

Davis, G.H., Reynolds, S.J., Kluth, C.F., 2011. Structural Geology of Rocks and Regions. *John Wiley & Sons*

Day-Lewis, F.D., Hsieh, P.A., Gorelick, S.M., 2000. Identifying fracture-zone geometry using simulated annealing and hydraulic-connection data. *Water Resour. Res.* 36, 1707-1721. https://doi.org/10.1029/2000WR900073.

De Marsily, G., Lavedan, G., Boucher, M., Fasamino, G., 1984. Interpretation of interference tests in a well field using geostatistical techniques to fit the permeability distribution in a reservoir model. In: Presented at the Geostatistics for Natural Resources Characterization. NATO Advanced Study Institute, pp. 831-849.

de Souza, F.M., Gomes, I.F., Nogueira, F.C.C., Vasconcelos, D.L., Canabarro, B., de Souza, J.A.B., Guimarães, L.J.do.N., Beserra, L.B.de.S., 2022. 2D modeling and simulation of deformation bands' effect on fluid flow: Implications for hydraulic properties in siliciclastic reservoirs. *J. Struct. Geol.* 158, 104581. https://doi.org/10.1016/j.jsg.2022.104581.

De Waele, J., Gutiérrez, F., Parise, M., Plan, L., 2011. Geomorphology and natural hazards in karst areas: a review. *Geomorphology* 134, 1-8. https://doi.org/10.1016/j geomorph.2011.08.001.

Delottier, H., Therrien, R., Young, N.L., Paradis, D., 2022. A hybrid approach for integrated surface and subsurface hydrologic simulation of baseflow with Iterative Ensemble Smoother. *J. Hydrol.* 606, 127406. https://doi.org/10.1016/j.jhydrol.2021.127406.

Doherty, J., 2003. Ground water model calibration using pilot points and regularization. *Groundwater* 41, 170-177. https://doi.org/10.1111/j.1745-6584.2003.tb02580.x.

Doherty, J., 2015. Calibration and uncertainty analysis for complex environmental models. *Groundwater* 53, 673-674. https://doi.org/10.1111/gwat.12360.

Doherty, J., Hunt, R., 2010. Approaches to highly parameterized inversion—A guide to using PEST for groundwater-model calibration: U.S. (Open-File Report No. 5169), Geological Survey Scientific Investigations Report 2010.

Doherty, J., Welter, D., 2010. A short exploration of structural noise. *Water Resour. Res.* 46. https://doi.org/10.1029/2009WR008377.

Doherty, J., 2011. Modeling: Picture Perfect or Abstract Art? Ground Water 49, 455-455. https://doi.org/10.1111/j.1745-6584.2011.00812.x.

Dörfliger, N., Fleury, P., Ladouche, B., 2009. Inverse modeling approach to allogenic karst system characterization. *Groundwater* 47, 414-426. https://doi.org/10.1111/j.1745-6584.2008.00517.x.

Drogue, C., 1972. Analyse statistique des hydrogrammes de decymes des sources karstiques statistical analysis of hydrographs of karstic springs. *J. Hydrol.* 15, 49-68. https://doi.org/10.1016/0022-1694(72)90075-3.

Essouayed, E., Ferré, T., Cohen, G., Guiserix, N., Atteia, O., 2021. Application of an iterative source localization strategy at a chlorinated solvent site. *Journal of Hydrology* X 13, 100111. https://doi.org/10.1016/j.hydroa.2021.100111.

Evans, J.P., Forster, C.B., Goddard, J.V., 1997. Permeability of fault-related rocks, and implications for hydraulic structure of fault zones. *J. Struct. Geol.* 19, 1393-1404. https://doi.org/10.1016/S0191-8141(97)00057-6.

Faoro, I., Niemeijer, A., Marone, C., Elsworth, D., 2009. Influence of shear and deviatoric stress on the evolution of permeability in fractured rock. *J. Geophys. Res. Solid Earth* 114. https://doi.org/10.1029/2007JB005372.

Faulkner, D.R., 2004. A model for the variation in permeability of clay-bearing fault gouge with depth in the brittle crust. *Geophys. Res. Lett.* 31. https://doi.org/10.1029/2004GL020736.

Faulkner, D.R., Jackson, C.A.L., Lunn, R.J., Schliche, R.W., Shipton, Z.K., Wibberley, C. A.J., Withjack, M.O., 2010. A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones. *Journal of Structural Geology, Fault Zones* 32, 1557-1575. https://doi.org/10.1016/j.jsg.2010.06.009.

Ferraro, F., Agosta, F., Prasad, M., Vinciguerra, S., Violay, M., Giorgioni, M., 2020. Pore space properties in carbonate fault rocks of peninsular Italy. *J. Struct. Geol.* 130, 103913. https://doi.org/10.1016/j.jsg.2019.103913.

Fischer, P., Jardani, A., Soueid Ahmed, A., Abbas, M., Wang, X., Jourde, H., Lecoq, N., 2017. Application of large-scale inversion algorithms to hydraulic tomography in an alluvial aquifer. *Groundwater* 55, 208-218. https://doi.org/10.1111/gwat.12457.

Fischer, P., Jardani, A., Lecoq, N., 2018. Hydraulic tomography of discrete networks of conduits and fractures in a karstic aquifer by using a deterministic inversion algorithm. *Adv. Water Resour.* 112, 83-94. https://doi.org/10.1016/j.adwatres.2017.11.029.

Fisher, Q.J., Haneef, J., Grattoni, C.A., Allshorn, S., Lorinci, P., 2018. Permeability of fault rocks in siliciclastic reservoirs: recent advances. *Mar. Pet. Geol.* 91, 29-42. https://doi.org/10.1016/j.marpetgeo.2017.12.019.

Ford, D.C., Ewers, R.O., 1978. The development of limestone cave systems in the dimensions of length and depth. *Can. J. Earth Sci.* 15, 1783-1798. https://doi.org/10.1139/e78-186.

Ford, D., Williams, P., 2007. Karst hydrogeology and geomorphology. *John Wiley & Sons* Ltd, West Sussex, England. https://doi.org/10.1002/9781118649849.

Fredman, N., Tveranger, J., Semshaug, S., Braathen, A., Sverdrup, E., 2007. Sensitivity of fluid flow to fault core architecture and petrophysical properties of fault rocks in siliciclastic reservoirs: A synthetic fault model study. *Petroleum Geoscience* - *Petrol Geosci.* 13, 305-320. https://doi.org/10.1144/1354-079306-721.

Furtado, C.P.Q., Medeiros, W.E., Borges, S.V., Lopes, J.A.G., Bezerra, F.H.R., Lima-Filho, F.P., Maia, R.P., Bertotti, G., Auler, A.S., Teixeira, W.L.E., 2022. The influence of subseismic-scale fracture interconnectivity on fluid flow in fracture corridors of the Brejões carbonate karst system. *Brazil. Mar. Pet. Geol.* 141, 105689. https://doi.org/10.1016/j.marpetgeo.2022.105689.

Gassiat, C., Gleeson, T., Lefebvre, R., McKenzie, J., 2013. Hydraulic fracturing in faulted sedimentary basins: numerical simulation of potential contamination of shallow aquifers over long time scales. *Water Resour. Res.* 49, 8310-8327. https://doi.org/10.1002/2013WR014287.

Geng, F., Wang, H., Hao, J., Gao, P., 2021. Internal Structure Characteristics and Formation Mechanism of Reverse Fault in the Carbonate Rock, A Case Study of Outcrops in Xike'er Area, Tarim Basin, Northwest China. *Frontiers in Earth Science* 9.

Géraud, Y., Diraison, M., Orellana, N., 2006. Fault zone geometry of a mature active normal fault: a potential high permeability channel (Pirgaki fault, Corinth rift,
Greece). Tectonophysics, Natural Laboratories on Seismogenic Faults 426, 61-76. https://doi.org/10.1016/j.tecto.2006.02.023.

Giese, M., Reimann, T., Bailly-Comte, V., Maréchal, J.-C., Sauter, M., Geyer, T., 2018. Turbulent and laminar flow in karst conduits under unsteady flow conditions: interpretation of pumping tests by discrete conduit-continuum modeling. Water Resour. Res. 54, 1918-1933. https://doi.org/10.1002/2017WR020658.

Haines, T.J., Michie, E.A.H., Neilson, J.E., Healy, D., 2016. Permeability evolution across carbonate hosted normal fault zones. Mar. Pet. Geol. 72, 62-82. https://doi.org/10.1016/j.marpetgeo.2016.01.008.

Hausegger, S., Kurz, W., Rabitsch, R., Kiechl, E., Brosch, F.-J., 2010. Analysis of the internal structure of a carbonate damage zone: Implications for the mechanisms of fault breccia formation and fluid flow. Journal of Structural Geology, Faulting and Fracturing of Carbonate Rocks: New Insights into Deformation Mechanisms, Petrophysics and Fluid Flow Properties 32, 1349-1362. https://doi.org/10.1016/j.jgsl.2009.04.014.

Hill, M.C., Kavetski, D., Clark, M., Ye, M., Arabi, M., Lu, D., Foglia, L., Mehl, S., 2016. Practical use of computationally frugal model analysis methods. Groundwater 54, 159-170. https://doi.org/10.1111/gwat.12330.

Hunt, R.J., Doherty, J., Tonkin, M.J., 2007. Are models too simple? arguments for increased parameterization. Groundwater 45, 254-262. https://doi.org/10.1111/j.1745-6584.2007.00316.x.

Hunt, R.J., White, J.T., Duncan, L.L., Haugh, C.J., Doherty, J., 2011. Evaluating lower computational burden approaches for calibration of large environmental models. Groundwater 59, 788-798. https://doi.org/10.1111/gwat.13106.

Ilman, W.A., 2014. Hydraulic tomography offers improved imaging of heterogeneity in fractured rocks. Groundwater 52, 659-684. https://doi.org/10.1111/gwat.12119.

Irving, J., Singha, K., 2010. Stochastic inversion of tracer test and electrical geophysical data to estimate hydraulic conductivities. Water Resour. Res. 46. https://doi.org/10.1029/2009WR008340.

Iskandar, I., Koike, K., 2011. Distinguishing potential sources of arsenic released to groundwater around a fault zone containing a mine site. Environ Earth Sci 63, 595-608. https://doi.org/10.1007/s12665-010-0727-8.

Jaquet, O., Siegel, P., Klubertanz, G., Benabderrhamane, H., 2004. Stochastic discrete model of karstic networks. Adv. Water Resour. 27, 751-760. https://doi.org/10.1016/j.adwatres.2004.03.007.

Jeannin, P.-Y., Artigue, G., Butscher, C., Chang, Y., Charlier, J.-B., Duran, L., Gill, L., Hartmann, A., Johannet, A., Jourde, V., Kavousi, A., Liesch, T., Liu, Y., Lüthi, M., Malard, A., Mazzilli, N., Pardo-Igúzquiza, E., Thiéry, D., Reimann, T., Schuler, P., Wöhling, T., Wunsch, A., 2021. Karst modelling challenge 1: results of hydrological modelling. J. Hydrol. 600, 126508. https://doi.org/10.1016/j.jhydrol.2021.126508.

Jones, T.A., Detwiler, R.L., 2016. Fracture sealing by mineral precipitation: the role of small-scale mineral heterogeneity. Geophys. Res. Lett. 43, 7564-7571. https://doi.org/10.1002/2016GL069598.

Jourde, H., Flodin, E., Aydin, A., Durlofsky, L., Wen, X., 2002. Computing Permeability of Fault Zones in Eolian Sandstone from Outcrop Measurements. https://doi.org/10.1306/61EEDC4-173E-11D7-8645000102C1865D.

Kavousi, A., Reimann, T., Liedl, R., Raeisi, E., 2020. Karst aquifer characterization by inverse application of MODFLOW-2005 CFPv2 discrete-continuum flow and transport model. J. Hydrol. 587, 124922. https://doi.org/10.1016/j.jhydrol.2020.124922.

Kitanidis, P.K., 1997. Introduction to Geostatistics: Applications in Hydrogeology. Cambridge University Press.

Korneva, I., Tondi, E., Agosta, F., Rustichelli, A., Spina, V., Bitonte, R., Di Cuia, R., 2014. Structural properties of fractured and faulted Cretaceous platform carbonates, Merge Plateau (southern Italy). Mar. Pet. Geol. 57, 312-326. https://doi.org/10.1016/j.marpetgeo.2014.05.004.

Langevin, C.D., Provost, A.M., Panday, S., Hughes, J.D., 2022. Documentation for the MODFLOW 6 Groundwater Transport Model (Report No. 6-A61), Techniques and Methods. Reston, VA. https://doi.org/10.3133/tm6A61.

Lapperre, R.E., Bense, V.F., Kasse, C., van Balen, R.T., 2022. Temporal and spatial variability of cross-fault groundwater-level differences: the impact of fault-induced permeability reduction, precipitation and evapotranspiration. Hydrogeol. J. 30, 1233-1257. https://doi.org/10.1007/s10040-022-02465-w.

Larsen, B., Grunnaleite, I., Gudmundsson, A., 2010. How fracture systems affect permeability development in shallow-water carbonate rocks: An example from the Gargano Peninsula, Italy. Journal of Structural Geology, Faulting and Fracturing of Carbonate Rocks: New Insights into Deformation Mechanisms, Petrophysics and Fluid Flow Properties 32, 1212-1230. https://doi.org/10.1016/j.jgsl.2009.05.009.

Leaf, A.T., Fienen, M.N., 2022. Floppy: The Python Interface for MODFLOW. Groundwater 60, 710-712. https://doi.org/10.1111/gwat.13259.

Leonardi, V., Massonnat, G., Darey, J., Planteblat, C., Gal, C., 2016. Karst genesis modelling of a regional Mediterranean aquifer (Lez, France). In: 43rd IAH International Congress “groundwater and Society: 60 Years of IAH”.

Leray, S., de Dreuzy, J.-R., Bour, O., Bresciani, E., 2013. Numerical modeling of the productivity of vertical to shallowly dipping fractured zones in crystalline rocks. J. Hydrol. 481, 64-75. https://doi.org/10.1016/j.jhydrol.2012.12.014.

Li, Y., Sun, J., Wei, H., Song, S., 2019. Architectural features of fault-controlled karst reservoirs in the Tahe oilfield. J. Pet. Sci. Eng. 181, 106208. https://doi.org/10.1016/j petro.2019.106208.

Lien, J.-M., Liu, G., Langevin, 2015. GRIDGEN version 1.0—A computer program for generating unstructured finite-volume grids: U.S. (Open-File Report No. 1109), Geological Survey Open-File Report.

Littva, J., Hök, J., Bella, P., 2015. Cavitonics: Using caves in active tectonic studies (Western Carpathians, case study). J. Struct. Geol. 80, 47-56. https://doi.org/10.1016/j.jgsl.2015.08.011.

Lockner, D.A., Tanaka, H., Ito, H., Ikeda, R., Omura, K., Naka, H., 2009. Geometry of the Nojima Fault at Nojima-Hirabayashi, Japan – I. A simple damage structure inferred from borehole core permeability. Pure Appl. Geophys. 166, 1649-1667. https://doi.org/10.1007/s00024-009-0515-0.

Loke, M.H., 1999. Electrical imaging surveys for environmental and engineering studies. A Practical Guide to 2, 70.

Lopes, J.A.G., Medeiros, W.E., La Bruna, V., de Lima, A., Bezerra, F.H.R., Schozer, D.J., 2022. Advancements towards DFKN modelling: Incorporating fracture enlargement resulting from karstic dissolution in discrete fracture networks. J. Petrol. Sci. Eng. 29, 109944. https://doi.org/10.1016/j.petrol.2021.109944.

Lu, S., Chen, N., Hu, B., Cheng, J., 2012. On the inverse problems for the coupled continuum pipe flow model for flows in karst aquifers. Inverse Prob. 28, 065003. https://doi.org/10.1088/0266-5611/28/6/065003.

Maqbool, A.-R., Moustafa, A.R., Dowidar, H., Yousef, M., 2016. Architecture of fault damage zones of normal faults, Gebel Ataqa area, Gulf of Suez rift Egypt. Marine and Petroleum Geology 77, 43-53. https://doi.org/10.1016/j.marpetgeo.2016.04.012.

Maslia, M.L., Prowell, D.C., 1990. Effect of faults on fluid flow and chloride contamination in a carbonate aquifer system. J. Hydrol. 115, 1-49. https://doi.org/10.1016/0022-1694(90)90196-5.

Matsumoto, N., Shigematsu, N., 2018. In-situ permeability of fault zones estimated by hydraulic tests and continuous groundwater-pressure observations. Earth Planets Space 70, 13. https://doi.org/10.1186/s40623-017-0765-5.

Mayolle, S., Soliva, R., Caniven, Y., Wibberley, C., Ballas, G., Milesi, G., Dominguez, S., 2019. Scaling of fault damage zones in carbonate rocks. J. Struct. Geol. 124, 35-50. https://doi.org/10.1016/j.jgsl.2019.03.007.

Micarelli, L., Moretti, I., Daniel, J.M., 2003. Structural properties of rift-related normal faults: the case study of the Gulf of Corinth, Greece. Journal of Geodynamics, Active Faults: Analysis, Processes and Monitoring 36, 275-303. https://doi.org/10.1016/S0264-3707(03)00051-6.

Micarelli, L., Benedicto, A., Wibberley, C.A.J., 2006. Structural evolution and permeability of normal fault zones in highly porous carbonate rocks. J. Struct. Geol. 28, 1214-1227. https://doi.org/10.1016/j.jgsl.2006.03.036.

Mitchell, T.M., Faulkner, D.R., 2009. The nature and origin of off-fault damage surrounding strike-slip fault zones with a wide range of displacements: a field study from the Atacama fault system, northern Chile. J. Struct. Geol. 31, 802-816. https://doi.org/10.1016/j.jgsl.2009.05.002.

Mitchell, T.M., Faulkner, D.R., 2012. Towards quantifying the matrix permeability of fault damage zones in low porosity rocks. Earth Planet. Sci. Lett. 339-340, 24-31. https://doi.org/10.1016/j.epsl.2012.05.014.

Nooruddin, H.A., Hossain, M.E., Al-Yousef, H., Okasha, T., 2014. Comparison of permeability models using mercury injection capillary pressure data on carbonate rock samples. J. Pet. Sci. Eng. 121, 9-22. https://doi.org/10.1016/j petrol.2014.06.032.

Paloc, H., 1964. Caracteristiques hydrogeologiques des dolomies de la region languedocienne.

Parise, M., Closson, D., Gutiérrez, F., Stevanović, Z., 2015. Anticipating and managing engineering problems in the complex karst environment. Environ Earth Sci 74, 7823-7835. https://doi.org/10.1007/s12665-015-4647-5.

Pollock, D.W., 2016. MODPATH: A particle-tracking model for MODFLOW. https://doi.org/10.5066/F70P0X5X.

Pryet, A., Matran, P., Cousquer, Y., Roubinet, D., 2022. Particle tracking as a vulnerability assessment tool for drinking water production. Frontiers in Earth Science 10.

Rath, A., Exner, U., Tschegg, C., Grasemann, B., Laner, R., Draganis, E., 2011. Diagenetic control of deformation mechanisms in deformation bands in a carbonate grainstone. AAPG Bull. 95, 1369-1381. https://doi.org/10.1306/01031110118.

Rawling, G.C., Goodwin, L.B., Wilson, J.L., 2001. Internal architecture, permeability structure, and hydrologic significance of contrasting fault-zone types. Geol 29, 43. https://doi.org/1130/0091-7613(2001)209<0043-IAPSAH>2.0.CO;2.

Robineau, T., Tognelli, A., Goblet, P., Renard, F., Schaper, L., 2018. A double medium approach to simulate groundwater level variations in a fissured karst aquifer. J. Hydrol. 565, 861-875. https://doi.org/10.1016/j.jhydrol.2018.09.002.

Romano, V., De'Haven Hyman, J., Karra, S., Valocchi, A.J., Battaglia, M., Bigi, S., 2017. Numerical modeling of fluid flow in a fault zone: a case of study from Maژella Mountain (Italy). Energy Procedia, European Geosciences Union General Assembly 2017, EGU Division Energy, Resources & Environment (ERE) 125, 556-560. https://doi.org/10.1016/j.egypro.2017.08.191.

Rotevatn, A., Bastesen, E., 2014. Fault linkage and damage zone architecture in tight carbonate rocks in the Suez Rift (Egypt): implications for permeability along segmented normal faults. Geol. Soc. Lond. Spec. Publ. 374, 79-95. https://doi.org/10.1144/SP374.12.

Rotevatn, A., Torabi, A., Fossen, H., Braathen, A., 2008. Slipped deformation bands: a new type of cataclastic deformation bands in Western Sinai, Suez rift Egypt. Journal of Structural Geology 30, 1317-1331. https://doi.org/10.1016/j.jgsl.2008.06.010.

Sauro, F., 2014. Structural and lithological guidance on speleogenesis in quartz-sandstone: evidence of the arenisation process. Geomorphology 226, 106-123. https://doi.org/10.1016/j geomorph.2014.07.033.

Savage, H.M., Brodsky, E.E., 2011. Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones. J. Geophys. Res. Solid Earth 116. https://doi.org/10.1029/2010JB007665.

Scanlon, B.R., Mace, R.E., Barrett, M.E., Smith, B., 2003. Can we simulate regional groundwater flow in a karst system using equivalent porous media models? Case study, Barton Springs Edwards aquifer, USA. J. Hydrol. 276, 137-158. https://doi.org/10.1016/S0022-1694(03)00064-7.
Shipton, Z.K., Cowie, P.A., 2001. Damage zone and slip-surface evolution over µm to km scales in high-porosity Navajo sandstone, Utah. J. Struct. Geol. 23, 1825-1844. https://doi.org/10.1016/S0191-8141(01)00035-9.

Sibson, R.H., 1977. Fault rocks and fault mechanisms. J. Geol. Soc. London 133, 191-213. https://doi.org/10.1144/gsjgs.133.3.0191.

Singhal, B.B.S., Gupta, R.P., 2010. Applied hydrogeology of fractured rocks. Springer Science & Business Media.

Sivelle, V., Jourde, H., Bittner, D., Mazzilli, N., Tramblay, Y., 2021. Assessment of the relative impacts of climate changes and anthropogenic forcing on spring discharge of a Mediterranean karst system. J. Hydrol. 598, 126396. https://doi.org/10.1016/j.jhydrol.2021.126396.

Taillefer, A., Milesi, G., Soliva, R., Monnier, L., Delorme, P., Guillou-Frottier, L., Le Goff, E., 2021. Polyphased brittle deformation around a crustal fault: a multi-scale approach based on remote sensing and field data on the mountains surrounding the Tet hydrothermal system (Eastern Pyrénées, France). Tectonophysics 804, 228710. https://doi.org/10.1016/j.tecto.2020.228710.

Tihonov, A.N., 1963. Solution of incorrectly formulated problems and the regularization method. Soviet Math. 4, 1035-1038.

Tonkin, M.J., Doherty, J., 2005. A hybrid regularized inversion methodology for highly parameterized environmental models. Water Resour. Res. 41. https://doi.org/10.1029/2005WR003995.

Tveranger, J., Braathen, A., Skar, T., Skauge, A., 2005. Centre for integrated petroleum research - research activities with emphasis on fluid flow in fault zones. Nor. J. Geol.

van Beynen, P.E., Niedzielski, M.A., Bialkowska-Jelinska, E., Alsharif, K., Matusick, J., 2012. Comparative study of specific groundwater vulnerability of a karst aquifer in central Florida. Appl. Geogr. 32, 868-877. https://doi.org/10.1016/j.apgeog.2011.09.005.

Wang, X., Jardani, A., Jourde, H., Lonergan, L., Cosgrove, J., Gosselin, O., Massonnat, G., 2016. Characterisation of the transmissivity field of a fractured and karstic aquifer, Southern France. Adv. Water Resour. 87, 106-121. https://doi.org/10.1016/j.advwatres.2015.10.014.

White, J.T., Lavenue, M., 2023. Advances in the pilot point inverse method: Où En Sommes-Nous maintenant? Comptes Rendus. Géoscience. https://doi.org/10.5802/ergeos.161.

White, J.T., Hunt, R.J., Fienen, M.N., Doherty, J.E., 2020a. Approaches to highly parameterized inversion: PEST++ Version 5, a software suite for parameter estimation, uncertainty analysis, management optimization and sensitivity analysis (USGS Numbered Series No. 7-C26), Approaches to highly parameterized inversion: PEST++ Version 5, a software suite for parameter estimation, uncertainty analysis, management optimization and sensitivity analysis, Techniques and Methods. U.S. Geological Survey, Reston, VA. https://doi.org/10.3133/tm7C26.

White, J.T., Fienen, M.N., Doherty, J.E., 2016. A python framework for environmental model uncertainty analysis. Environ. Model. Softw. 85, 217-228. https://doi.org/10.1016/j.envsoft.2016.08.017.

White, J.T., Foster, L.K., Fienen, M.N., Knowling, M.J., Hemmings, B., Winterle, J.R., 2020. Toward reproducible environmental modeling for decision support: a worked example. Front. Earth Sci. 8. https://doi.org/10.3389/feart.2020.00050.

White, J.T., 2018. A model-independent iterative ensemble smoother for efficient history-matching and uncertainty quantification in very high dimensions. Environmental Modelling and Software 11.

Wibberley, C.A.J., Shimamoto, T., 2003. Internal structure and permeability of major strike-slip fault zones: the Median Tectonic Line in Mie Prefecture, Southwest Japan. J. Struct. Geol. 25, 59-78. https://doi.org/10.1016/S0191-8141(02)00014-7.

Wibberley, C.A.J., Yielding, G., Di Toro, G., 2008. Recent advances in the understanding of fault zone internal structure: a review. Geol. Soc. Lond. Spec. Publ. 299, 5-33. https://doi.org/10.1144/SP299.2.

Williams, A.J., Crossey, L.J., Karlstrom, K.E., Newell, D., Person, M., Woolsey, E., 2013. Hydrogeochemistry of the Middle Rio Grande aquifer system — fluid mixing and salinization of the Rio Grande due to fault inputs. Chem. Geol. 351, 281-298. https://doi.org/10.1016/j.chemgeo.2013.05.029.

Wilson, J.E., Chester, J.S., Chester, F.M., 2003. Microfracture analysis of fault growth and wear processes, Punchbowl Fault, San Andreas system, California. J. Struct. Geol. 25, 1855-1873. https://doi.org/10.1016/S0191-8141(03)00036-1.

Wu, G., Zhao, K., Qu, H., Scarselli, N., Zhang, Y., Han, J., Xu, Y., 2020. Permeability distribution and scaling in multi-stages carbonate damage zones: Insight from strike-slip fault zones in the Tarim Basin NW China. Marine and Petroleum Geology 114, 104208. https://doi.org/10.1016/j.marpetgeo.2019.104208.

Yang, J., Nie, Y., Chen, H., Wang, S., Wang, K., 2016. Hydraulic properties of karst fractures filled with soils and regolith materials: Implication for their ecohydrological functions. Geoderma 276, 93-101. https://doi.org/10.1016/j.geoderma.2016.04.024.

Zha, Y., Yeh, T.-C.-J., Illman, W.A., Tanaka, T., Bruines, P., Onoe, H., Saegusa, H., 2015. What does hydraulic tomography tell us about fractured geological media? A field study and synthetic experiments. Journal of Hydrology, Groundwater Flow and Transport in Aquifers: Insights from Modeling and Characterization at the Field Scale 531, 17-30. https://doi.org/10.1016/j.jhydrol.2015.06.013.

Zhou, H., Gómez-Hernández, J.J., Li, L., 2014. Inverse methods in hydrogeology: evolution and recent trends. Adv. Water Resour. 63, 22-37. https://doi.org/10.1016/j.advwatres.2013.10.014.

Zhuang, C., Zhou, Z., Illman, W., Guo, Q., Wang, J., 2017. Estimating hydraulic parameters of a heterogeneous aquitard using long-term multi-extensometer and groundwater level data. Hydrogeol. J. 25. https://doi.org/10.1007/s10040-017-1596-y.

Zimmerman, R., Main, I., 2004. Hydromechanical behavior of fractured rocks. Int. Geophys. Ser. 363-432.