1]Sandia National Laboratories, Nuclear Waste Disposal Research & Analysis, PO
Box 5800, Mail Stop 0747, Albuquerque NM, 87185, USA

# Generalized Solution for Double-Porosity Flow through a Graded Excavation
Damaged Zone

Kristopher L. Kuhlman<EMAIL_ADDRESS>[

###### Abstract

Prediction of flow to boreholes or excavations in fractured low-permeability
rocks is important for resource extraction and disposal or sequestration
activities. Analytical solutions for fluid pressure and flowrate, when
available, are powerful, insightful, and efficient tools enabling parameter
estimation and uncertainty quantification. A flexible porous media flow
solution for arbitrary physical dimension is derived and extended to double
porosity for converging radial flow when permeability and porosity decrease
radially as a power law away from a borehole or opening. This distribution can
arise from damage accumulation due to stress relief associated with drilling
or mining. The single-porosity graded conductivity solution was initially
found for heat conduction, the arbitrary dimension flow solution comes from
hydrology, and the solution with both arbitrary dimension and graded
permeability distribution appeared in reservoir engineering. These existing
solutions are here combined and extended to two implementations of the double-
porosity conceptual model, for both a simpler thin-film mass transfer and more
physically realistic diffusion between fracture and matrix. This work presents
a new specified-flowrate solution with wellbore storage for the simpler
double-porosity model, and a new more physically realistic solution for any
wellbore boundary condition. A new closed-form expression is derived for the
matrix diffusion solution (applicable to both homogeneous and graded
problems), improving on previous infinite series expressions.

###### keywords:

double porosity, excavation damaged zone, porous media flow, fractured rock,
salt, analytical solution

##  1 Introduction

Fluid flow through damage-induced fracture networks in otherwise low-
permeability crystalline rocks (e.g., granite, argillite or halite) is of
interest to geothermal energy production [73], radioactive waste disposal
[75], hydrogen storage [1], and compressed air energy storage [43]. Rock
damage around an excavation (i.e., the Excavation Damaged Zone, EDZ; Davies
and Bernier [23]) increases the connected porosity, and leads to increased
permeability. Fractured rock often has higher porosity and permeability than
intact rock. Damage near a borehole or excavation will decrease the relative
contribution from flow in the lower-permeability far-field, and will confound
the estimation of hydrologic properties using approaches that assume uniform
homogeneous distributions of permeability and porosity. There is a need for a
flexible analytical solution for flow to a borehole or excavation in the
presence of damage, that includes wellbore storage, double-porosity flow, and
variable flow dimension. This is most evident in a mechanically weak, low-
permeability medium like salt, but should also apply to other low-permeability
fractured rocks like granite or shale.

In salt, the far-field (i.e., undamaged) permeability is unmeasurably low [8]
due to salt’s tendency to creep shut any unsupported openings. The
permeability around a borehole in salt is derived from accumulated damage due
to stress redistribution around the excavation itself [76, 71, 20, 39, 48].

Stormont et al [71] presented brine and gas permeability data measured in salt
for packer-isolated intervals of small boreholes before and after a central
1-meter diameter borehole was drilled (i.e., a mine-by experiment). Figure 1
shows these data support the conceptual model of permeability and porosity
decaying away from an excavation. Cosenza [20] proposed the power-law model
for permeability and porosity plotted in the figure. These data show porosity
and permeability decrease with distance from the central excavation. Two lines
are shown with to the data; one is a monomial power-law, the other includes an
additive background term. The two curves differ primarily away from the
excavation ($r/r_{w}\geq 3$), where larger uncertainties in estimated porosity
and permeability exist, for three reasons. First, the access drift EDZ (test
conducted in the floor of a 5-m wide room) is superimposed on the 1-m borehole
EDZ. Second, the small-diameter (2.5-cm) measurement boreholes themselves each
have a small EDZ overprinted on the 1-m borehole EDZ. Lastly, the apparent
background permeability may represent the measurement limit of the packer
system used (i.e., compliance of the packer inflation elements and working
fluid). Especially in salt, the undisturbed background permeability is near
zero, and is difficult to measure consistently in the field [8]. The power-law
distribution of permeability matches the more certain near-field permeability
distribution, and is conceptually more elegant than a finite domain or a flow
domain with piece-wise heterogeneous properties (i.e., a higher-permeability
EDZ adjacent to lower-permeability intact rock). Other investigations have
also shown porosity and permeability decaying away with distance from an
excavation in crystalline rocks [68, 18, 33] and sedimentary rocks [65, 64].

Figure 1: Permeability and porosity observations around a 1-m borehole (radial
distance scaled by excavation radius) in salt from small-scale mine-by
experiment (data from Stormont et al [71])

Salt permeability has been related to both the confining and shear stresses
[66, 53, 70, 3]. Confining stresses reduce fracture aperture and bulk
permeability, while shear stresses are associated with increased bulk
permeability. Aydan et al [4] present solutions for radial and tangential
plane stress and strain (i.e., dilatation or a change in porosity) around a
circular excavation. Strain is proportional to $r_{D}^{-2}$ or $r_{D}^{-3}$
(where $r_{D}$ is radial distance into the formation scaled by the excavation
size), depending on whether the region is experiencing elastic (exponent $2$)
or plastic (exponent $\approx 3$) deformation. These relationships illustrate
a possible behavior of rock in the EDZ. The true extent of the EDZ depends on
drilling or excavation method, borehole or tunnel geometry, state of stress,
and rock mechanical properties [40]. Softer or weaker sedimentary rocks like
argillite or halite typically have a larger EDZ than stiffer or stronger rocks
like granite.

There are several well-known empirical power-law relationships between
porosity and permeability in fractured or granular media [e.g., 44, 14] and
many studies have discussed their applicability [22, 50]. Permeability in
fractured rocks is more sensitive to small changes in porosity than granular
rocks (i.e., fractured rocks have higher pore compressibility resulting in
larger exponents in porosity-permeability relationships).

Based on evidence from these observations, graded dimensionless porosity is
assumed to follow

$n({r})=n_{0}\left(\frac{{r}}{r_{w}}\right)^{-\eta},$ (1)

where $r_{w}$ is the borehole or excavation radius [$\mathrm{m}$],
$n_{0}=n(r_{w})$ is maximum porosity at the borehole wall, and $\eta$ is a
dimensionless exponent (see Table 1 for a list of physical variables and
notation). Using the same form, the graded permeability can be represented
with the form

$k({r})=k_{0}\left(\frac{{r}}{r_{w}}\right)^{-\kappa},$ (2)

where $k_{0}=k(r_{w})$ is the maximum permeability [$\mathrm{m^{2}}$] at the
borehole wall and $\kappa$ is another dimensionless exponent. Based on lab
measurements on fractured granite, the empirical relationship $\kappa\approx
3\eta$ has been proposed [45, 22]. The Stormont et al [71] salt data (Figure
1) support $\eta=4.5$ and $\kappa=17$, which shows a somewhat faster-decaying
permeability ($\kappa=3.8\eta$) than seen in granitic rocks.

The power-law permeability and porosity distribution conceptual model
presented here is an alternative to flow models using wellbore skin [72, 63],
finite domain [32, 55], or low-permeability non-Darcy flow with a threshold
gradient [56, 57]. These three conceptualizations all lead to reduced
contributions of flow from the far field, but only borehole skin can account
for observed distributions of higher porosity or permeability near the
excavation, which are important when analyzing pressure or flowrate data at
early time. The contribution from lower permeability in the far field are more
important at late time. Finite domains and skin can have analytical flow
solutions, but low-permeability non-Darcy flow does not typically lend itself
to analytical solutions.

Barker [7] developed a generalized solution for converging flow to a borehole
with variable non-integer dimension, $D$. This conceptualization has been used
to characterize flow in fractured systems, where lower-dimension (i.e., $D<3$)
results associated with discrete fractures are more common than higher
dimension results [9, 54, 11, 31]. Doe [30] extended the solution of Barker
[7] to the conceptualization where permeability varies with radial distance,
through analogy with the heat conduction literature [15].

A single-porosity flow solution is derived here with power-law variable
properties, like the approach of Doe [30] (who did not present a derivation).
The single-porosity solution is then readily extended to a double-porosity
conceptualization, using first the approach of Warren and Root [78] for thin-
film mass transfer between fractures and matrix, then the more physically
realistic matrix diffusion approach of Kazemi [42].

Double-porosity flow is a common and efficient conceptualization in fractured
rocks [2, 34, 21]. The medium is conceptualized as two communicating
physically overlapping continua including fractures with high permeability
(but little to no storage) and matrix or intact rock with significant storage
(but little to no flow) [5, 6]. Many extensions to the basic double-porosity
conceptual model exist, including multiple matrix or fracture porosities, and
different assumptions about the geometry or underlying physics governing flow
in the fractures or matrix [17, 49]. The Warren and Root [78] solution
simplifies the exchange between matrix and fractures to a mass-transfer thin-
film approximation, leading to numerous analytical solutions [2, 17]. It is
commonly used for this reason, even though it is well-known that spatial
pressure gradients in matrix blocks are important, as the matrix is low-
permeability and would therefore be expected to experience steep, slow-
changing gradients. A series representation of the Kazemi [42] solution is
used here, an extension of the multirate mass transfer model to double-
porosity flow [51]. The more physically correct (but more difficult to solve)
solution can be represented by an infinite series of porosities, which can
either represent an infinite number of Warren-Root type matrix porosities, or
if the coefficients are chosen specifically, a single Kazemi-type matrix
diffusion porosity. More recently, Wang et al [77] has developed a semi-
analytical solution for flow in a double-porosity formation, for the case when
non-Darcian flow is significant. Moutsopoulos et al [62] have provided
analytical and semi-analytical solutions for two classical problems in flow of
unconfined double-porosity aquifers, based on Moutsopoulos [61]. De-Smedt [26]
presented an analytical solution for flow in double-porosity media for
fractional flow dimensions, which is a generalization of De-Smedt [25]. Hayek
et al [37] presented a semi-analytical solution for flow due to pumping a
double-porosity aquifer via a constant-pressure boundary condition (without
wellbore storage) where permeability varied as a power law.

The fractal reservoir flow problem [16] is also analogous to the radially
variable properties approach presented here, but the governing equations of
the two problems are only equivalent when the spectral exponent ($\theta$ in
Chang and Yortsos [16]) in the fractal problem is zero. The fractal reservoir
governing equation is typically solved approximately, since the additional
terms due to non-zero spectral exponent in the governing equation do not
readily allow closed-form analytical solution.

In the next section, the governing equations and boundary conditions are
developed for the variable-dimension single-porosity flow problem [30]. This
solution is mapped onto the modified Bessel equation, allowing solution for
flow to both specified pressure (type-I) and specified flowrate with wellbore
storage (type-III). These more general single-porosity solutions are shown to
degenerate down to several well-known cases. The single-porosity solutions are
then extended to a simpler Warren-Root type double-porosity model for type-I
[37] and type-III (new) and then a new Kazemi type double-porosity model. The
Kazemi series solution approach is then summed analytically to arrive at a new
closed-form expression for the response in Laplace space, a solution that is
new for both graded and homogeneous domains. Finally, a summary and discussion
of limitations is given for the new solutions.

The approach taken here, representing the porosity and permeability of
fractured rocks as power-law distributions, was first developed by [28], and
first pursued by the author for applications in deep ($>3$ km) borehole
disposal of radioactive waste in basement rock [13, 52]. The approach is also
applicable to flow in salt surrounding excavations, like those in mine-by
experiments [71].

##  2 Development of Flow Problem

To introduce and contrast with the dual-porosity solution, the single-porosity
solution is developed first. To make a single solution for Cartesian linear,
cylindrical, and spherical geometries, a variable-dimension approach like
Barker [7] is used, including variable permeability and porosity, like Doe
[30]. The governing equation for slightly compressible time-dependent change
in pressure $p$ [$\mathrm{Pa}$] in a general 1D coordinate [7] is

$n(r)c\frac{\partial p}{\partial t}=\frac{1}{r^{m}}\frac{\partial}{\partial
r}\left[\frac{k(r)r^{m}}{\mu}\frac{\partial p}{\partial r}\right],$ (3)

where $c$ is bulk compressibility $\mathrm{[1/Pa]}$ and the dimensionless
parameter $m$ is 0 for a Cartesian strip, 1 for a cylinder, and 2 for a sphere
(i.e., $m=D-1$, where $D$ is the dimension). The derivative of the bracketed
term in (3) is expanded via chain rule; starting from (2),
$\frac{\mathrm{d}k}{\mathrm{d}r}=-\kappa k(r)/r$ is substituted with the
definitions of $k(r)$ and $n(r)$, to get

$n_{0}c\left(\frac{{r}}{r_{w}}\right)^{-\eta}\frac{\partial p}{\partial
t}=\frac{k_{0}}{\mu}\left(\frac{{r}}{r_{w}}\right)^{-\kappa}\left[\frac{m-\kappa}{r}\frac{\partial
p}{\partial r}+\frac{\partial^{2}p}{\partial r^{2}}\right].$ (4)

For converging radial flow in a semi-infinite domain, the relevant wellbore
boundary conditions are constant-pressure (type-I), constant-flux (type-II),
or constant-flux with wellbore storage (type-III in Laplace space). The
initial, far-field, and source borehole boundary conditions for a borehole in
an infinite symmetric domain are

$\displaystyle\mathrm{initial}\quad p({r},t=0)$ $\displaystyle=0$
$\displaystyle\mathrm{far-field}\quad p({r}\rightarrow\infty,t)$
$\displaystyle<\infty$ $\displaystyle\mathrm{wellbore\;type-I}\quad
p^{I}(r=r_{w},t)$ $\displaystyle=p_{1}(t);\quad\mathrm{or}$ (5)
$\displaystyle\mathrm{wellbore\;type-
II}\quad\left.\frac{A_{m}k_{0}}{\mu}\frac{\partial p^{II}(t)}{\partial
r}\right|_{r=r_{w}}$ $\displaystyle=Q(t);\quad\mathrm{or}$
$\displaystyle\mathrm{wellbore\;type-
III}\quad\frac{A_{m}k_{0}}{\mu}\left.\frac{\partial p^{III}(t)}{\partial
r}\right|_{r=r_{w}}$ $\displaystyle=Q(t)+\frac{A_{c}}{\rho g}\frac{\partial
p_{w}(t)}{\partial t},$

respectively. See Appendix A for definition of source borehole boundary
condition terms. These boundary conditions represent a homogeneous uniform
initial condition, a requirement that the solution remains finite at large
distance, and a specified pressure or pressure gradient at the source
($r=r_{w}$).

The Type-II boundary condition (specified flowrate) is a special case
($\sigma=0$) of the wellbore storage boundary condition (flowrate linearly
proportional to change in pressure), so it is not developed further.

###  2.1 Dimensional Analysis

A solution is derived for equation (4), using the approach of Doe [30], which
was based on analogy with the heat conduction literature [15]. Reducing the
governing equation (4) to dimensionless form using characteristic time,
$T_{c}=n_{0}cL_{c}^{2}\mu/k_{0}$, and characteristic length, $L_{c}=r_{w}$,
leads to

$r_{D}^{\kappa-\eta}\frac{\partial p_{D}}{\partial
t_{D}}=\frac{m-\kappa}{r_{D}}\frac{\partial p_{D}}{\partial
r_{D}}+\frac{\partial^{2}p_{D}}{\partial r_{D}^{2}},$ (6)

where the dimensionless quantities $r_{D}=r/L_{c}$, $t_{D}=t/T_{c}$, and
$p^{\\{I,III\\}}_{D}=p/p^{\\{I,III\\}}_{c}$ are used (see Table 2 for a
summary of dimensionless quantities).

The characteristic pressure change is given by $p^{I}_{c}=\hat{p}_{1}$, where
$p_{1}(t)=\hat{p}_{1}f_{t}$ separates the time-dependent specified pressure
into a constant characteristic pressure and a dimensionless variable time
behavior (for a constant specified pressure, $f_{t}=1$). The dimensionless
type-I initial and boundary conditions are

$\displaystyle p_{D}(r_{D},t_{D}=0)$ $\displaystyle=0$ $\displaystyle
p_{D}(r_{D}\rightarrow\infty,t_{D})$ $\displaystyle<\infty$ (7) $\displaystyle
p^{I}_{D}(r_{D}=1,t_{D})$ $\displaystyle=f_{t}.$

Using $p^{III}_{c}=\frac{r_{w}\hat{Q}\mu}{A_{m}k_{0}}$, where
$Q(t)=\hat{Q}f_{t}$ similarly separates the time-dependent volumetric flowrate
into a constant characteristic flowrate and a dimensionless time behavior. The
dimensionless type-III source borehole boundary condition is

$\left.\frac{\partial p^{III}_{D}}{\partial
r_{D}}\right|_{r_{D}=1}=f_{t}+\sigma\frac{\partial p^{III}_{D}}{\partial t},$
(8)

where $\sigma$ is a dimensionless wellbore storage coefficient (see Appendix
A) and the same initial and far-field conditions apply as the type-I case.

###  2.2 Laplace Transform

Taking the dimensionless Laplace transform
$\left(\bar{f}(s)=\int_{0}^{\infty}e^{-st_{D}}f(t_{D})\;\mathrm{d}t_{D}\right)$
of the governing partial differential equation (6) (without loss of generality
assuming zero initial condition) leads to the ordinary differential equation

$\frac{\mathrm{d}^{2}\bar{p}_{D}}{\mathrm{d}r_{D}^{2}}+\frac{m-\kappa}{r_{D}}\frac{\mathrm{d}\bar{p}_{D}}{\mathrm{d}r_{D}}-s\bar{p}_{D}r_{D}^{\kappa-\eta}=0,$
(9)

assuming $\kappa$, $\eta$, and $m$ are not functions of time, and $s$ is the
dimensionless Laplace transform parameter.

The transformed type-I and far-field boundary conditions (2.1) are

$\displaystyle\bar{p}_{D}(r_{D}\rightarrow\infty)$ $\displaystyle<\infty$ (10)
$\displaystyle\bar{p}^{I}_{D}(r_{D}=1)$ $\displaystyle=\bar{f}_{t},$

where $\bar{f}_{t}$ represents the Laplace transform of the boundary
condition’s time behavior. For a unit step change at $t=0$ (where $f_{t}=1$, a
typical assumption), $\bar{f}_{t}=\frac{1}{s}$. Other temporal behaviors are
simply handled, including a step change at a non-zero time, an exponentially
decaying source term, an arbitrary piecewise-constant or piecewise-linear
behavior, or a sinusoidal source term [46, 60].

The transformed wellbore-storage boundary condition is

$\left.\frac{\mathrm{d}\bar{p}^{III}_{D}}{\mathrm{d}r_{D}}\right|_{r_{D}=1}=\bar{f}_{t}+\sigma
s\bar{p}^{III}_{D},$ (11)

which now more clearly resembles a Type-III boundary condition.

###  2.3 Numerical Inverse Laplace Transform

The governing equations and associated boundary conditions are solved exactly
in Laplace space, then numerically inverted back to the time domain using one
of several viable approaches [47]. The equations were rapidly prototyped and
inverted using the Python library mpmath [41], which provides arbitrary
precision special functions and numerical inverse Laplace transform
algorithms. A Fortran program was also developed to facilitate plotting and
parameter estimation, implementing the inversion algorithm of de Hoog et al
[24]. Python and Fortran implementations of the solution are available at
https://github.com/klkuhlm/graded.

##  3 Solution of Flow Problem

###  3.1 Mapping onto Modified Bessel Equation

The governing ordinary differential equation (9) can be made equivalent to a
form of the modified Bessel equation after a change of variables first used by
Lommel [58] for the standard Bessel equation. Appendix B illustrates an
analogous change of variables to the modified Bessel equation. Comparing (9)
to this scaled version of the modified Bessel equation (41), they are
equivalent given the following correspondences

$\displaystyle\alpha=$ $\displaystyle\frac{1}{2}\left(\kappa-m+1\right)$
$\displaystyle\gamma=$ $\displaystyle\frac{1}{2}\left(\kappa-\eta+2\right)$
(12) $\displaystyle\nu=$
$\displaystyle\sqrt{\frac{\alpha^{2}}{\gamma^{2}}}=\frac{\kappa-m+1}{\kappa-\eta+2}$
$\displaystyle\beta$
$\displaystyle=\sqrt{\frac{s}{\gamma^{2}}}=\sqrt{\frac{4s}{\left(\kappa-\eta+2\right)^{2}}}.$

The transformed modified Bessel equation has the general solution (37)

$y=z^{\alpha}\left[A\mathrm{I}_{\nu}\left(\beta
z^{\gamma}\right)+B\mathrm{K}_{\nu}\left(\beta
z^{\gamma}\right)\right],\qquad\left(\gamma\neq 0\right),$ (13)

where $A$ and $B$ are constants determined by the boundary conditions and
$\mathrm{I}_{\nu}(z)$ and $\mathrm{K}_{\nu}(z)$ are the first- and second-kind
modified Bessel functions of non-integer order and real argument [59, 12, 69,
29].

The finiteness boundary condition (10) requires $A=0$ to keep the solution
finite as $r_{D}\rightarrow\infty$, since the first-kind modified Bessel
function grows exponentially with increasing real argument, leaving

$\bar{p}_{D}\left(r_{D}\right)=r_{D}^{\alpha}B\mathrm{K}_{\nu}\left(\beta
r_{D}^{\gamma}\right),$ (14)

which is not defined for $\gamma=0$ (i.e., $\kappa-\eta=-2$, which is
unrealistic because $\kappa$ is larger than $\eta$ for physical reasons), and
$B$ is determined by the Laplace-space source borehole boundary conditions.

###  3.2 Constant-Pressure (Type-I) at Borehole

The borehole boundary condition ($r_{D}=1$) for specified change in pressure
leads to the solution (the [78] double porosity solution for this wellbore
boundary condition is equivalent to Hayek et al [37])

$\bar{p}^{I}_{D}(r_{D})=\bar{f}_{t}r_{D}^{\alpha}\frac{\mathrm{K}_{\nu}\left(\beta
r_{D}^{\gamma}\right)}{\mathrm{K}_{\nu}\left(\beta\right)}$ (15)

and its radial gradient (i.e., proportional to flow of fluid into the
borehole)

$\frac{\mathrm{d}\bar{p}^{I}_{D}}{\mathrm{d}r_{D}}=\bar{f}_{t}r_{D}^{\alpha-1}\left[\left(\alpha-\gamma\nu\right)\frac{\mathrm{K}_{\nu}\left(\beta
r_{D}^{\gamma}\right)}{\mathrm{K}_{\nu}\left(\beta\right)}+\beta\gamma
r_{D}^{\gamma}\frac{\mathrm{K}_{\nu-1}\left(\beta
r_{D}^{\gamma}\right)}{\mathrm{K}_{\nu}\left(\beta\right)}\right],$ (16)

using a recurrence relationship for the derivative of the Bessel function in
terms of Bessel functions of adjacent orders [29, §10.29.2].

Restricting $\kappa\geq\eta$ (i.e., permeability decreases as fast or faster
than porosity), then $\gamma>0$ and $\alpha=\gamma\nu$ (for $\gamma<0$,
$\alpha-\gamma\nu=2\alpha$). This physically motivated restriction on
parameters simplifies (16) to

$\frac{\mathrm{d}\bar{p}^{I}_{D}}{\mathrm{d}r_{D}}=\sqrt{s}\bar{f}_{t}r_{D}^{\alpha+\gamma-1}\frac{\mathrm{K}_{\nu-1}\left(\beta
r_{D}^{\gamma}\right)}{\mathrm{K}_{\nu}\left(\beta\right)},$ (17)

since $\beta\gamma=\sqrt{s}$ for $\gamma>0$. When evaluated in the source
borehole ($r_{D}=1$), the solution simplifies further.

Figure 2 shows plots of the predicted pressure gradient at $r_{D}=1$ due to a
constant-pressure condition there (top row) and the predicted decrease in
pressure radially away from the boundary (values of $\eta$, $\kappa$, and $m$
for each simulation are listed in the caption and title of each figure). Both
rows of plots show the variability with the porosity exponent ($\eta$, given
by the line color) and the permeability exponent ($\kappa=\eta\tau$, given by
the line type). The same results are shown for Cartesian linear ($m=0$),
cylindrical ($m=1$), and spherical ($m=2$) geometries in three columns.

For a given set of parameters, a higher-dimensional domain (larger $m$) leads
to a slower drop in produced fluids at any time. The highest sustained
flowrate for all dimensions is achieved with constant properties in space
(i.e., the red curve $\eta=\kappa=0$). More negative exponents in the porosity
and permeability power-laws lead to more rapid decrease in flowrate, as the
contribution to flow from large radius vanishes when the exponent increases in
magnitude. These types of responses might be mis-interpreted as being
associated with lower permeability (which would also lead to a faster decrease
in flowrate) using a model with constant properties and a fixed dimension.

In the source well (top row of subplots), the effect of $\kappa$ is different
and are predicted to reverse between dimensions. For $\eta=3$ (black lines),
the $\kappa=\left\\{3,6,9\right\\}$ cases are swapped between $m=1$ and $m=2$.
For $\eta=2$ (blue lines), the $\kappa$ cases are swapped between $m=0$ and
$m=1$.

The bottom row of figures shows the predicted pressure with distance at
$t_{D}=10$. At locations away from the source well ($r_{D}>1$), changes in the
porosity exponent, $\eta$, have relatively less impact than changes in the
permeability exponent, $\kappa$ (different colored solid lines are close
together, while colored lines of different line type are widely separated).
The dimensionality ($m$) has a smaller effect at locations away from the
source borehole than it had on the gradient predicted at the source borehole.

Figure 2: Type-I flowrate (top row at $r_{D}=1$) and pressure (bottom row at
$r_{D}>1$ and $t_{D}=10$) solution at borehole for $m=0,1,2$ (Cartesian,
cylindrical, and spherical) and at different radial distances. Line color
indicates $\eta$; line type indicates $\kappa/\eta$. Line segments in top row
illustrate slopes of $1/2$, $1$, and $3/2$.

###  3.3 Constant-Flowrate with Wellbore Storage (Type-III)

The wellbore-storage boundary condition for the specified flowrate solution at
$r_{D}=1$ results in the general solution (that is new for any double-porosity
solution with power-law variation in material properties)

$\bar{p}^{III}_{D}(r_{D})=\bar{f}_{t}r_{D}^{\alpha}\frac{\mathrm{K}_{\nu}\left(\beta
r_{D}^{\gamma}\right)}{\left(\alpha-\gamma\nu+\sigma
s\right)\mathrm{K}_{\nu}\left(\beta\right)+\beta\gamma\mathrm{K}_{\nu-1}\left(\beta\right)},$
(18)

which can be simplified using $\alpha=\gamma\nu$ and $\beta\gamma=\sqrt{s}$ to

$\bar{p}^{III}_{D}(r_{D})=\bar{f}_{t}r_{D}^{\alpha}\frac{\mathrm{K}_{\nu}\left(\beta
r_{D}^{\gamma}\right)}{\sqrt{s}\mathrm{K}_{\nu-1}\left(\beta\right)+\sigma
s\mathrm{K}_{\nu}\left(\beta\right)}.$ (19)

Analogous to the results for the Type-I solution but only showing the $m=1$
and $m=2$ cases, Figure 3 shows the predicted pressure through time at the
boundary for a specified flowrate at the boundary. Figure 3 results are for no
wellbore storage ($\sigma=0$), while Figure 4 shows the same results with non-
zero wellbore storage (all model parameters listed in caption or title of each
figure). Wellbore storage is important at early time, leading to a smaller
predicted change in pressure, with the predicted response giving a
characteristic $1:1$ slope on log-log plots before formation storage
contributes significantly to the flow (i.e., pumping in a bathtub). Wellbore
storage makes more of a difference (i.e., shows a larger deviation from
$\sigma=0$ case) for larger $\eta$ (and $\kappa$, since $\kappa=2\eta$).

Figure 3: Type-II solution (Type-III with $\sigma=0$) at borehole for $m=1,2$
(cylindrical and spherical). Line color indicates $\eta$; line type indicates
$\kappa/\eta$.
Figure 4: Type-III solution at borehole ($r_{D}=1$), for $m=1,2$ (cylindrical
and spherical). Line color indicates $\eta$; line type indicates $\sigma$. All
curves for $\kappa/\eta=2$.

###  3.4 Parameter Combinations Yielding Simpler Solutions

When $\eta=\kappa=0$, permeability and porosity are constant in space; in this
case (9) simplifies to

$\frac{\mathrm{d}^{2}\bar{p}_{D}}{\mathrm{d}r_{D}^{2}}+\frac{m}{r_{D}}\frac{\mathrm{d}\bar{p}_{D}}{\mathrm{d}r_{D}}-s\bar{p}_{D}=0,$
(20)

which is the dimensionless form of the equation solved by Barker [7]. In this
case $\gamma=1$, $\alpha=(1-m)/2$, $\nu=\alpha$, and $\beta=\sqrt{s}$. The
solution in Laplace-space under these conditions becomes

$\bar{p}_{D}\left(r_{D}\right)=r_{D}^{\nu}B\mathrm{K}_{\nu}\left(\sqrt{s}r_{D}\right),$
(21)

which was found by Barker [7, Eqn. 15].

When $\eta=\kappa=m=0$ the time-domain solution simplifies to
$p_{D}(t)=1/\sqrt{\pi t}$, because $\nu=1/2$ and $\nu-1=-1/2$, the numerator
and denominator of (17) are equal since
$\mathrm{K}_{\nu}(z)\equiv\mathrm{K}_{-\nu}(z)$.

Another simplification occurs when $m=\kappa=\eta$, not necessarily zero. In
this case, the permeability and porosity decrease at the same rate radially
that the surface area of the domain grows in size ($A_{0}\propto 1$,
$A_{1}\propto r_{D}$, $A_{2}\propto r_{D}^{2}$), resulting in an equivalent
Cartesian coordinate system,

$\frac{\mathrm{d}^{2}\bar{p}_{D}}{\mathrm{d}r_{D}^{2}}-s\bar{p}_{D}=0,$ (22)

which has a solution in terms of $\sin(\sqrt{s}r_{D})$ and
$\cos(\sqrt{s}r_{D})$ or $\exp(\pm\sqrt{s}r_{D})$ and typically has an
explicit inverse Laplace transform. In this case $\alpha=\nu=1/2$, $\gamma=0$,
and $\beta=\sqrt{s}$.

When $\nu=n\pm\frac{1}{2}$ (for $n$ integer), the modified Bessel functions
become modified spherical Bessel functions [29, §10.47], and when
$\nu=\pm\frac{1}{3}$, they become Airy functions [29, §9.6]. These additional
special cases are not handled differently here (i.e., the more general
solution in terms of modified Bessel functions is still valid), since in the
case given here $\nu$ varies with $\kappa$, $\eta$, and $m$ (12).

##  4 Extension of Solution to Double Porosity

###  4.1 Mass-Transfer Coefficient Approximation

Beginning with the Warren and Root [78] formulation for double-porosity (i.e.,
high-conductance fractures and high-capacity matrix), the power-law
permeability and porosity distributions are incorporated. The equations for
double-porosity flow in the fractures and matrix are

$\displaystyle\frac{1}{r^{m}}\frac{\partial}{\partial
r}\left[\frac{k_{f}}{\mu}\frac{\partial p_{f}}{\partial r}\right]$
$\displaystyle=n_{r}c_{r}\frac{\partial p_{r}}{\partial
t}+n_{f}c_{f}\frac{\partial p_{f}}{\partial t}$
$\displaystyle\frac{\hat{\alpha}k_{r}}{\mu}\left(p_{f}-p_{r}\right)$
$\displaystyle=n_{r}c_{r}\frac{\partial p_{r}}{\partial t}$ (23)

where $\hat{\alpha}$ is the shape factor $\mathrm{[1/m^{2}]}$ of Warren and
Root [78], subscript $f$ indicates fracture, and subscript $r$ indicates
matrix (rock). The matrix equation does not involve a spatial gradient of
pressure, nor a matching of pressure and flux at the boundary, but simply a
difference between the fracture and matrix pressure (i.e., the mass transfer
coefficient approximation often used for heat transfer across thin films).
This behavior is sometimes referred to in the petroleum engineering literature
as “steady-state” flow between the fracture and matrix [21], but it also
represents one-dimensional diffusion in the matrix with a thin-film mass-
transfer approximation between the fracture and matrix reservoirs, analogous
to Newton’s law of cooling.

Substituting the permeability
$k_{i}=k_{i0}\left(\frac{r}{r_{w}}\right)^{-\kappa_{i}}$ and porosity
$n_{i}=n_{i0}\left(\frac{r}{r_{w}}\right)^{-\eta_{i}}$ ($i\in\\{f,r\\}$), then
converting to dimensionless form using an analogous approach to Warren and
Root [78], where $\omega=n_{f0}c_{f}/\left(n_{r0}c_{r}+n_{f0}c_{f}\right)$ is
the dimensionless fracture storage coefficient and
$\lambda=\hat{\alpha}k_{r}r_{w}^{2}/k_{f}$ is the dimensionless interporosity
exchange coefficient. Finally, taking the Laplace transform of both equations
results in the pair of ordinary differential equations

$\displaystyle\left[\frac{\mathrm{d}^{2}\bar{p}_{fD}}{\mathrm{d}r_{D}^{2}}+\frac{m-\kappa_{f}}{r_{D}}\frac{\mathrm{d}\bar{p}_{fD}}{\mathrm{d}r_{D}}\right]r^{-\kappa_{f}}$
$\displaystyle=(1-\omega)r_{D}^{-\eta_{r}}\bar{p}_{mD}s+\omega
r_{D}^{-\eta_{f}}\bar{p}_{fD}s$
$\displaystyle\lambda\left(\bar{p}_{fD}-\bar{p}_{rD}\right)r_{D}^{-\kappa_{r}}$
$\displaystyle=(1-\omega)r_{D}^{-\eta_{r}}\bar{p}_{rD}s.$ (24)

Solving for matrix pressure in the matrix equation,
$\bar{p}_{rD}=\bar{p}_{fD}\lambda
r_{D}^{-\kappa_{r}}/\left[(1-\omega)sr_{D}^{-\eta_{r}}+\lambda
r_{D}^{-\kappa_{r}}\right]$, and substituting this into the fracture equation
leads to a single equation solely in terms of dimensionless Laplace-domain
fracture pressure

$\displaystyle\left[\frac{\mathrm{d}^{2}\bar{p}_{fD}}{\mathrm{d}r_{D}^{2}}+\frac{m-\kappa_{f}}{r_{D}}\frac{\mathrm{d}\bar{p}_{fD}}{\mathrm{d}r_{D}}\right]r^{-\kappa_{f}}$
$\displaystyle=r_{D}^{-\eta_{r}}\bar{p}_{fD}\left\\{\frac{(1-\omega)sr_{D}^{-\kappa_{r}}\lambda}{(1-\omega)sr_{D}^{-\eta_{r}}+\lambda
r_{D}^{-\kappa_{r}}}\right\\}+\omega r_{D}^{-\eta_{f}}\bar{p}_{fD}s.$ (25)

To force the term in curly brackets in (25) to be independent of $r_{D}$,
$\kappa_{r}=\eta_{r}$ is assumed. Setting $\kappa_{r}$ and $\eta_{r}$ equal to
$\eta_{f}$ allows $r_{D}$ and $\bar{p}_{fD}$ to be similar form to previous
solutions. Simplifying the subsequent notation $\kappa_{f}\rightarrow\kappa$,
$\eta_{r}\rightarrow\eta$, and $\bar{p}_{fD}\rightarrow\bar{p}_{D}$ results in

$\displaystyle\frac{\mathrm{d}^{2}\bar{p}_{D}}{\mathrm{d}r_{D}^{2}}+\frac{m-\kappa}{r_{D}}\frac{\mathrm{d}\bar{p}_{D}}{\mathrm{d}r_{D}}$
$\displaystyle=r_{D}^{\kappa-\eta}\bar{p}_{D}\left[\frac{(1-\omega)s\lambda}{(1-\omega)s+\lambda}+\omega
s\right],$ (26)

which is the same form as (9). This solution corresponds to the same scaled
Bessel equation, with only the definition of $\beta$ changing to

$\beta_{WR}=\sqrt{\left[\frac{\lambda}{\lambda/(1-\omega)+s}+\omega\right]\frac{s}{\gamma^{2}}}.$
(27)

Any more general spatial behavior of matrix properties (e.g.,
$\eta_{r}\neq\kappa_{r}$) would not be solvable with the same approach. This
limitation still makes physical sense, as the the most important terms to vary
with space are the fracture permeability and the matrix storage. Setting
$\kappa=\eta=0$ and $m=1$ results in the Warren and Root [78] solution.

Figure 5 shows typical solution behaviors for the cylindrical ($m=1$) case for
Type-I and Type-II wellbore boundary conditions, for $\eta=3$ and $\kappa=6$.
Figure 6 shows behavior from the “middle” curve in Figure 5 ($\lambda=10^{-5}$
and $\omega=10^{-4}$), for a range of porosity and permeability exponents
similar to those shown in [78], listed in the figure caption.

Figure 5: Type-I flowrate solution at borehole (left) and Type-II solution for
pressure ($\sigma=0$, right), for $m=1$ (cylindrical). Line color indicates
$\lambda$; line type indicates $\omega$.

Figure 6: Type-I flowrate solution at borehole (left) and Type-II solution for
pressure ($\sigma=0$, right), for $m=1$ (cylindrical). All curves are for
$\lambda=10^{-5}$ and $\omega=10^{-4}$ (middle curves shown in Figure 5). Line
color indicates $\eta$; line type indicates $\kappa/\eta$.

###  4.2 Matrix Diffusion

The matrix diffusion problem of Kazemi [42] is more physically realistic [2,
21], but it is typically solved numerically or via late-time approximations
[27], rather than analytically like Warren and Root [78]. The series approach
of Kuhlman et al [51] is used here to represent matrix diffusion in a single
matrix continuum through the sum of an infinite series of Warren-Root matrix
continua, and the infinite sum is then analytically summed.

The generalization of (4.1) to multiple matrix continua starts with

$\displaystyle\frac{1}{r^{m}}\frac{\partial}{\partial
r}\left[\frac{k_{f}}{\mu}\frac{\partial p_{f}}{\partial r}\right]$
$\displaystyle=\sum_{j=1}^{N}n_{j}c_{j}\frac{\partial p_{j}}{\partial
t}+n_{f}c_{f}\frac{\partial p_{f}}{\partial t}$
$\displaystyle\frac{\hat{\alpha_{j}}k_{j}}{\mu}\left(p_{f}-p_{j}\right)$
$\displaystyle=n_{j}c_{j}\frac{\partial p_{j}}{\partial t}\qquad j=1,2,\dots
N,$ (28)

where $N$ is the number of matrix continua (one additional equation for each
continuum). Similarly taking the Laplace transform of this set of equations,
solving for $\bar{p}_{f}$, substituting the matrix equations into the fracture
equation, and simplifying the notation leads to

$\displaystyle\frac{\mathrm{d}^{2}\bar{p}_{D}}{\mathrm{d}r_{D}^{2}}+\frac{m-\kappa}{r_{D}}\frac{\mathrm{d}\bar{p}_{D}}{\mathrm{d}r_{D}}$
$\displaystyle=r_{D}^{\kappa-\eta}\bar{p}_{D}\omega s(1+\bar{g}),$ (29)

where

$\bar{g}=\sum_{j=1}^{N}\frac{\hat{\xi}_{j}u_{j}}{s+u_{j}}$ (30)

is a matrix memory kernel [35], $\hat{\xi}_{j}$ is related to the storage
properties of each matrix continuum (analogous to $\omega$ of Warren and Root
[78]), and $u_{j}$ is related to the interporosity flow coefficient of each
matrix continuum (analogous to $\lambda$ of Warren and Root [78]). The
Laplace-space memory kernel approach is flexible, and is used elsewhere in
hydrology and reservoir engineering [38, 36, 67]. Equation (29) can be
simplified to Warren and Root [78] with a particular choice of $\bar{g}$ and
$N=1$, and to the solution for a triple-porosity reservoir [19] with a
different choice of $\bar{g}$ and $N=2$ [51].

When $N\rightarrow\infty$ in (30), the it is more convenient to specify the
mean and variance of the parameter distributions than the individual
parameters associated with each porosity. Several different distributions are
possible [35]. In the form presented by [51], the parameters are specified as
the infinite series

$u_{j}=\frac{(2j-1)^{2}\pi^{2}\lambda}{4(1-\omega)}\qquad\hat{\xi}_{j}=\frac{8(1-\omega)}{(2j-1)^{2}\omega\pi^{2}}\qquad
j=1,2,\dots N\rightarrow\infty$ (31)

which leads to the Kazemi [42] solution for matrix diffusion. The parameters
$\lambda$ and $\omega$ have the same definitions as in Warren and Root [78].

Setting $\kappa=\eta=0$ results in the solution of Kuhlman et al [51]. The new
governing equation is the same form and the modified Bessel function solution,
only requiring re-definition of $\beta$ as

$\displaystyle\beta_{KZ}=\sqrt{\left[\sum_{j=1}^{N}\frac{\omega\hat{\xi}_{j}u_{j}}{u_{j}+s}+\omega\right]\frac{s}{\gamma^{2}}},\qquad
N\rightarrow\infty.$ (32)

Substituting the definitions of $u$ and $\hat{\xi}$ from (31) and simplifying
leads to

$\displaystyle\beta_{KZ}=\sqrt{\left[\sum_{j=1}^{N}\frac{2\lambda}{W^{2}_{j}\lambda/(1-\omega)+s}+\omega\right]\frac{s}{\gamma^{2}}},\qquad
N\rightarrow\infty,$ (33)

where $W_{j}=\pi(2j-1)/2$. This is similar in form to (27) but the term in the
denominator grows as the index increases, illustrating how the series solution
approximates the Kazemi [42] solution through an infinite series of modified
Warren and Root [78] matrix porosities.

Further simplifying the approach of Kuhlman et al [51], the infinite series in
(33) can be evaluated in closed form using residue methods [79], resulting in

$\beta_{KZ}=\sqrt{\left[\sqrt{\frac{\lambda(1-\omega)}{s}}\tanh\left(\sqrt{\frac{s(1-\omega)}{\lambda}}\right)+\omega\right]\frac{s}{\gamma^{2}}},$
(34)

where $\tanh(\cdot)$ is the hyperbolic tangent. This closed-form expression
derived here is more accurate and numerically more efficient than truncating
or accelerating the infinite series in (32), which is an improvement over the
series presented in Kuhlman et al [51] for graded or homogeneous domains.

Figure 7: Type-II solution for pressure at source borehole ($\sigma=0$), for
$m=1$ (cylindrical) for different number of terms. All curves are for
$\lambda=10^{-5}$, $\omega=10^{-4}$, $\kappa=\eta=0$.

Figure 7 illustrates the transition from the Warren and Root [78] ($N=1$) to
the Kazemi [42] series approximation for increasing terms
($N=\left\\{2,10,100,1000\right\\}$, heavy colored solid lines) and the
expression for the infinite sum (34) (heavy black dashed line) for flow to a
specified flux (type-II, $\sigma=0$) cylindrical ($m=1$) borehole of constant
material properties ($\kappa=\eta=0$). The bounding Theis [74] behavior is
shown for the fracture and matrix compressibilities (thin red dashed lines).

##  5 Applications and Limitations

A general converging radial flow solution for specified flowrate or specified
wellhead pressure was derived for domains with power-law variability in
porosity and permeability due to damage. The single-porosity version has
already been presented by Doe [30], and a solution for constant-pressure
condition without wellbore storage was derived by [37], but the specified-
flowrate double-porosity solution with wellbore storage presented here is new.
The infinite series approximation to Kazemi was summed analytically, resulting
in a new closed-form expression of the series presented in Kuhlman et al [51],
which is an improvement for both graded and homogeneous properties. The newly
developed analytical solutions are more general (i.e., several existing
solutions are special cases of the new solution) and include more behaviors
typical in well-test solutions (i.e., wellbore storage, positive skin, double
porosity), while still being straightforward and parsimonious (i.e., as few
free parameters as possible) in their implementation.

The basic flow solution assumes linear single-phase flow of a fluid in a
slightly compressible formation. The double-porosity solution assumes the
fractures are high permeability, with low storage capacity, while the matrix
(i.e., intact rock between fractures) is high storage capacity with low
permeability. These assumptions are representative for analytical solutions to
subsurface porous media flow problems in the hydrology and petroleum
engineering literature, and are shared by the solutions of [7], [30], [78],
[42], and [51].

To apply this analytical solution to observed data, either observed data would
be transformed into dimensionless space, or the analytical solution could be
transformed to dimensional space, then a parameter estimation routine would be
used to minimize the model-data misfit, and possibly explore the uncertainty
or uniqueness of the solution. The solution method developed to solve these
solutions uses numerical inverse Laplace transforms and runs quickly enough to
be used in parameter estimation (e.g., Monte Carlo methods that require
hundreds of thousands of evaluations).

The analytical solution might be of most use with parameter estimation to fit
observations, but the non-uniqueness of the curves may make estimation of
unique physical parameters difficult, without further physical or site-
specific constraints. Realistically, the parameters in the Bessel equation may
be estimable (i.e., $\alpha$, $\beta$, $\gamma$, and $\nu$ defined in (12)),
but without defining the flow dimension ($m$) or the relationship between the
porosity and permeability exponents ($\tau=\kappa/\eta$), it may be difficult
to identify all the parameters from data alone, since many the curves have
similar shapes, unlike classical Type curves [10].

$A_{c}$ | borehole cross-sectional area | $\mathrm{m^{2}}$
---|---|---
$A_{m}$ | borehole cylindrical surface area | $\mathrm{m^{2}}$
$c$ | bulk compressibility | $\mathrm{1/Pa}$
$f_{t}$ | time variability | $\mathrm{-}$
$g$ | gravitational acceleration | $\mathrm{m/s^{2}}$
$h$ | hydraulic head | $\mathrm{m}$
$k$ | permeability | $\mathrm{m^{2}}$
$L_{c}$ | characteristic length $(r_{w})$ | $\mathrm{m}$
$m$ | dimension $(D-1)$ | $\mathrm{-}$
$n$ | porosity | $\mathrm{-}$
$p$ | change in pressure | $\mathrm{Pa}$
$s$ | Laplace transform parameter | $\mathrm{-}$
$Q$ | volumetric flowrate | $\mathrm{m^{3}/s}$
$r$ | distance coordinate | $\mathrm{m}$
$r_{w}$ | borehole or excavation radius | $\mathrm{m}$
$\hat{\alpha}$ |  Warren and Root [78] shape factor | $\mathrm{1/m^{2}}$
$\eta$ | porosity power-law exponent | $\mathrm{-}$
$\kappa$ | permeability power-law exponent | $\mathrm{-}$
$\rho$ | fluid density | $\mathrm{kg/m^{3}}$
$\mu$ | fluid viscosity | $\mathrm{Pa\cdot s}$
Table 1: Physical Properties and Parameters $p_{D}$ | scaled pressure | $p/p_{c}$
---|---|---
$t_{D}$ | scaled time | $tk_{0}/n_{0}cL_{c}^{2}\mu$
$r_{D}$ | scaled distance | $r/L_{c}$
$\lambda$ | interporosity exchange coefficient | $\hat{\alpha}k_{r}r_{w}^{2}/k_{f}$
$\sigma$ | wellbore storage coefficient | $A_{c}/(r_{w}n_{0}c\rho gA_{m})$
$\omega$ | fracture storage coefficient | $n_{f0}c_{f}/(n_{r0}c_{r}+n_{f0}c_{f})$
Table 2: Dimensionless Quantities

## Statements and Declarations

### Funding

The author thanks the U.S. Department of Energy Office of Nuclear Energy’s
Spent Fuel and Waste Science and Technology program for funding.

### Conflicts of Interest

The author has no competing interests to declare.

### Availability of Data and Material

No data or materials were used by the author in the preparation of the
manuscript.

### Code Availability

The source code of Fortran and Python implementations of the program are
available from the author upon request.

## Acknowledgments

This paper describes objective technical results and analysis. Any subjective
views or opinions that might be expressed in the paper do not necessarily
represent the views of the U.S. Department of Energy or the United States
Government.

This article has been authored by an employee of National Technology &
Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the
U.S. Department of Energy (DOE). The employee owns all right, title and
interest in and to the article and is solely responsible for its contents. The
United States Government retains and the publisher, by accepting the article
for publication, acknowledges that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this article or allow others to do so, for United States
Government purposes. The DOE will provide public access to these results of
federally sponsored research in accordance with the DOE Public Access Plan
https://www.energy.gov/downloads/doe-public-access-plan.

The author thanks Tara LaForce from Sandia for technically reviewing the
manuscript.

##  6 Appendix A: Wellbore Storage Boundary Condition

The wellbore-storage boundary condition accounts for the storage in the finite
borehole arising from the mass balance
$Q_{\mathrm{in}}-Q_{\mathrm{out}}=A_{c}\frac{\partial h_{w}}{\partial t}$.
$Q_{\mathrm{in}}$ $\mathrm{[m^{3}/s]}$ is volumetric flow into the borehole
from the formation, $Q_{\mathrm{out}}$ is possibly time-variable flow out of
the well through the pump ($Q(t)$ $\mathrm{[m^{3}/s]}$), and $\frac{\partial
h_{w}}{\partial t}$ is the change in hydraulic head [m]
($h_{w}=\frac{p_{w}}{\rho g}+z$) of water standing in the borehole through
time, $p_{w}$ is change in pressure [Pa] of water in the borehole, $\rho$ is
fluid density $\mathrm{[kg/m^{3}]}$, $z$ is an elevation datum [m], and $g$ is
gravitational acceleration $\mathrm{[m/s^{2}]}$. $A_{c}$ is the cross-
sectional surface area of the pipe, sphere or box providing storage (it may be
a constant or a function of elevation); for a typical pipe, it becomes
$A_{c}=\pi r_{c}^{2}$, where $r_{c}$ is the radius of the casing where the
water level is changing. The mass balance is then

$\frac{A_{m}k_{0}}{\mu}\left.\frac{\partial p}{\partial
r}\right|_{r=r_{w}}-Q(t)=\frac{A_{c}}{\rho g}\frac{\partial p_{w}}{\partial
t},$ (35)

where $A_{m}$ is the area of the borehole communicating with the formation.
For the integer $m$ considered here these are $A_{0}=b^{2}$, $A_{1}=2\pi
r_{w}b$, $A_{2}=4\pi r_{w}^{2}$ ($b$ is a length independent of the borehole
radius).

Assuming the change in water level in the borehole ($h_{w}=p_{w}/\left(\rho
g\right)$) is equal to the change in formation water level ($h=p/\left(\rho
g\right)$), this can be converted into dimensionless form as

$\left.\frac{\partial p_{D}}{\partial
r_{D}}\right|_{r_{D}=1}-f_{t}=\sigma\frac{\partial p_{D}}{\partial t},$ (36)

where $\sigma=A_{c}/\left(r_{w}n_{0}c\rho gA_{m}\right)$ is a dimensionless
ratio of formation to wellbore storage; $\sigma\rightarrow 0$ is an
infinitesimally small well with only formation response, while
$\sigma\rightarrow\infty$ is a well with no formation response (i.e., a
bathtub).

##  7 Appendix B: Transformation of Modified Bessel Equation

Following the approach of Bowman [12], alternative forms of the Bessel
equation are found, this approach is a simplification of the original approach
of Lommel [58]. An analogous approach is applied here to “back into” the
desired modified Bessel equation. The equation satisfied by the pair of
functions

$\displaystyle y_{1}$ $\displaystyle=x^{\alpha}\mathrm{I}_{\nu}\left(\beta
x^{\gamma}\right),$ $\displaystyle y_{2}$
$\displaystyle=x^{\alpha}\mathrm{K}_{\nu}\left(\beta x^{\gamma}\right)$ (37)

is sought, where $\alpha$, $\beta$, $\gamma$, and $\nu$ are constants. Using
the substitutions $\zeta=yx^{-\alpha}$ and $\xi=\beta x^{\gamma}$ gives
$\zeta_{1}=\mathrm{I}_{\nu}\left(\xi\right)$ and
$\zeta_{2}=\mathrm{K}_{\nu}\left(\xi\right)$, which are the two solutions to
the modified Bessel equation [29, §10.25.1],

$\xi\frac{\mathrm{d}}{\mathrm{d}\xi}\left(\xi\frac{\mathrm{d}\zeta}{\mathrm{d}\xi}\right)-(\xi^{2}+\nu)\zeta=0.$
(38)

Given

$\xi\frac{\mathrm{d}}{\mathrm{d}\xi}\left(\xi\frac{\mathrm{d}\zeta}{\mathrm{d}\xi}\right)=\frac{x}{\gamma^{2}}\frac{\mathrm{d}}{\mathrm{d}x}\left(x\frac{\mathrm{d}\zeta}{\mathrm{d}x}\right),$
(39)

and

$x\frac{\mathrm{d}}{\mathrm{d}x}\left(x\frac{\mathrm{d}\zeta}{\mathrm{d}x}\right)=\frac{y^{\prime\prime}}{x^{\alpha-2}}-\frac{\left(2\alpha-1\right)y^{\prime}}{x^{\alpha-1}}+\frac{\alpha^{2}y}{x^{\alpha}},$
(40)

the standard-form equation satisfied by $y$ is

$y^{\prime\prime}+\left(1-2\alpha\right)y^{\prime}+\frac{\alpha^{2}y}{x^{\alpha}}-\left(\beta^{2}\gamma^{2}x^{2\gamma-2}-\frac{\alpha^{2}-\nu^{2}\gamma^{2}}{x^{2}}\right)y=0.$
(41)

This equation can be compared to the Laplace-space ordinary differential
equation (9), allowing direct use of the product of powers and modified Bessel
function (37) as solutions (13).

## References

  *   * AbuAisha and Billiotte [2021] AbuAisha M, Billiotte J (2021) A discussion on hydrogen migration in rock salt for tight underground storage with an insight into a laboratory setup. Journal of Energy Storage 38:102958
  * Aguilera [1980] Aguilera R (1980) Naturally Fractured Reservoirs. PennWell Books
  * Alkan [2009] Alkan H (2009) Percolation model for dilatancy-induced permeability of the excavation damaged zone in rock salt. International Journal of Rock Mechanics & Mining Sciences 46(4):716–724
  * Aydan et al [1993] Aydan Ö, Akagi T, Kawamoto T (1993) The squeezing potential of rocks around tunnels; theory and prediction. Rock Mechanics and Rock Engineering 26(2):137–163
  * Barenblatt and Zheltov [1960] Barenblatt GI, Zheltov YuP (1960) Basic flow equations for homogeneous fluids in naturally fractured rocks. Doklady Akademii Nauk SSSR 132(3):545–548
  * Barenblatt et al [1960] Barenblatt GI, Zheltov YuP, Kochina IN (1960) Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata]. Prikladnaya Matematika i Mekhanika 24(5):852–864
  * Barker [1988] Barker JA (1988) A generalized radial flow model for hydraulic tests in fractured rock. Water Resources Research 24(10):1796–1804
  * Beauheim and Roberts [2002] Beauheim RL, Roberts RM (2002) Hydrology and hydraulic properties of a bedded evaporite formation. Journal of Hydrology 259(1):66–88
  * Beauheim et al [2004] Beauheim RL, Roberts RM, Avis JD (2004) Well testing in fractured media: Flow dimension and diagnostic plots. Journal of Hydraulic Research 42:69–76
  * Bourdet et al [1989] Bourdet D, Ayoub JA, Pirard YM (1989) Use of pressure derivative in well-test interpretation. SPE Formation Evaluation 4(2):293–302
  * Bowman et al [2013] Bowman DO, Roberts RM, Holt RM (2013) Generalized radial flow in synthetic flow systems. Groundwater 51(5):768–774
  * Bowman [1958] Bowman F (1958) Introduction to Bessel Functions. Dover
  * Brady et al [2017] Brady PV, Freeze GA, Kuhlman KL, et al (2017) Deep borehole disposal of nuclear waste: US perspective. In: Geological Repository Systems for Safe Disposal of Spent Nuclear Fuels and Radioactive Waste. Elsevier, p 89–112
  * Carman [1937] Carman PC (1937) Fluid flow through granular beds. Transactions of the Institution of Chemichal Engineers 15:150–167
  * Carslaw and Jaeger [1959] Carslaw HS, Jaeger JC (1959) Conduction of Heat in Solids, 2nd edn. Oxford
  * Chang and Yortsos [1990] Chang J, Yortsos YC (1990) Pressure-transient analysis of fractured reservoirs. SPE Formation Evaluation 5(1):31–38
  * Chen [1989] Chen ZX (1989) Transient flow of slightly compressible fluids through double-porosity, double-permeability systems – a state-of-the-art review. Transport in Porous Media 4(2):147–184
  * Cho et al [2013] Cho WJ, Kim JS, Lee C, et al (2013) Gas permeability in the excavation damaged zone at KURT. Engineering Geology 164(17):222–229
  * Clossman [1975] Clossman PJ (1975) An aquifer model for fissured reservoirs. Society of Petroleum Engineers Journal 15:385–398
  * Cosenza [1996] Cosenza Ph (1996) Coupled effects between mechanical behavior and mass transfer phenomena in rock salt. PhD thesis, Ecole Polytechnique
  * Da Prat [1990] Da Prat G (1990) Well Test Analysis for Fractured Reservoir Evaluation. No. 27 in Developments in Petroleum Science, Elsevier
  * David et al [1994] David C, Wong TF, Zhu W, et al (1994) Laboratory measurement of compaction-induced permeability change in porous rocks: Implications for the generation and maintenance of pore pressure excess in the crust. PAGEOPH 143(1/2/3):425–456
  * Davies and Bernier [2005] Davies C, Bernier F (eds) (2005) Impact of the Excavation Disturbed or Damaged Zone (EDZ) on the Performance of Radioactive Waste Geological Repositories, EUR 21028 EN, Brussels, Belgium
  * de Hoog et al [1982] de Hoog FR, Knight JH, Stokes AN (1982) An improved method for numerical inversion of Laplace transforms. SIAM Journal of Scientific and Statistical Computing 3(3):357–366
  * De-Smedt [2011] De-Smedt F (2011) Analytical solution for constant-rate pumping test in fissured porous media with double-porosity behavior. Transport in Porous Media 88:479–489
  * De-Smedt [2022] De-Smedt F (2022) Analytical solution for fractional well flow in a double-porosity aquifer with fraction transient exchange between matrix and fractures. Water 14(456)
  * De Swaan [1976] De Swaan AO (1976) Analytic solutions for determining naturally fractured reservoir properties by well testing. Society of Petroleum Engineers Journal 16(03):117–122
  * Delay et al [2007] Delay F, Kaczmaryk A, Ackerer P (2007) Inversion of interference hydraulic pumping tests in both homogeneous and fractal dual media. Advances in Water Resources 30:314–334
  * DLMF [2023] DLMF (2023) NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.1.12 of 2023-12-15, URL http://dlmf.nist.gov/, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * Doe [1991] Doe TW (1991) Fractional dimension analysis of constant-pressure well tests. In: SPE Annual Technical Conference and Exhibition, pp 461–467
  * Ferroud et al [2018] Ferroud A, Chesnaux R, Rafini S (2018) Insights on pumping well interpretation from flow dimension analysis: The learnings of a multi-context field database. Journal of Hydrology 556:449–474
  * Gelbard [1992] Gelbard F (1992) Exact analysis of a two-dimensional model for brine flow to a borehole in a disturbed rock zone. Tech. Rep. SAND92–1303, Sandia National Laboratories, Albuquerque, NM
  * Ghazvinian [2015] Ghazvinian E (2015) Fracture initiation and propagation in low porosity crystalline rocks: Implications for excavation damage zone (EDZ) mechanics. PhD thesis, Queen’s University, Kingston, Ontario, Canada
  * van Golf-Racht [1982] van Golf-Racht TD (1982) Fundamentals of Fractured Reservoir Engineering. No. 12 in Developments in Petroleum Science, Elsevier
  * Haggerty and Gorelick [1995] Haggerty R, Gorelick SM (1995) Multiple-rate mass-transfer for modeling diffusion and surface-reactions in media with pore-scale heterogeneity. Water Resources Research 31(10):2383–2400
  * Haggerty et al [2000] Haggerty R, McKenna SA, Meigs LC (2000) On the late-time behavior of tracer test breakthrough curves. Water Resources Research 36(12):3467–3479
  * Hayek et al [2018] Hayek M, Younes A, Zouali J, et al (2018) Analytical solution and Bayesian inference for interference pumping tests in fractal dual-porosity media. Computational Geosciences 22:413–421
  * Herrera and Yates [1977] Herrera I, Yates R (1977) Integrodifferential equations for systems of leaky aquifers and applications 3. a numerical method of unlimited applicability. Water Resources Research 13(4):725–732
  * Hou [2003] Hou Z (2003) Mechanical and hydraulic behavior of rock salt in the excavation disturbed zone around underground facilities. International Journal of Rock Mechanics & Mining Sciences 40(5):725–738
  * Hudson et al [2009] Hudson JA, Bäckström A, Rutqvist J, et al (2009) Characterising and modelling the excavation damaged zone in crystalline rock in the context of radioactive waste disposal. Environmental Geology 57:1275–1297
  * Johansson et al [2017] Johansson F, Steinberg V, Kirpichev SB, et al (2017) mpmath: a Python library for arbitrary-precision floating-point arithmetic. 10.5281/zenodo.1476882, URL https://doi.org/10.5281/zenodo.1476882
  * Kazemi [1969] Kazemi H (1969) Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution. Society of Petroleum Engineers Journal 9(04):451–462
  * Kim et al [2012] Kim HM, Rutqvist J, Jeong JH, et al (2012) Characterizing excavation damaged zone and stability of pressurized lined rock caverns for underground compressed air energy storage. Rock Mechanics and Rock Engineering 46:1113–1124
  * Kozeny [1927] Kozeny J (1927) Uber kapillare leitung der wasser in boden. Royal Academy of Science, Vienna, Proc Class I 136:271–306
  * Kranz et al [1979] Kranz RL, Frankel AD, Engelder T, et al (1979) The permeability of whole and jointed Barre granite. International Journal of Rock Mechanics Mineral Sciences and Geomechanical Abstracts 16:225–234
  * Kruseman and de Ridder [1994] Kruseman GP, de Ridder NA (1994) Analysis and evaluation of pumping test data. Tech. Rep. Publication 47, International Institute for Land Reclamation and Improvement, Wageningen, The Netherlands
  * Kuhlman [2013] Kuhlman KL (2013) Review of inverse Laplace transform algorithms for Laplace-space numerical approaches. Numerical Algorithms 63(2):339–355
  * Kuhlman [2014] Kuhlman KL (2014) Summary results for brine migration modeling performed by LANL, LBNL, and SNL for the UFD program. Tech. Rep. SAND2014–18217R, Sandia National Laboratories, Albuquerque, NM, 10.2172/1163122
  * Kuhlman and Heath [2021] Kuhlman KL, Heath JE (2021) Multicontinuum flow models for assessing two-phase flow in containment science. Tech. Rep. SAND2021-7191, Sandia National Laboratories, Albuquerque, NM, 10.2172/1809129
  * Kuhlman and Matteo [2018] Kuhlman KL, Matteo EN (2018) Porosity and permeability: Literature review and summary. In: SaltMechIX. Federal Institute for Geosciences and Natural Resources (BGR), Hannover, Germany, pp 1–13
  * Kuhlman et al [2015] Kuhlman KL, Malama B, Heath JE (2015) Multiporosity flow in fractured low-permeability rocks. Water Resources Research 51(2):848–860
  * Kuhlman et al [2019] Kuhlman KL, Hardin EL, Rigali MJ (2019) Deep borehole laboratory and borehole testing strategy: Generic drilling and testing plan. Tech. Rep. SAND2019–1896, Sandia National Laboratories, Albuquerque, NM, 10.2172/1497220
  * Lai [1971] Lai CS (1971) Fluid flow through rock salt under various stress states. PhD thesis, Michigan State University, East Lansing, MI
  * Le Borgne et al [2004] Le Borgne T, Bour O, De Dreuzy J, et al (2004) Equivalent mean flow models for fractured aquifers: Insights from a pumping tests scaling interpretation. Water Resources Research 40(3)
  * Lin et al [2016] Lin YC, Yang SY, Fen CS, et al (2016) A general analytical model for pumping tests in radial finite two-zone confined aquifers with Robin-type outer boundary. Journal of Hydrology 540:1162–1175
  * Liu [2014] Liu HH (2014) Non-Darcian flow in low-permeability media: key issues related to geological disposal of high-level nuclear waste in shale formations. Hydrogeology Journal 22(7):1525–1534
  * Liu [2017] Liu HH (2017) Fluid Flow in the Subsurface. Springer, New York, NY
  * Lommel [1868] Lommel E (1868) Studien Über die Bessel’schen Functionen. von B.G. Teubner
  * McLachlan [1955] McLachlan NW (1955) Bessel Functions for Engineers, 2nd edn. Oxford
  * Mishra et al [2013] Mishra PK, Vessilinov V, Gupta H (2013) On simulation and analysis of variable-rate pumping tests. Groundwater 51(3):469–473
  * Moutsopoulos [2021] Moutsopoulos KN (2021) A simple model for the simulation of the flow behavior in unconfined double porosity aquifers. Journal of Hydrology 596(126076)
  * Moutsopoulos et al [2022] Moutsopoulos KN, Papaspyros JN, Fahs M (2022) Approximate solutions for flow in unconfined double porosity aquifers. Journal of Hydrology 615(128679)
  * Pasandi et al [2008] Pasandi M, Samani N, Barry DA (2008) Effect of wellbore storage and finite thickness skin on flow to a partially penetrating well in a preatic aquifer. Advances in Water Resources 31:383–398
  * Perras and Diederichs [2016] Perras MA, Diederichs MS (2016) Predicting excvation damage zone depths in brittle rocks. Journal of Rock Mechanics and Geotechnical Engineering 8(1):60–74
  * Perras et al [2010] Perras MA, Diederichs MS, Lam T (2010) A review of excavation damage zones in sedimentary rocks with emphasis on numerical modelling for EDZ definition. In: Proceedings of the 63rd Canadian Geotechnical Conference, Calgary, Canada, pp 742–750
  * Reynolds and Gloyna [1960] Reynolds TD, Gloyna EF (1960) Reactor fuel waste disposal project: Permeability of rock salt and creep of underground salt cavities. Tech. Rep. TID-12383, Atomic Energy Commission Division of Technical Information, Oak Ridge, TN
  * Schumer et al [2003] Schumer R, Benson DA, Meerschaertm MM, et al (2003) Fractal mobile/immobile solute transport. Water Resources Research 39(10):1296
  * Shen et al [2011] Shen B, Stephansson O, Rinne M, et al (2011) FRACOD modeling of rock fracturing and permeability change in excavation-damaged zones. International Journal of Geomechanics 11(4):302–313
  * Spanier and Oldham [1987] Spanier J, Oldham KB (1987) An Atlas of Functions. Hemisphere
  * Stormont and Fuenkajorn [1994] Stormont JC, Fuenkajorn K (1994) Dilation-induced permeability changes in rock salt. In: International Conference on Computer Methods and Advances in Geomechanics, Morgantown, West Virginia, pp 1–5
  * Stormont et al [1991] Stormont JC, Howard CL, Daemen JJK (1991) Changes in rock salt permeability due to nearby excavation. In: 32nd US Symposium on Rock Mechanics. OnePetro, Norman, Oklahoma, pp 1–10
  * Streltsova [1988] Streltsova TD (1988) Well Testing in Heterogeneous Formations. Exxon Monograph, John Wiley & Sons
  * Tao et al [2021] Tao J, Shi AC, Li HT, et al (2021) Thermal-mechanical modelling of rock response and damage eveolution during excavation in prestressed goethermal deposits. International Journal of Rock Mechanics and Mining Sciences 147:104913
  * Theis [1935] Theis CV (1935) The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using ground-water storage. Transactions, American Geophysical Union 16:519–524
  * Tsang et al [2005] Tsang CF, Bernier F, Davies C (2005) Geohydromechanical processes in the excavation damaged zone in crystalline rock, rock salt, and indurated and plastic clays–in the context of radioactive waste disposal. International Journal of Rock Mechanics and Mining Sciences 42(1):109–125
  * Wallace et al [1990] Wallace M, Pietz J, Case J, et al (1990) Coupled fluid-flow modeling of brines flowing through deforming salt around the excavations for the Waste Isolation Pilot Plant (WIPP), in the Permian Salado Formation. In: Waste Management 90, Tucson, AZ, pp 1–8
  * Wang et al [2021] Wang Y, Zhan H, Huang K, et al (2021) Identification of non-Darcian flow effect in double-porosity fracture aquifer based on multi-well pumping test. Journal of Hydrology 600(126541)
  * Warren and Root [1963] Warren J, Root P (1963) The behavior of naturally fractured reservoirs. Society of Petroleum Engineers Journal 3(03):245–255
  * Wolfram Research, Inc. [2021] Wolfram Research, Inc. (2021) Mathematica, Version 12.3. URL https://www.wolfram.com/mathematica, champaign, IL