# Analytical solution and Bayesian inference for interference pumping tests in fractal dual-porosity media

Mohamed Hayek, Anis Younes, Jabran Zouali, Noura Fajraoui, Marwan Fahs

► To cite this version:

Mohamed Hayek, Anis Younes, Jabran Zouali, Noura Fajraoui, Marwan Fahs. Analytical solution and Bayesian inference for interference pumping tests in fractal dual-porosity media. Computational Geosciences, 2018, 22 (1), pp.413-421. 10.1007/s10596-017-9701-9. hal-03527143

HAL Id: hal-03527143

https://hal.science/hal-03527143

Submitted on 15 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2 Analytical solution and Bayesian inference for interference pumping tests in fractal du-

3 al-porosity media

4

5 Mohamed Hayek<fnref n="1" />, Anis Younes<fnref n="2" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, Jabran Zouali<fnref n="2" />, Noura Fajraoui<fnref n="5" />, Marwan Fahs<fnref n="2" />

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

Submitted to Computational Geosciences

<footnote n="1">¹AF-Consult, Groundwater Protection and Waste Disposal, Täfernstrasse 26, CH-5405 Baden, Switzerland</footnote>

<footnote n="2">²LHyGES, Université de Strasbourg/EOST, CNRS, 1 rue Blessig, 67084 Strasbourg, France</footnote>

<footnote n="3">³IRD UMR LISAH, F-92761 Montpellier, France</footnote>

<footnote n="4">⁴ LMHE, Ecole Nationale d'Ingénieurs de Tunis, Tunisie</footnote>

<footnote n="5">⁵ Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland</footnote>

<footnote>* Contact person: younes@unistra.fr</footnote>
21

22 **Abstract**

23 A new analytical solution is developed for interference hydraulic pumping tests in fractured
24 porous media using the dual-porosity concept. Heterogeneous fractured reservoirs are consid-
25 ered with hydrodynamic parameters assumed to follow power-law functions in radial dis-
26 tance. The developed analytical solution is verified by comparison against a finite volume
27 numerical solution. The comparison shows that the numerical solution converges toward the
28 analytical one when the size of the time step decreases. The applicability of the fractal dual-
29 porosity model is then assessed by investigating the identifiability of the parameters from a
30 synthetic interference pumping test with a set of noisy data using Bayesian parameter infer-
31 ence. The results show that if the storage coefficient in the matrix is fixed, the rest of the pa-
32 rameters can be appropriately inferred, otherwise, the identification of the parameters is faced
33 with convergence problems because of equifinality issues.

34 **Keywords:** Fractured porous media, interference pumping test, dual-porosity, fractal media,
35 Bayesian inversion.

36
37 **1 INTRODUCTION**

Simulation of interference pumping tests in porous media allows testing different prospective exploitation scenarios at low costs. These tests consist in pumping a well at a controlled flow rate and measuring the water level response (drawdown) in one or more surrounding observation wells. Pumping tests in fractured porous media have been discussed by several authors (see [1-8], among others). Fractured media are often addressed using the dual-porosity model suggested by Barenblatt et al. [9]. The main idea of this model is to consider two overlapping continua: (i) the fractures which carry the major part of water to the pumping well and (ii) the low permeability rock matrix which feeds fluid to the fractures. At each point of the domain, water flow between the fractures and matrix are coupled via a linear exchange term proportional to the difference in hydraulic head between the two continua. As stated by De Smedt [8], the simulation of interference pumping tests is often based on numerical models which cannot insure full control of numerical accuracy [10]. Therefore, analytical solutions are highly sought-after because (i) they provide valuable insights into fundamental physics, as they are free of numerical errors and (ii), although limited to simplified situations; they allow validation of more complex numerical codes.

In the context of interference pumping tests in fractured porous media, an analytical solution has been developed for homogenous aquifers by De Smedt [8]. Furthermore, Delay et al. [7] developed a numerical solution using the dual-medium approach to describe radial convergent flow in heterogeneous fractal media where power laws in space are used for the hydrodynam- ic parameters. As a consequence, parameter values decrease with the distance between the pumped well and the observed ones. This type of scaling laws inherits from theoretical works on fractals [11] and has been used by several authors for single porous media ([12-15], among others). Hence, the main objective of this paper is to develop, for the first time, an analytical solution for the mathematical model describing interference pumping tests in heterogeneous
62 fractal dual media. Applicability of the model is then assessed by investigating identifiability
63 of the different parameters from a synthetic pumping test with a set of noisy data using Statis-
64 tical Calibration (SC) [16]. The latter means refining the prior distribution of uncertain pa-
65 rameters based on matching simulation outputs with data [17]. In this work, the SC is per-
66 formed with the DREAM(ZS) software [16] based on the Markov Chain Monte Carlo process
67 (MCMC). DREAM(ZS) computes multiple sub-chains in parallel to thoroughly explore the
68 parameter space. Taking the last 25% of individuals of the MCMC (when the chains have
69 converged) yields multiple sets of parameters that adequately fit the model onto observations.
70 These sets can then be used to estimate the updated parameter distributions and therefore the
71 best estimates of the parameters and the associate confidence intervals.

## 2 MATHEMATICAL MODEL

73 We assume a horizontal, isotropic, confined and infinite aquifer with a uniform thickness. The
74 pumping well only draws water from the fractures, has a constant rate and fully penetrates the
75 aquifer. Using the dual-porosity concept, the flow in fractured porous media can be described
76 by the following equations [18, 19]:

$$S_f \frac{\partial h_f}{\partial t} = \nabla \cdot (K_f \nabla h_f) + \alpha (h_m - h_f) \quad (1)$$

$$S_m \frac{\partial h_m}{\partial t} = \alpha (h_f - h_m) \quad (2)$$

79 where indexes *m* and *f* refer to the matrix and fracture continua, respectively, $S[L^{-1}]$ is the
80 specific storage, $h[L]$ is the hydraulic head, $K_f[LT^{-1}]$ is the hydraulic conductivity of the
81 fractures, $\alpha[L^{-1}T^{-1}]$ is the exchange rate coefficient between the fractures and the matrix.
82 Following Warren and Root [18], equation (2) assumes pseudo-steady flow in the matrix con-
tinuum (term $K_m$ is removed). Assuming radial symmetry, equations (1) and (2) can be for-
mulated as:

$$S_f(r) \frac{\partial h_f(r,t)}{\partial t} = \frac{1}{r} \frac{\partial}{\partial r} \left[ r K_f(r) \frac{\partial h_f(r,t)}{\partial r} \right] + \alpha(r) [h_m(r,t) - h_f(r,t)] \quad (3)$$

$$S_m(r) \frac{\partial h_m(r,t)}{\partial t} = \alpha(r) [h_f(r,t) - h_m(r,t)] \quad (4)$$

where $r[L]$ is the radial coordinate with the center of the pumping well as origin.

The drawdown in an observation well intercepting the fractures can be expressed as follows

$$S_f(r) \frac{\partial H_f(r,t)}{\partial t} = \frac{1}{r} \frac{\partial}{\partial r} \left( r K_f(r) \frac{\partial H_f(r,t)}{\partial r} \right) + \alpha(r) [H_m(r,t) - H_f(r,t)] \quad (5)$$

$$S_m(r) \frac{\partial H_m(r,t)}{\partial t} = \alpha(r) [H_f(r,t) - H_m(r,t)] \quad (6)$$

where $H_f = h_0 - h_f$ and $H_m = h_0 - h_m$ are respectively the drawdown in the fracture and in the rock matrix. Initially, a uniform hydraulic head ($h_0$) is assumed in both the fracture and matrix continua. Therefore, the drawdowns satisfy the following initial conditions:

$$H_f(r,0) = H_m(r,0) = 0 \quad (7)$$

Far away from the pumping well, the drawdowns are assumed to be zero

$$\lim_{r \to \infty} H_f(r,t) = \lim_{r \to \infty} H_m(r,t) = 0 \quad (8)$$

The limit where the radius approaches zero corresponds to the well withdrawal rate

$$\lim_{r \to 0} \left[ -2\pi r B K_f(r) \frac{\partial H_f(r,t)}{\partial r} \right] = Q_0, \quad t > 0 \quad (9)$$

where $B$ is the thickness of the aquifer and $Q_0$ is the pumping flow rate.

In the sequel, we limit our investigation for the semi-analytical solution of equations (5-9) to the case of a confined reservoir of a fractal structure with a single well in an unbounded me-
dium as described in [7, 12]. Note that more complex fractal approaches have been developed in the case of many wells with variable rates [20] and by considering transient flow from block to fissure and skin between the fissure and the block [21].

In this work, following Delay et al. [7], we consider dual porosity media in which the hydrodynamic parameters are power-law functions in space (decreasing functions of the radial distance $r$ between the pumping well and the observed one). Moreover, for sake of simplicity, we assume the same functions for the exchange rate coefficient $\alpha$ and for the storage coefficients $S_f$ and $S_m$. Hence, the values of the hydraulic parameters vary with respect to radial distance according to the following power-law functions:

$$ 
\begin{aligned}
K_f(r) &= K_{f0} r^{-a}; & S_f(r) &= S_{f0} r^{-b} \\
S_m(r) &= S_{m0} r^{-b}; & \alpha(r) &= \alpha_0 r^{-b}
\end{aligned}
\qquad (10) $$

where the power-law exponents $a$ and $b$ are assumed to be positive.

## 3 ANALYTICAL SOLUTION

### 3.1 Analytical solution in the Laplace domain

In this section, we follow the strategy developed by De Smedt [8] to derive the analytical solution of the system formed by equations (5)-(10). The solution is derived using the Laplace transform technique with respect to time variable $t$. The Laplace transform of a time-dependent function $y$ is

$$ \bar{y}(p) = L[y(t); t \to p] = \int_{0}^{\infty} y(t)e^{-pt}dt \qquad (11) $$

Applying the Laplace transform to (5)-(6) and using (10) leads to

$$ pS_{f0}r^{-b}\bar{H}_f(r, p) = \frac{K_{f0}}{r} \frac{\partial}{\partial r} \left[ r^{1-a} \frac{\partial \bar{H}_f(r, p)}{\partial r} \right] + \alpha_0 r^{-b} \left[ \bar{H}_m(r, p) - \bar{H}_f(r, p) \right] \qquad (12) $$

$$ pS_{m0}\bar{H}_m(r, p) = \alpha_0[\bar{H}_f(r, p) - \bar{H}_m(r, p)] \qquad (13) $$
Equation (13) can be rearranged

$$
124 \qquad \bar{H}_m(r,p) = \frac{\alpha_0 / S_{m_0}}{p + \alpha_0 / S_{m_0}} \bar{H}_f(r,p) \qquad (14)
$$

125 Substituting (14) into (12) yields

$$
126 \qquad r^2 \frac{\partial^2 \bar{H}_f}{\partial r^2} + (1-a) r \frac{\partial \bar{H}_f}{\partial r} - \frac{S_{f_0}}{K_{f_0}} \left[ p + \frac{(\alpha_0 / S_{f_0}) p}{p + \alpha_0 / S_{m_0}} \right] r^{a-b+2} \bar{H}_f = 0 \qquad (15)
$$

127 Equation (15) is a second order ordinary differential equation; its general solution is of the
128 form (cf. Polyanin and Zaitsev [22], page 159, equation (127))

$$
129 \qquad \bar{H}_f(r,p) = r^{\frac{a}{2}} \left\{ C_1 J_v \left[ \frac{2iS_{f_0}^{1/2} r^{n/2}}{nK_{f_0}^{1/2}} \sqrt{A(p)} \right] + C_2 Y_v \left[ \frac{2iS_{f_0}^{1/2} r^{n/2}}{nK_{f_0}^{1/2}} \sqrt{A(p)} \right] \right\} \quad (16)
$$

130 where $n = a - b + 2$, $v = \frac{a}{n}$, $A(p) = p + \frac{(\alpha_0 / S_{f_0}) p}{p + \alpha_0 / S_{m_0}}$, $J_v(z)$ and $Y_v(z)$ are the Bessel functions
131 of the first and second kinds, respectively. The parameters $C_1$ and $C_2$ are constant that can be
132 calculated using boundary conditions.
133

In the Laplace domain, the boundary conditions become

$$
134 \qquad \lim_{r \to \infty} \bar{H}_f(r,p) = 0 \tag{17}
$$

135 and

$$
136 \qquad \lim_{r \to 0} \left[ -2\pi B K_{f_0} r^{1-a} \frac{\partial \bar{H}_f(r,p)}{\partial r} \right] = \frac{Q_0}{p} \qquad (18)
$$

137 The solution of (15), taking into account the boundary conditions (17) and (18) is

$$
138 \qquad \bar{H}_f(r,p) = \frac{2r^{\frac{a}{2}} \bar{Q}(p) K_v \left[ \frac{2S_{f_0}^{1/2} r^{n/2}}{nK_{f_0}^{1/2}} \sqrt{p + \frac{(\alpha_0 / S_{f_0}) p}{p + \alpha_0 / S_{m_0}}} \right]}{n^{1-\nu} \pi \Gamma(1-\nu) B K_{f_0}^{1-\nu/2} S_{f_0}^{\nu/2} \left[ \sqrt{p + \frac{(\alpha_0 / S_{f_0}) p}{p + \alpha_0 / S_{m_0}}} \right]^\nu} \quad (19)
$$
where $K_v$ is the modified Bessel function of the second kind of order $v$. More details about the derivation of the solution (19) are given in the Appendix.

141

## 3.2 Analytical solution in the time domain

The solution in the time domain $H_f(r,t)$ is sought using the convolution theorem of the Laplace transform [23]

$$L\left[\int_{0}^{t} y(\tau, t-\tau) d\tau; t \rightarrow p\right] = L\{L[y(t_1, t_2); t_2 \rightarrow p]; t_1 \rightarrow p\} \quad (20)$$

Hence, equation (19) is written as

$$\bar{H}_f(r, p) = \bar{H}_f(r; p_1, p_2) = \frac{2Q_0 r^{\frac{a}{2}} K_v \left[ \frac{2S_{f_0}^{1/2} r^{n/2}}{nK_{f_0}^{1/2}} \sqrt{p_1 + \frac{(\alpha_0 / S_{f_0}) p_2}{p_2 + \alpha_0 / S_{m_0}}} \right]}{n^{1-v} \pi \Gamma(1-v) B K_{f_0}^{1-v/2} S_{f_0}^{v/2} p_2 \left[ \sqrt{p_1 + \frac{(\alpha_0 / S_{f_0}) p_2}{p_2 + \alpha_0 / S_{m_0}}} \right]^v} \quad (21)$$

where $p_1 = p_2 = p$ but different indexes are used to indicate that the inverse transformation is performed in two steps.

The different indexes used in (21) indicate that the inverse Laplace transform is developed in two steps. In a first step, the inverse Laplace transform with respect to the variable $p_1$ is expressed as follows

$$L^{-1}\left[\bar{H}_f(r; p_1, p_2); p_1 \rightarrow t_1\right] = \frac{n^\nu Q_0 e^{-\frac{S_{f_0} r^n}{2n^2 K_{f_0} t_1}} \cdot W_{\frac{1-\nu}{2}, \frac{\nu}{2}}\left(\frac{S_{f_0} r^n}{n^2 K_{f_0} t_1}\right)}{\pi \Gamma(1-\nu) B K_{f_0}^{\frac{1-\nu}{2}} S_{f_0}^{\frac{1+\nu}{2}} r^{\frac{2-b}{2}}} \cdot e^{-\frac{(\alpha_0 / S_{f_0}) t_1 p_2}{p_2 + \alpha_0 / S_{m_0}}} \quad (22)$$

where $t_1$ is the time variable obtained by the inverse Laplace transform with respect to the variable $p_1$ and $W_{\kappa, \mu}(z)$ is the Whittaker function (Abramowitz and Stegun [24] p. 505).
156 In the second step, the inverse Laplace transform with respect to $p_2$ is obtained using the inverse transform of the J-function [25] based on the formula

157
158

$$
L^{-1} \left[ \frac{1}{p} e^{-\frac{\alpha p}{p+\beta}} ; p \to t \right] = J(\alpha, \beta t)
$$

which yields

$$
159
$$
$$
\bar{H}_f (r; p_1, p_2) = \frac{n^\nu Q_0 e^{-\frac{S_{f_0} r^n}{2n^2 K_{f_0} t_1}} \cdot t_1^{\frac{\nu-1}{2}} \cdot W_{\frac{1-\nu}{2}, \frac{\nu}{2}} \left( \frac{S_{f_0} r^n}{n^2 K_{f_0} t_1} \right) \cdot J \left( \frac{\alpha_0 t_1}{S_{f_0}}, \frac{\alpha_0 t_2}{S_{m_0}} \right)}{2\pi\Gamma(1-\nu) B K_{f_0}^{\frac{1-\nu}{2}} S_{f_0}^{\frac{1+\nu}{2}} r^{\frac{2-b}{2}}} \quad (23)
$$

160 The Bessel integral J-function is defined as follows [8, 26]

$$
161
$$
$$
J(a,b) = 1 - e^{-b} \int_{0}^{a} e^{-\lambda} I_{0}(2\sqrt{b\lambda}) d\lambda = e^{-a-b} \sum_{n=0}^{\infty} \frac{b^n}{n!} \sum_{m=0}^{n} \frac{a^m}{m!} \quad (24)
$$

162 In view of (20) and (23), the analytical expression of $H_f(r,t)$ is

$$
163
$$
$$
H_f(r,t) = \frac{n^\nu Q_0}{2\pi\Gamma(1-\nu)B K_{f_0}^{\frac{1-\nu}{2}} S_{f_0}^{\frac{1-\nu}{2}} r^{\frac{1-\nu}{2}}} \int_0^t \tau^{\frac{1-\nu}{2}-1} e^{-\frac{S_{f_0} r^n}{2n^2 K_{f_0} \tau}}} \cdot W_{\frac{1-\nu}{2}, \frac{\nu}{2}} \left( \frac{S_{f_0} r^n}{n^2 K_{f_0} \tau} \right) \cdot J \left( \frac{\alpha_0 \tau}{S_{f_0}}, \frac{\alpha_0 (t-\tau)}{S_{m_0}} \right) d\tau \quad (25)
$$

165 A similar expression can be obtained for $H_m(r,t)$ using the convolution theorem of the Laplace transform for equation (14) and using equation (19).
166

167 In the case of homogeneous media ($a=b=0$), we have $W_{\frac{1}{2},0}(z) = z^2 e^{-\frac{z}{2}}$ and (25) reduces to
168 the De Smedt [8] solution:

$$
169
$$
$$
H_f(r,t) = \frac{Q_0}{4\pi B K_{f_0}} \int_0^t \tau^{-1} e^{-\frac{S_{f_0} r^2}{4K_{f_0} \tau}} \cdot J\left(\frac{\alpha_0 \tau}{S_{f_0}}, \frac{\alpha_0 (t-\tau)}{S_{m_0}}\right) d\tau \quad (26)
$$

170

171

172
173 **4 COMPARISON TO NUMERICAL SOLUTIONS**

174 We numerically solve two interference pumping tests dealing with homogeneous and fractal
175 dual media. For each test, the numerical solution is compared against the corresponding ana-
176 lytical solution using (25). The numerical solution of the system (5)-(10) is developed using
177 the standard Finite Volume (FV) method. The domain of size L=10,000m is discretized us-
178 ing a uniform spatial discretization of size Δr = 5 m. A similar solution is obtained using a
179 finer spatial discretization of size Δr = 1 m which shows the weak sensitivity of the numerical
180 solution to the spatial discretization. The duration of the simulation is 2×10⁶ s. Three time
181 step sizes are investigated with Δt = 10⁵ s, Δt = 5×10⁴ s and Δt = 10⁴ s. The pumping flow rate
182 is Q₀ = 5 m³/h. The initial head in the aquifer is h₀ = 100 m and the aquifer thickness is
183 B = 100 m. The specific storage of the fractures and rock matrix are respectively Sf₀ = 5×10⁻⁶
184 m⁻¹ and Sm₀ = 10⁻⁴ m⁻¹. The hydraulic conductivity of the fractures is Kf₀ = 10⁴ m/s and the
185 rate of exchange between the matrix and the fractures is α₀ = 5×10⁻¹¹ m⁻¹s⁻¹.

186 Figure 1 depicts the drawdown in both the fracture and matrix at 10m and 100m. The draw-
187 down in the fissure starts at early times because storage of the fissures is immediately effec-
188 tive at the start of the pumping test, while storage of the matrix needs much more time to be-
189 come noticeable. A flattening of the drawdown curve in the fissure can be observed at inter-
190 immediate pumping time which is typical for double-porosity behavior because of the transition
191 from fracture to matrix water storage release.

192 The numerical FV solution is compared to the analytical solution (25) in the case of homoge-
193 neous (a = b = 0) and fractal (a = 0.6 and b = 0.4) dual media. The comparison is performed
194 at two observation wells located respectively at distances r₁=10m and r₂=100m from the
195 pumping well. Significant discrepancies can be observed between analytical and numerical
solutions when the latter is evaluated using a large time step size (Figure 2). The drawdown in the case of fractal dual media is larger than in homogeneous dual media because of the reduction of the hydrodynamic parameters with the distance $r$ from the pumping well. In fact, in the case of fractal dual media, the permeability of the fractures is proportional to $r^{-0.6}$ whereas the exchange rate and the specific storage of both the fractures and the matrix are proportional to $r^{-0.4}$. Thereby, the discrepancy between analytical and numerical solutions (Figure 2) is more pronounced with dual fractal media than with homogeneous dual media. For both homogeneous and fractal dual media, the numerical solution converges toward the analytical solution (25) when the size of the time step decreases (Figure 2). Note that the reduction of the permeability of the fractures with respect to the distance $r$ is also responsible, via the Whittaker function in equation (25), of the delay for the stabilisation of the solution in the fractal case (Figure 2 b) compared to the homogeneous case (Figure 2a).

## 5 BAYESIAN INVERSION OF A PUMPING TEST IN FRACTAL DUAL MEDIA

In practical applications, most of the hydraulic parameters cannot be measured directly and should be estimated using an inversion procedure. In the case of fractal dual media, the inversion can be challenging because of model equifinality, implying that different parameter combinations can fit the model responses to data. To assess the applicability of the fractal dual porosity model, the identifiability of the parameters of the analytical solution (25) is assessed using a synthetic interference pumping test with a set of noisy data. Parameter identification is performed in a Bayesian framework where the prior knowledge about the parameters and the observed data are merged to define the joint posterior probability distribution function (pdf) of the parameters. In this work, the pdf is performed using the DREAM(ZS) software [16] based on the Markov Chain Monte Carlo (MCMC) sampler. DREAM(ZS) generates random sequences of parameter sets that asymptotically converge toward the target joint posterior dis-
tribution [27]. Thus, if the number of runs is sufficiently high, the generated samples allows exploring the entire parameter space of the posterior distribution of the parameters and provides the pairwise parameter correlations and the uncertainty of model predictions.

The Bayes theorem states that the probability density function of the model parameters conditioned onto data can be expressed as:

$$p(\xi | \mathbf{y}_{mes}) \propto p(\mathbf{y}_{mes} | \xi) p(\xi), \qquad (27)$$

where $p(\xi | \mathbf{y}_{mes})$ is the likelihood function measuring how well the model fits the observations $\mathbf{y}_{mes}$, and $p(\xi)$ is the prior information about the parameter before the observations are made. In this work, a Gaussian distribution defines the likelihood function because the observations are simulated and corrupted with Gaussian errors. Hence, the parameter posterior distribution is expressed as:

$$p(\xi | \mathbf{y}_{mes}) \propto \sigma_h^{-N_h} \exp \left( -\frac{SS_h(\xi)}{2\sigma_h^2} \right), \qquad (28)$$

where $SS_h(\xi)$, are the sums of the squared differences between the observed and modeled drawdowns. For instance, $SS_h(\xi) = \sum_{k=1}^{N_h} (h_{mes}^{(k)} - h_{mod}^{(k)}(\xi))^2$, which includes the observed $h_{mes}^{(k)}$ and predicted $h_{mod}^{(k)}$ drawdowns at time $t_k$ for the number of observations $N_h$.

A synthetic pumping test is generated by running the analytical solution using the reference parameter values given in Table 1 for $K_{f0}$, $S_{f0}$, $\alpha_0$, $a$ and $b$. As in Fahs et al. [28], the storage coefficient in the matrix $S_{m0}$ is fixed in this study to be $10^{-4}$ m$^{-1}$.

Delay et al [7] showed that results of calibration improved when several drawdown curves at several locations are inverted at once. In the sequel, the analytical solution is run for a long period of $2 \times 10^7$ s to calculate the drawdown at four observation wells located respectively at $r_1 = 1$ m, $r_2 = 10$ m, $r_3 = 100$ m and $r_4 = 500$ m from the pumping well. Four datasets of 2000
242 values each (corresponding to the drawdown at the observation wells measured each $10^4$ s)
243 are used as target responses for model inversion.

244 Gaussian noise with a standard deviation of 0.01m and a mean of zero was added to the target
245 responses. We assume that all hydraulic parameters have uniform prior distributions over the
246 ranges given in Table 1. As recommended in Vrugt et al. [29], we consider that the posterior
247 distribution is stationary if the Gelman and Ruban [30] criterion is less than 1.2 and that the
248 chains are not autocorrelated.

249 The MCMC method was terminated after 15000 model runs. The convergence was reached at
250 around 10000 model runs. Figure 3 shows the results of the identification based on the last
251 25% of individuals of the MCMC sampler (when the chains have converged). The "on-
252 diagonal" plots in these figures display the inferred parameter distributions, whereas the "off-
253 diagonal" plots represent the pairwise correlations in the MCMC sample. If the drawings are
254 independent, non-sloping scatterplots should be observed. However, if a good value of a giv-
255 en parameter is conditioned by the value of another parameter, then their pairwise scatterplot
256 should show a narrow sloping stripe. The results of Figure 3 show that the five parameters
257 $K_{f0}$, $S_{f0}$, $\alpha_0$, $a$ and $b$ are appropriately estimated; they have almost symmetric bell-shaped
258 posterior distributions and have strongly narrowed their prior intervals. The scatter plots in
259 this figure show a moderate correlation between the parameters $a$ and $b$ (a correlation factor
260 of 0.61), a high correlation between the parameters $S_{f0}$ and $b$ (a correlation factor of 0.79)
261 and a very high correlation between the parameters $K_{f0}$ and $a$ (a correlation factor of 0.98).
262 The mean and 95% confidence intervals of the samples that adequately fit the model onto
263 observations are reported in Table 2. From this table, we can see that the parameters $K_{f0}$, $a$
264 and $b$ are well identified with respectively a mean value of 0.99×10⁻⁵ ms⁻¹, 0.299 and 0.2.
265 These values are almost identical to the reference values of Table 1. The uncertainty on these
parameters is small since the 95% confidence intervals are respectively [0.9-1.1]×10⁻⁵ m⁻¹,
[0.28-0.32] and [0.17-0.23]. The parameters S<sub>f₀</sub> and α₀ are less well identified. Indeed, alt-
hough their mean values (respectively 5.03×10⁻⁵ m⁻¹ and 3.06×10⁻¹¹ m⁻¹s⁻¹) are very close to
the reference values, their posterior uncertainty remains quite large since their confidence
intervals are respectively [3.2-7.6]×10⁻⁵ m⁻¹ and [2.1-4.5]×10⁻¹¹ m⁻¹s⁻¹.

Finally note that the identification procedure does not converge if the storage coefficient in
the matrix S<sub>m₀</sub> is involved in the inversion procedure (not fixed). In that case, the identifica-
tion of the six parameters K<sub>f₀</sub>, S<sub>f₀</sub>, S<sub>m₀</sub>, α₀, a and b was not possible. Indeed, in this case,
the MCMC algorithm was faced to convergence problem because of **equifinality** issues.

276 **6 CONCLUSION**

A new analytical solution for interference hydraulic pumping tests in fractal fractured porous
media is developed using the dual-porosity concept and power laws in space for all hydrody-
namic parameters. The developed analytical solution is compared to a numerical one obtained
with the FV method using different time step sizes. Discrepancies are observed between ana-
lytical and numerical solutions when the latter is used with large time steps. These discrepan-
cies are more pronounced with dual fractal media than with homogeneous dual media. When
the size of the time step decreases, the numerical solution converges toward the analytical one
for both homogeneous and fractal dual media.

The applicability of the fractal dual-porosity model was then assessed by investigating identi-
fiability of the hydraulic parameters from a synthetic interference pumping test with a set of
noisy data using the MCMC sampler. The results show that if the storage coefficient in the
matrix is fixed, the rest of the parameters can be appropriately identified, whereas conver-
gence problems can be encountered if all the hydraulic parameters are to be estimated.

290

291 Acknowledgements

292 This work was partially supported by the GdR MoMaS (PACEN/CNRS, ANDRA, BRGM,
293 CEA, EDF, IRSN) France and by the French National Research Agency (ANR) through the
294 program AAP Blanc - SIMI 6 project RESAIN (n°ANR-12-BS06-0010-02).

295
## Appendix: Derivation of the Analytical Solution in the Laplace Domain

Using the boundary condition at infinity (17) and the following asymptotic expansions for large arguments ($|z| \to \infty$) of $J_v$ and $Y_v$ [22], $J_v(z) = \sqrt{2/\pi z} \{P(\nu, z)\cos\chi - Q(\nu, z)\sin\chi\}$ and $Y_v(z) = \sqrt{2/\pi z} \{P(\nu, z)\cos\chi + Q(\nu, z)\sin\chi\}$ for $|\arg z| < \pi$ with $\chi = z - \left(\frac{1}{2}\nu + \frac{1}{4}\right)\pi$, $P(\nu, z) = \sum_{k=0}^{\infty} (-1)^k \frac{(\nu, 2k)}{(2z)^{2k}}$ and $Q(\nu, z) = \sum_{k=0}^{\infty} (-1)^k \frac{(\nu, 2k+1)}{(2z)^{2k+1}}$, equation (16) gives $C_1 + iC_2 = 0$. Therefore, this latter can be written as

$$ \bar{H}_f(r, p) = C_1 r^2 H_v^{(1)} \left[ \frac{2i S_{f_0}^{1/2} r^{n/2}}{n K_{f_0}^{1/2}} \sqrt{A(p)} \right] \quad (A.1) $$

where $H_v^{(1)}$ is the Hankel function of the first kind defined by $H_v^{(1)}(z) = J_v(z) + iY_v(z)$ for any complex number $z$.

The boundary condition at the well (18) is then used to find $C_1$. Note that for $\nu > 0$ the Hankel function $H_v^{(1)} \sim -\frac{i}{\pi}\Gamma(\nu)\left(\frac{1}{2}z\right)^{-\nu}$ as $z$ approaches zero [22], where $\Gamma$ is the Gamma function. Using this expansion we can write

$$ r^{1-a} \frac{\partial H_f(r,t)}{\partial r} \sim \frac{i^{1-\nu} S_{f_0}^{\nu/2} C_1}{2\pi K_{f_0}^{\nu/2}} \left\{ -r^{-a} n^{\nu} \left[ \sqrt{A(p)} \right]^{-\nu} \{a\Gamma(\nu) - n\Gamma(1+\nu)\} + n^{1-\nu} \left[ \sqrt{A(p)} \right]^{\nu} \Gamma(1-\nu) \right\} \quad (A.2) $$
312 as $r$ approaches zero. Using the property of the Gamma function $\Gamma(1+\nu)=\nu\Gamma(\nu)$ for $\nu > 0$,

313 the first term of right hand side of equation (A.2) vanishes since $n\nu = a$. Substituting (A.2)

314 into (18) we get $C_1 = - \frac{\bar{Q}(p)}{i^{1-\nu} n^{1-\nu} \Gamma(1-\nu) B K_{f_0}^{1-\nu/2} S_{f_0}^{\nu/2} [\sqrt{A(p)}]^\nu}$.

315 Finally, inserting the value of $C_1$ into (A.1) and using the relationship

316 $K_v(z) = \frac{1}{2} \pi i e^{\frac{1}{2} \nu \pi i} H_v^{(1)}(ze^{\frac{1}{2} \pi i})$ for $-\pi < \arg z \le \frac{1}{2} \pi$, where $K_v$ is the modified Bessel func-

317 tion of the second kind of order $\nu$, we get the solution in the Laplace domain as defined by

318 (19).

319
321 REFERENCES

322 [1] Hsieh, P.A.: A brief survey of hydraulic tests in fractured rocks. In Faybishenko B, 323 Witherspoon PA, Benson SM (eds.) Dynamics of fluids in fractured rocks, American 324 Geophysical Union Monograph 122:59-66 (2000).

325 [2] Kruseman, G.P., de Ridder, N.A.: Analysis and evaluation of pumping test data. Inter- 326 national Institute for Land Reclamation and Improvement Publication, Wageningen 327 (1991).

328 [3] Batu, V.: Aquifer Hydraulics: a comprehensive guide to hydrogeologic data analysis. 329 John Wiley & Sons Inc. New York (1998).

330 [4] Cheng, A.H.D.: Multilayered aquifer systems fundamentals and applications. Marcel 331 Dekker, New York/Basel (2000).

332 [5] Walton, W.C.: Aquifer test modelling. CRC press, Taylor & Francis Group, Boca Raton 333 (2006).

334 [6] Nielsen, K.A.: Fractured aquifers formation evaluation by Well Testing. Trafford Pub- 335 lishing, Canada (2007).

336 [7] Delay, F., Kaczmaryk, A., Ackerer, P.: Inversion of interference hydraulic pumping 337 tests in both homogeneous and fractal dual media. Adv. in Water Resour. 30(3):314- 338 334. (2007) doi:10.1016/j.advwatres.2006.06.008

339 [8] De Smedt, F.: Analytical solution Analytical Solution for Constant-Rate Pumping Test 340 in Fissured Porous Media with Double-Porosity Behaviour. Transp Porous Med. 341 88:479-489 (2011). doi:10.1007/s11242-011-9750-9.

342 [9] Barenblatt, G.E., Zheltov, I.P., Kochina, I.N.: Basis concepts in the theory of seepage of 343 homogeneous liquids in fissured rocks. J. Appl. Math. Mech. Engl. Transl. 24(5):1286- 344 1303 (1960). doi:10.1016/0021-8928(60)90107-6.
[10] Zhan, H., Wen, Z., Gao, G.: An analytical solution of two-dimensional reactive solute transport in an aquifer aquitard system. Water Resour. Res. 45, W10501. (2009). doi:10.1029/2008WR007479

[11] O'Shaughnessy, B., Procaccia, I.: Diffusion in fractals. Phys Rev A. 32:3073 (1985). doi:10.1103/PhysRevA.32.3073

[12] Acuna, J.A., Yortsos, Y.C.: Application of fractal geometry to the study of networks of fractures and their pressure transient. Water Resour Res; 31(3):527-40 (1995). doi: 10.1029/94WR02260

[13] Chang, J., Yortsos, Y.C.: Pressure-transient analysis of fractal reservoirs. Soc Petrol Eng J 18710:631-43 (1990) doi:10.2118/18170-PA.

[14] Delay, F., Porel, G., Bernard, S.: Analytical 2D model to invert hydraulic pumping tests in fractured rocks with fractal behavior. Geophys Res Lett. 31(16) (2004). doi:10.1029/2004GL020500.

[15] Le Borgne, T., Bour, O., de Dreuzy, J.R., Davy, P., Touchard, F.: Equivalent mean flow models for fractured aquifers: insights for a pumping test scaling interpretation. Water Resour Res;40(3) (2004). doi:10.1029/2003WR002436.

[16] Laloy, E., Vrugt, J.A.: High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing, Water Resour. Res., 48 W01526. (2012). doi:10.1029/2011WR010608.

[17] Campbell, K.: Statistical calibration of computer simulations, Reliability Engineering and System Safety 91 1358-1363 (2006). doi:10.1016/j.ress.2005.11.032

[18] Warren, J.E., Root, P.J.: The behavior of naturally fractured reservoirs. Soc Petrol Eng 3:245-55. (1963). doi:10.2118/426-PA.

[19] Moench, A.F.: Double-porosity models for a fissured groundwater reservoir with fracture skin. Water Resour. Res. 20(7):831-846. (1984). doi:10.1029/WR020i007p00831.

[20] Chang, J., Yortsos, Y. C. A Note On Pressure-Transient Analysis Of Fractal Reservoirs. Society of Petroleum Engineers. (1993). doi:10.2118/25296-PA
372 [21] Hamm, S.Y., Bidaux, P. Dual-porosity fractal models for transient flow analysis in fis-
sured rocks: Water Resources Research, v. 32, no. 9, p. 2733-2745 (1996).
373

374 [22] Polyanin, A.D., Zaitsev, V.F. Handbook of Exact Solutions for Ordinary Differential
Equations, 2nd ed., 787 pp., Chapman and Hall/CRC, Boca Raton, Fla (2003).
375

376 [23] Sneddon, I.H.: The Use of Integral Transforms, McGraw-Hill, New York (1972).

377 [24] Abramowitz, M., Stegun, I.A. Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables, Government Printing Off., Washington, D. C. (1972).
378
379

380 [25] Goldstein, F.R.S.: On the mathematics of exchange processes in fixed columns I: Math-
ematical solutions and asymptotic expansions. Proc. R. Soc. London 219: 151-185.
381
382
(1953).

383 [26] De Smedt, F., Wierenga, P.J.: A generalized solution for solute flow in soils with mobile
384
385 and immobile water. Water Resour. Res. 15(5): 1137-1141. (1979). doi:
10.1029/WR015i005p01137.

386 [27] Gelman, A., Carlin, J.B., Stren, H.S., Rubin, D.B. Bayesian data analysis (1997), Chap-
387
mann and Hall, London.

388 [28] Fahs, H., Hayek, M., Fahs, M., Younes, A.: An efficient numerical model for hydrody-
namic parameterization in 2D fractured dual-porosity media. Adv. Water Resour. 63:
389
390
179-193. (2013). doi: 10.1016/j.advwatres.2013.11.008

391 [29] Vrugt, J.A., Gupta, H.V., Bouten, W., Sorooshian, S. A shuffled complex evolution Me-
tropolis algorithm for optimization and uncertainty assessment for hydrologic model pa-
392
393
parameters. Water Resour. Res. 39(8):1201 (2003), doi:10.1029/2002WR001642.

394 [30] Gelman, A., Rubin, D.B. Inference from iterative simulation using multiple sequences.
395
Stat. Sci. 7:457-472 (1992). doi: 10.1214/ss/1177011136

396
397

398

399 **List of Tables**

400 Table 1: Reference values, lower and upper bounds of the parameters for a synthetic pumping test case in fractal dual media.

402 Table 2: Mean values and 95% confidence intervals for the parameters estimated with the MCMC sampler.

404

405

406 **List of Figures**

407 Figure 1: Drawdown versus time curves for a pumping test in homogeneous and fractal dual media.

409 Figure 2: Comparison between analytical and numerical solutions in the case of (i) homogeneous and (ii) fractal dual media.

411 Figure 3: MCMC solutions of the interference pumping test in fractal dual media. The diagonal plots represent the inferred posterior probability distribution of model parameters. The off-diagonal scatterplots represent the value of a parameter α versus that of a parameter β in the best solutions inferred by MCMC.

415

416

417
418

419

420

421

<caption>422 Table 1</caption>

423

424

<table><thead><tr><td>Parameters</td><td>Lower bounds</td><td>Upper bounds</td><td>Reference values</td></tr></thead><tbody><tr><td>K<sub>f0</sub> [ms<sup>-1</sup>]</td><td>5×10<sup>-7</sup></td><td>5×10<sup>-4</sup></td><td>1×10<sup>-5</sup></td></tr><tr><td>S<sub>f0</sub> [m<sup>-1</sup>]</td><td>10<sup>-6</sup></td><td>10<sup>-3</sup></td><td>5×10<sup>-5</sup></td></tr><tr><td>α<sub>0</sub> [m<sup>-1</sup>s<sup>-1</sup>]</td><td>10<sup>-12</sup></td><td>10<sup>-9</sup></td><td>3×10<sup>-11</sup></td></tr><tr><td>a [-]</td><td>0.1</td><td>0.5</td><td>0.3</td></tr><tr><td>b [-]</td><td>0.1</td><td>0.5</td><td>0.2</td></tr></tbody></table>

425
426

427

428

<caption>Table 2</caption>

430

<table><thead><tr><td>Parameters</td><td>Mean</td><td>95% Conf. Int.</td></tr></thead><tbody><tr><td>K<sub>f0</sub> [ms<sup>-1</sup>]</td><td>0.99×10<sup>-5</sup></td><td>[0.9-1.1]×10<sup>-5</sup></td></tr><tr><td>S<sub>f0</sub> [m<sup>-1</sup>]</td><td>5.03×10<sup>-5</sup></td><td>[3.2-7.6]×10<sup>-5</sup></td></tr><tr><td>α<sub>0</sub> [m<sup>-1</sup>s<sup>-1</sup>]</td><td>3.06×10<sup>-11</sup></td><td>[2.1-4.5]×10<sup>-11</sup></td></tr><tr><td>a [-]</td><td>0.299</td><td>[0.28-0.32]</td></tr><tr><td>b [-]</td><td>0.2</td><td>[0.17-0.23]</td></tr></tbody></table>

431

432

433

434

435

436

437

438

439

440

441
442 **Figure 1**

443

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

444

445
446

<figure><img src="image_3.png" /><figcaption>447 Figure 2</figcaption></figure>

448

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

449

450
<figure><img src="image_5.png" /><figcaption>Figure 3</figcaption></figure>