The Theis solution for subdiffusive flow in rocks

The central contribution of this work is the development of a “master” solution similar to the Theis solution to evaluate well responses under subdiffusive flow. Models based on subdiffusion employ fractional constitutive laws, a redefinition of Darcy’s law. Subdiffusive models discussed here are particularly useful to address situations where the internal architecture of the geological medium, such as fluvial and fractured systems, matters and where the existence of topological, geometrical and spatial influences result in distorted flow paths and a loss in connectivity. The developed solution provides the means for addressing these ends.


Introduction
Forecasts, assessments and evaluations of subsurface resources whether they be hydrocarbon production or geothermal extraction (Arps, 1945;Fetkovich, 1980;Whiting and Ramey, 1969), sequestration and waste disposal technologies (Bodvarsson et al. 1999;IPCC, 2005) or groundwater systems require discernment of flow and transport mechanisms. The Theis (1935) solution has been the linchpin to assess the flow properties of the reservoir rock (aquifer). Both its adequacy and its limits have been of interest and tested over the years by statistical techniques because it ignores the influence of spatial variability within its support volume (Romeu and Noetinger, 1995;Sanchez-Vila et al., 2006;Schad and Teutsch, 1994;Warren and Price, 1961). Ignoring internal architecture, geologic structure and internal imprint such as stratification, intercalations and barriers that produce distorted flow paths, however, results in a serious underestimation of the resource (Raghavan, 2004(Raghavan, , 2009bThomas et al., 2005). Other approaches that account for spatial variability have represented transient diffusion by representing the porous medium as a fractal object (Camacho-Velázquez et al., 2008;Chang and Yortsos, 1990;Raghavan, 2009a). This approach is particularly fruitful because it has been suggested (Kim et al., 2015;Raghavan, 2011) and shown (Benson et al., 2004;Chu et al., 2017Chu et al., , 2018 that there is an interdependence between the topological and geometric properties of the reservoir rock and the physics of fluid movement particulary if the geology is complex. This hypothesis leads to the consideration that diffusion is anomalous and a reconsideration of Darcy's law; see for example, Nigmatullin (1984Nigmatullin ( , 1986, Le Mẽhautẽ (1984), Dassas and Duby (1995), Henry et al. (2010) to name a few of the many such citations. Anomalous diffusion implies a nonlinear scaling of the mean square displacement of the transient front with time. Well responses often display power-law behavior (Chang and Yortsos, 1990;Raghavan and Chen, 2013a).
The goal of this study is to explore the development and application of interference tests for subdiffusive flows. The first step involves examination of the Theis (1935) solution for subdiffusive flow. Although such a solution has been available for some time, the solution is rarely recognized probably because it was presented as an aside to the study of fractured wells (Raghavan, 2012). We flesh out the details of the solution as they are not yet available and more importantly present a correlation in the manner of Theis (1935) as a ''working'' solution. In addition, we examine new developments that have appeared in the literature and indicate adjustments that need to be made if the goal is to model subdiffusive flow (Su et al., 2015). We then examine pressure distributions by injecting through a vertical fracture along the lines proposed in Uraiet et al. (1977) as the active (pumping) well is often stimulated by a vertical fracture.
The model discussed in this work is a viable option to evaluate well tests in fluvial (and fractured) systems that display near power-law behaviors; indeed, they avoid conjuring up geological options where we would tend to look askance at analyses.

Problem statement
We proceed in the usual way by considering the flow of a slightly compressible fluid of compressibility, c, and viscosity, l, to a well of radius, r w , located at the origin of a reservoir of constant thickness, h, that is at initial pressure, p i , and infinite in its areal extent. The well center is taken to be the origin and the fluid properties are assumed to be constant. Again, as usual, we formulate the problem by combining the conservation equation, the flux law for subdiffusive flow and the equation of state for a slightly compressible fluid. The solution is developed in terms of the Laplace transformation in the cylindrical coordinate system presuming symmetry under the assumptions that the well penetrates the reservoir completely and is produced (pumped) at a constant rate, q. Thus, except for the flux law, a fractional form of Darcy's law, we consider the standard problem. We work in terms of the Laplace transformation and denote the Laplace variable by the symbol, s. As a result, the ideas of Warren and Root (1963) or de Swaan-O (1976) for naturally fractured reservoirs may be incorporated in the formulations that follow through a kernel, F ðsÞ, in the arguments of the Bessel functions; see Noetinger and Estebenet (2000), Noetinger et al. (2001Noetinger et al. ( , 2016, Albinali et al. (2016) and Raghavan and Chen (2017) for illustrations that incorporate subdiffusive mechanisms in the matrix system and in the fracture system if needed, and for the definition ofF ðsÞ. Another direct benefit is the option to include the wellbore storage effect (see Appendix).
Over the past few years, particularly since 2011, a number of articles pertaining to the behavior of wells in subdiffusive systems has appeared in the literature, and different choices have been made with respect to the development of the mathematical model. The choices made in this work are dictated by our experiences in working on shale reservoirs and the need to address reserves and forecast performance. A comprehensive discussion of the choices may be found in Raghavan and Chen (2013a) where we discuss modifications to Beier's (1994) model.

The Darcy law
In our opinion, one of the best expositions of the flux law for subdiffusive flow is the one given in Henry et al. (2010). They show that for subdiffusive flows the velocity, vðx; tÞ, at a location, x, at time, t, is given by where the symbol, pðx; tÞ is the pressure at a point in time, andk The fractional derivative, o a f ðtÞ=ot a , for 0 < a < 1 used in the Caputo (1967) sense is given by where CðÁÞ is the Gamma function. For a ¼ 1, equation (1) reverts to the classical Darcy equation. For reservoir rocks, to our knowledge, only the recent study by Raghavan and Chen (2018) has directly addressed estimates of a; for shales in the Permian basin they report the following estimates: 0:56 a f 0:91 and 0:77 a m 0:94, where the subscripts ''f'' and ''m'' stand for the fissure and matrix systems, respectively. Estimates of a may also be deduced from other works using the guidelines in Beier (1990) and Metzler et al. (1994). Values in the range 0:6 < a < 0:86 are reported in Acuna et al. (1995), while Flamenco-Lopez and Camacho-Velázquez (2001) report that a was in the range 0:47 < a < 0:67. Also from Le Borgne et al. (2004Borgne et al. ( , 2007, we obtain a to be 0:71. The Laplace transformation of the Caputo fractional operator, c 0 D a t c f ðtÞ, is given by This expression in equation (1) implies that the flux at time, t, may not be obtained directly from the instantaneous pressure distribution, pðx; tÞ.

Formulation
Considering cylindrical coordinates and assuming the well axis passes through the origin, the application of the conservation of mass principle to a control volume for a slightly compressible fluid yields 1 r where / is the porosity of the medium, and r is the distance from the well center. On substituting the right-hand side of equation (1) for vðr; tÞ and integrating, we obtain the partial differential equation for transient diffusion under subdiffusive flow to be whereg, the ''diffusivity'' of the medium, is given as We seek a solution to equation (6) subject to the following initial and boundary conditions: Equation (10) specifies the condition for a well of finite radius, r w . The Theis solution, however, requires that we determine the solution when the well radius is vanishingly small; that is, a line source or the solution as lim r w ! 0. As the finite well-radius solution is more general, we consider this general solution first and then extract the solution corresponding to a well that is a line-source.

General solution
The procedure for obtaining the solution we desire is straightforward and may be found in many texts. Working in terms of the pressure difference, Áp ðr; tÞ [p i À p ðr; tÞ, with respect to the initial pressure, p i , and denoting Áp to be the Laplace transform of Áp, the Laplace transform of equation (6) becomes a modified Bessel equation of order 0, the general solution of which is with I 0 ðxÞ and K 0 ðxÞ representing the modified Bessel functions of the first and second kind of order 0, respectively, and If the outer-boundary condition given by equation (9) is to be satisfied, then B ¼ 0, because I 0 ðxÞ ! 1 as x ! 1. Application of the Laplace transform to equation (10) and further simplification provides the constant A; thus, the pressure distribution is 3.1 Solutions for the situation xK 1 (x) as x ! 0 This approximation leads to the consideration of solutions for two important situations. First, if we consider long enough times, then the approximation provides the longtime response for a well of radius, r w ; second, if we were to consider situations where r w ! 0, then the approximation provides the solution for a line-source well, and the analog of the Theis (1935) solution for subdiffusive flow is obtained. These observations also suggest that at long enough times the line-source solution becomes asymptotic to the finite-well-radius solution similar to the classical case; see Mueller and Witherspoon (1965). Because xK 1 ðxÞ ! 1 as x ! 0, the expression for the pressure distribution around a line-source well under subdiffusion turns out to be Saxena et al. (2006) show that where the symbol H m;n p;q ½xj ða 1 ;A 1 Þ;:::;ðap;ApÞ ðb 1 ;B 1 Þ;:::;ðbq ;Bq Þ is the Fox function or the H -function, Rð.Þ > 0, Rðz 2 Þ > 0, RðsÞ > 0. To our knowledge, the Fox function is best computed with the aid of packages such as Mathematica; see Weisstein (2018); see also Mainardi et al. (2005). The simplest alternative to compute Áp r; t ð Þ is by inverting the expression in equation (13) by the Stehfest algorithm (1970a, b).
On defining dimensionless variables corresponding to distance, r D , time, t D , and pressure, p D ðr D ; t D Þ, respectively, by we obtain As noted earlier, an expression identical to that in equation (19) derived by the method of sources and sinks is given in Raghavan (2012) but the details we are about to discuss were not considered, as the goal there was to discuss the performance of fractured wells. First and foremost we should note that for a ¼ 1, equation (19) does reduce to the Theis (1935) solution where ÀEi Àu ð Þ is the Exponential Integral and u ¼ r 2 D =ð4t D Þ. Second, it is important to recognize that the time term in the coefficient of the Fox function arises out of the fractional Darcy law, equation (10), much in the same way as the exponent a on the right hand side of equation (6) arises in our development. There are many studies where the fractional form of Darcy's law is not used (Atangana and Bildik, 2013), and solutions are presented using equation (6) with the classical definition of Darcy's Law; we are yet to understand the physics behind such approaches to study fluid movement in heterogeneous environments for it is not obvious as to how one may arrive at the differential equation of the form used here without the use a fractional form of Darcy's law (or some other similar consideration). Third and most importantly, the expression on the right hand side of equation (19) suggests that a working solution for analyzing responses under subdiffusive flow may be obtained by plotting p D r D ; t D ð Þ=r m D vs t D =r 2 D with m ¼ 2ð1 À aÞ=ðaÞ; this aspect is considered in Section 3.2 on Computational Results. This observation may be advantageous if measurements at several wells are addressed simultaneously during an interference test.

An asymptotic solution of equation (13)
Considering s ! 0, and using the result that for small values of x, K 0 ðxÞ is given by (Carslaw and Jaeger, 1959) where c is Euler's constant, then after ignoring all but the first two terms on the right-hand side of equation (21), we may readily show from that the long-term approximation of equation (14) (or the long-time approximation of equation (13)) is where wðÁÞ is the Digamma function. The logarithmic derivative of p D ðr D ; t D Þ corresponding to equation (22) is, of course, given by Again for a = 1, this result corresponds to the semilogarithmic approximation of the Theis (1935) solution, namely the Cooper and Jacob (1946) approximation. As we discuss below in a number of situations, like that discussed in Thomas et al. (2005), active well responses follow near power-law trends similar to those suggested by equation (23).

The instantaneous, line-source solution corresponding to equation (14)
Using the ideas in Carslaw and Jaeger (1959) it is rather easy to obtain the expression for the pressure distribution, " fðx; yÞ, at a point, ðx; yÞ, caused by an instantaneous, line-source at t ¼ 0 located at the point ðx i ; y i Þ. The expression is The three dimensional version of equation (24) may be found in Mendes et al. (2005). The line-source or pointsource solution is particularly convenient for determining the nonlinear scaling of the space-time behavior of the transient front with the second moment, hr 2 i, of the transient scaling in the form For purposes of understanding well behavior, point-or line-source solutions are particularly convenient for developing pressure distributions in reservoirs produced through complex wellbores; see, for example, Raghavan and Ozkan (1994).

Vertically-fractured wells
The solution for a vertically fractured well may be obtained by suitably integrating the right-hand side of equation (24). Assuming the center of the fracture is at ðx w ; y w Þ and that the fracture extends from ðx w À L f Þ to ðx w þ L f Þ with its plane parallel to y ¼ 0, the pressure distribution in terms of the Laplace transformation is with l being the reference length,qðx; tÞ being the flux distribution, andx wD defined bŷ In the following we consider the uniform-flux solutions; that is,q is independent of x and t. If this were the case, then equation (26) becomes Although we address only the uniform-flux case, equation (26) forms the basis for addressing finiteconductivity solutions; see Cinco-Ley and Meng (1988).
We may evaluate the integral in equation (28) for ðy 6 ¼ y w Þ along the lines indicated in Raghavan and Ozkan (1994). Denoting the integral in equation (28) by I, and the upper and lower limits of I by b and a, respectively, we may write for x D ! b, and for a x D b. The situation for x D a may be obtained by symmetry through equation (29). The conclusions we present in the following apply qualitatively to situations when the well is produced through a finite-conductivity fracture (Meehan et al., 1989;Mousli et al., 1982).

Computational results
The principal goal here is to demonstrate the characteristics of the responses under subdiffusive conditions in addition to demonstrating the solutions are computable and accurate. All results are obtained through the Stehfest algorithm (1970a, b). Figure 1 presents responses at the observation well for the situation where r D ¼ 1000. Both the pressure response and its logarithmic derivative are shown as unbroken and dashed lines respectively for three values of the subdiffusive coefficient, a. In general, the characteristics of the curves are similar to the classical solutions, but the shapes of the derivative curves, in particular, are significantly different at later times because of the existence of the t 1Àa=a D term in equation (22). If this term had not existed, all derivative curves would be much like the classical case. This result is a direct consequence of our use of the fractional form of the Darcy law, equation (1), which accounts for the topological and geometrical influences of the porous rock on hydraulic movement. As we shall see below, this characteristic of the solutions is important. The circles and squares represent the asymptotic expressions given in equations (22) and (23), respectively; agreement with the rigorous solutions is excellent. Figure 2 tests the proposition suggested on the basis of equation (19) that responses under subdiffusive flow may be correlated in terms of p D r D ; t D ð Þ=r m D vs t D =r 2 D where m ¼ 2ð1 À aÞ=a. Results are shown for three values of r D in the range 1 r D 10 3 for an a of 0:8; similar results are obtained for other values of the exponent a. This method of viewing results essentially removes the reference length, r w , and considers responses in terms of the intrinsic properties of the rock. As we shall see, this approach is also advantageous in considering fractured wells. The principal advantage is that it simplifies analysis allowing a working solution which depends solely on the exponent a as a parameter of interest to be prepared. Such a working solution is shown in Figure 3 for analyzing responses influenced by subdiffusion in terms of the exponent a. Derivative responses were correlated in a similar manner and are not shown principally for clarity.
We now turn our attention to the consideration of the influence of a finite wellbore; see Figure 4 where we consider the influence of r D for a ¼ 0:9. The r D ¼ 1 is the well response for the finite-well-radius case. At long enough times, the pressure response follows the trend suggested in equation (22). Furthermore, at distances r D > 20 the results in Figure 4 suggest that the influence of the well radius will be negligibly small, a conclusion similar to that  (22) and (23), respectively. suggested by Mueller and Witherspoon (1965) for the classical case, a ¼ 1. Figure 5 considers pressure distributions around a well produced through a vertical fracture, equation (26). Results are based on assuming the characteristic length to be the fracture half-length, L f ; r D now given by is assumed to be 1:2. The symbol h, the parameter of interest, is the orientation of the line on the horizontal plane (y ¼ 0) joining the center of the fracture ðx w ; y w Þ to the point r ðx À x w ; y À y w Þ with respect to the fracture plane; that is, The results shown here are typical of all values of r D except for the fact that the influence of h diminishes as r D increases and is essentially negligible if r D were greater than 2. This observation, of course, is independent of the subdiffusion coefficient, a. Interestingly, the h ¼ 90 solution corresponds to the line-source solution, equation (19), if responses are plotted in terms of p D r D ; t D ð Þ=r m D vs t D =r 2 D . This is one of the reasons why the line-source solution has been often successfully used to evaluate fractured wells. Considering the alignment of the curves, if determining the orientation, h, is the goal, then choosing a location where h 45 should be preferable. At long enough times, behavior similar to that given in equation (23) will be evident for all orientations, h, and distances, r D . This point may be illustrated by substituting the right-hand side of equation (21) in equation (28) for K 0 ðzÞ and integrating with respect to x; see Raghavan and Ozkan (1994). All of these observations in a qualitative sense also apply to  (19). Both pressure and derivative responses are considered for values of the distance, r D , in the range 1 r D 10 3 .   situations where the fracture conductivity is finite (Meehan et al., 1989;Mousli et al., 1982).

Commentary
Historically, subdiffusive flows in porous rocks have been modeled by considering transient diffusion on a fractal object where the second moment of the distance travelled by the front, hr 2 i, in terms of the anomalous diffusion coefficient or random walk dimension, d w ð> 2Þ, is given by If d w were equal to 2, then classical diffusion results. Equation (33) is based on the presumption of symmetry and that the diffusivity term, g, in the transient diffusivity equation be of the form (Gefen et al., 1983) where h > 0 is given by The stipulation of Gefen et al. (1983) translates to the requirement that both porosity and permeability be dependent on distance, r, namely; see Chang and Yortsos (1990): The symbols d f , and d in equation (36) represent the fractal (Hausdorff) and Euclidean (d = 1, 2, or 3) dimensions respectively. For 2D problems in subsurface rocks, Beier (1994) recommends that we use d f ¼ 2.
In terms of applications to subsurface flow, the requirements of both symmetry and power-law structure specificized in equation (36) impose restrictions if one were to consider issues such as well spacing, fracture spacing (when wells are stimulated by multiple hydraulic fractures), estimation of reserves, location of boundaries, anisotropy, and similar considerations. Beier (1994) recognizes the limitation symmetry imposes on the study of 2D problems as his primary interest was in the application of equation (36) to fractured wells. Nevertheless he proceeds by using equation (36) with the observation that although his solution does ''... not strictly carry over to the vertical-fracture geometry'' and is thus not valid for all times, his solution is exact during the linear and the pseudoradial flow periods. The analog of Beier's solution in terms of the Laplace transformation is discussed in Raghavan and Chen (2013a). In addition, to address the matter of symmetry they propose adopting Ball et al. (1987) and using the expression instead of equation (33). Here d wi is the random walk dimension in the direction, i. The analog of the line-source solution corresponding to equation (14) and reflecting equation (37) with the power-law behavior being preserved. This expression, however, still does not address the power-law dependence of permeability and porosity on distance. Our approach addresses this aspect and still preserves the subdiffusive characteristics we desire for a number of problems of interest; see Raghavan and Chen (2013a, b, 2017, 2018, Albinali et al. (2016); it is similar to that used in Benson et al. (2004) in another context.

Discussion and concluding remarks
The primary contribution of this work is to provide a ''working'' solution to evaluate well responses under subdiffusive flow as it is difficult, if not virtually impossible, to address the appearance of the loss of connectivity that is often evident even in environments other than naturally fractured reservoirs; see, for example, derivative responses in Figure 6 taken from Thomas et al. (2005) at a number of wells producing from a variety of geological environments which display near power-law behavior suggesting the existence of spatial scaling of properties. For the responses shown, based on the model discussed here, the slopes of these lines suggest that 0:6 a 0:8. The responses are for wells producing distinct fluvial deposits in a number of geographical locations in an exploration context. Consequently, it becomes untenable to tap into fractured-well models if we are circumscribed by the requirement that the classical option be used. We may consider the option that the wells produce commingled deposits where the lateral extent of some of the zones is limited in that these zones deplete during the producing phase and are fed by the more extensive zones during the buildup period (differential depletion, backflow); see Gao et al. (1994). Again, it seems that this option would not be viable in each of these tests when the geological setting does not suggest that this is the case. Models of the kind discussed here enable us to eschew the option of ''simple physics and simple geology'' when indeed the responses suggest that the geology is complex. The concept of subdiffusion offers a viable option to evaluate well tests which display near power-law behaviors; as already mentioned, classical methods to address tests of the kind shown in Figure 6 require the conjuring up of counterfactual geological options or ignoring heterogeneity altogether. Models incorporating subdiffusive effects have been shown to be a promising alternative in explaining the performance of shale reservoirs; see Chu et al. (2017Chu et al. ( , 2018.