Partially Penetrated Well Solution of Fractal Single-Porosity Naturally Fractured Reservoirs

In the oil industry, many reservoirs produce from partially penetrated wells, either to postpone the arrival of undesirable fluids or to avoid problems during drilling operations. The majority of these reservoirs are heterogeneous and anisotropic, such as naturally fractured reservoirs. The analysis of pressure-transient tests is a very useful method to dynamically characterize both the heterogeneity and anisotropy existing in the reservoir. In this paper, a new analytical solution for a partially penetrated well based on a fractal approach to capture the distribution and connectivity of the fracture network is presented. This solution represents the complexity of the flow lines better than the traditional Euclidean flow models for single-porosity fractured reservoirs, i.e., for a tight matrix. The proposed solution takes into consideration the variations in fracture density throughout the reservoir, which have a direct influence on the porosity, permeability, and the size distribution of the matrix blocks as a result of the fracturing process. This solution generalizes previous solutions to model the pressure-transient behavior of partially penetrated wells as proposed in the technical literature for the classical Euclidean formulation, which considers a uniform distribution of fractures that are fully connected. Several synthetic cases obtained with the proposed solution are shown to illustrate the influence of different variables, including fractal parameters.


Introduction
In the literature, several analytical solutions for modeling the behavior of pressure-transient tests of partially penetrated wells have been proposed [1][2][3][4][5][6][7][8][9][10].Some of these works have proposed the use of point and line source solutions derived in the Laplace space, considering finite and infinite systems, with homogeneous and naturally fractured reservoirs [2,5,6,8,10].Other studies considered gas anisotropic reservoirs using a uniform flow solution [9].All of these works assumed reservoirs with Euclidean geometry, that is, they used traditional mass conservation and flow equations.
Starting from mass conservation and flow equations with fractal characteristics, the authors of [11][12][13][14] analyzed the behavior of the pressure-transient tests of single and double porosity reservoirs with fractal geometry.These studies established the existence of a power-law behavior during the transient period instead of the classical semi-logarithmic behavior that exists in reservoirs with Euclidean geometry.It has been demonstrated that the radial flow regime is a special case of more general fractal behavior.All these studies considered vertical fully penetrated wells.Up to date, no study has been presented that considers the pressure-transient behavior of partially penetrated wells produced from anisotropic heterogeneous reservoirs with fractal properties.
In this study, a single-porosity system was considered, which can be represented by a naturally fractured reservoir with a tight matrix, where the porosity and permeability of the system are due to the fracture network.Additionally, it was considered that there was a folding where the density of fractures was greater at the top of the anticline and decreased toward the flanks.Thus, there was a heterogeneous and anisotropic reservoir where the radial and vertical permeabilities were functions of the radial and vertical position, respectively.Due to the complexity of this fracture network, it was convenient to consider fractal geometry, instead of assuming a uniform distribution of fractures, and all fractures as being interconnected, as is considered in the traditional formulation with Euclidean geometry.
The purpose of this work was to obtain an analytical solution that represented the behavior of pressure-transient tests in vertical wells partially penetrating heterogeneous and anisotropic reservoirs with fractal geometry.The heterogeneity and anisotropy were due to a fracture network caused by the thrust of a salt dome.

Problem Statement
The solution proposed in this study considered a closed cylindrical reservoir with a single porosity, i.e., a network of fractures may exist, but the matrix is compact and does not contribute to the reservoir response.The well was produced from a restricted interval of the formation.In the reservoir, there were fractal distributions of permeability and porosity in the radial and vertical directions, that is, it was a heterogeneous and anisotropic reservoir.Using the continuity equation in cylindrical coordinates: considering a distribution of permeability in the fracture network like that existing in an anticline, where the radial permeability decreases as the radial distance from the center of the anticline increases, and the vertical permeability also decreases with the increment of vertical depth from the top of the anticline.Thus, the fractal distribution of permeability in the radial and vertical directions are given as follows: where k rw and k zw represent the radial permeability at the center and the vertical permeability at the top of the anticline, respectively.D r = 2 and D z = 1 are the Euclidean dimensions in the horizontal and vertical directions, respectively.The fracture density is represented by the fractal dimensions d f r and d f z , in the radial and vertical directions, respectively.θ r and θ z represent the connectivity indexes of the fracture network in the radial and vertical directions, respectively.The definition of radial permeability is similar to that used in References [11][12][13][14].
The porosity of the fracture network is also a function of the radial distance from the center of the anticline and the vertical position from the top of the anticline.Thus, using the fractal definition of porosity proposed by Cossio et al. [15] in 2D (r and z), the fracture porosity is given by: where φ 0 represents the average porosity in the near wellbore region at the top of the reservoir.In the following, we use φ 0 = φ 0 /2.
Assuming Darcy´s Law for the velocities in the radial and vertical directions and considering Equations ( 2)-(4) into Equation (1), the following equation can be obtained: It can be noted that instead of fractal derivatives, fractal definitions of the petrophysical properties are used in the derivation of this equation, following a similar path to that proposed in References [11][12][13][14][15][16][17].Some applications of the use of fractional derivatives on the fluid flow in porous media are presented elsewhere [18][19][20][21].Using the values of the Euclidean dimensions in the horizontal and vertical directions, D r = 2 and D z = 1, we obtained the following equation: Using the following definitions of dimensionless variables: Considering a slightly compressible fluid of constant viscosity (µ), and small pressure gradients, we obtained: In Figure 1, a diagram of the problem to be solved in cylindrical coordinates is shown.Applying Newman's method according to Razminia et al. [10], "the instantaneous Green function is equal to the product of the instantaneous Green functions in one and/or two directions", in our case: With the above, Equation ( 19) will be solved for the two directions independently.Applying Newman's method according to Razminia et al. [10], "the instantaneous Green function is equal to the product of the instantaneous Green functions in one and/or two directions", in our case: With the above, Equation ( 19) will be solved for the two directions independently.

Analytical Solution of the Problem
The solution was deduced by applying the methods of the Laplace's transform, separation of variables, and Newman's product using instantaneous source functions.In Appendix A, the procedure for obtaining the solution in the radial direction for total penetration, Equation (A.7), can be found.This solution was used together with the solution in the vertical direction, Equation (B.18), obtained in Appendix B, to acquire the solution for a partially penetrated well through the use of the Newman's product.Thus, Equation (B.26) is written as follows: where: If this expression is evaluated for ℎ = 1,  = 0, and ℎ = 1 , we obtain the fully-penetrated well solution:

Analytical Solution of the Problem
The solution was deduced by applying the methods of the Laplace's transform, separation of variables, and Newman's product using instantaneous source functions.In Appendix A, the procedure for obtaining the solution in the radial direction for total penetration, Equation (A7), can be found.This solution was used together with the solution in the vertical direction, Equation (A25), obtained in Appendix B, to acquire the solution for a partially penetrated well through the use of the Newman's product.Thus, Equation (A33) is written as follows: where: If this expression is evaluated for h pD = 1, z wD = 0, and h wD = 1, we obtain the fully-penetrated well solution: where λ n are the characteristic values given by the roots of Equation (A17), and: The second term of Equation ( 23) represents the pseudo-skin due to partial penetration considering fractal behavior in both radial and vertical directions.
To include wellbore storage and mechanical skin effects, the following expression, given by Van Everdingen and Hurst [22], is applied: where p D (s) is given by Equation ( 23).

Results
In this section, some results are presented with the proposed analytical solution given by Equations ( 23) and (25) in the case of wellbore storage and skin effects using Stehfest's algorithm [23].
Figures 2-5 show the solution for a Euclidean isotropic case (d fr = 2.0, θ r = 0, d fz = 1.0, θ z = 0), where the upper part of the formation is open to production. Figure 2 shows results without mechanical skin damage, S = 0, where only the thickness of the formation varies.The dashed lines in Figures 2-5 correspond to the pressure and pressure derivative given by Razminia et al. [10] for some Euclidian cases.In all cases, the agreement is excellent, so the proposed solution, Equation (23), is able to reproduce the Euclidian results as particular cases.
considering fractal behavior in both radial and vertical directions.
To include wellbore storage and mechanical skin effects, the following expression, given by Van Everdingen and Hurst [22], is applied: where ̅ () is given by Equation ( 23).

Results
In this section, some results are presented with the proposed analytical solution given by Equations ( 23) and (25) in the case of wellbore storage and skin effects using Stehfest's algorithm [23].
Figures 2-5 show the solution for a Euclidean isotropic case (dfr = 2.0, θr = 0, dfz = 1.0, θz = 0), where the upper part of the formation is open to production. Figure 2 shows results without mechanical skin damage,  = 0, where only the thickness of the formation varies.The dashed lines in Figures 2-5 correspond to the pressure and pressure derivative given by Razminia et al. [10] for some Euclidian cases.In all cases, the agreement is excellent, so the proposed solution, Equation (23), is able to reproduce the Euclidian results as particular cases.In Figure 3, the magnitude of the open interval varies, including the case of the fully penetrated well, keeping the thickness of the formation constant.In Figures 4 and 5, the mechanical skin damage and wellbore storage vary, respectively, keeping the thickness of the formation and the open interval constant.All these cases are Euclidean and serve to evaluate the accuracy of the fractal analytic solution proposed for these cases.
Fractal Fract.2019, 3, 23 7 of 17 The cases with fractal geometry are shown below.Figure 6 shows a case where the fractal dimension in the radial direction is varying, dfr ≤ 2, where the value of 2 represents the Euclidean case (θr = 0).Thus, the traditional Euclidean case is a special case of the fractal case.In the Euclidean case, the classical spherical flow with a slope of −0.5, before the radial period, is present.It can be observed that this period of flow is not present for the fractal cases, where instead of the semi-logarithmic period, a power-law behavior can be observed in both the pressure drop and its derivative at late times during the transient period.
In Figure 3, the magnitude of the open interval varies, including the case of the fully penetrated well, keeping the thickness of the formation constant.In Figures 4 and 5, the mechanical skin damage and wellbore storage vary, respectively, keeping the thickness of the formation and the open interval constant.All these cases are Euclidean and serve to evaluate the accuracy of the fractal analytic solution proposed for these cases.
The cases with fractal geometry are shown below.Figure 6 shows a case where the fractal dimension in the radial direction is varying, d fr ≤ 2, where the value of 2 represents the Euclidean case (θ r = 0).Thus, the traditional Euclidean case is a special case of the fractal case.In the Euclidean case, the classical spherical flow with a slope of −0.5, before the radial period, is present.It can be observed that this period of flow is not present for the fractal cases, where instead of the semi-logarithmic period, a power-law behavior can be observed in both the pressure drop and its derivative at late times during the transient period.The cases with fractal geometry are shown below.Figure 6 shows a case where the fractal dimension in the radial direction is varying, dfr ≤ 2, where the value of 2 represents the Euclidean case (θr = 0).Thus, the traditional Euclidean case is a special case of the fractal case.In the Euclidean case, the classical spherical flow with a slope of −0.5, before the radial period, is present.It can be observed that this period of flow is not present for the fractal cases, where instead of the semi-logarithmic period, a power-law behavior can be observed in both the pressure drop and its derivative at late times during the transient period.
Figure 7 shows fractal cases where the connectivity index in the radial direction (θ r ) varies, and now all other parameters are kept constant, including the fractal dimension d fr = 2. Again, the Euclidean case occurs when θ r = 0, i.e., radial flow exists at late times during the transient period, and when θ r > 0, the power-law response is present at these times.
Fractal Fract.2019, 3, 23 8 of 18 Figure 7 shows fractal cases where the connectivity index in the radial direction (θr) varies, and now all other parameters are kept constant, including the fractal dimension dfr = 2. Again, the Euclidean case occurs when θr = 0, i.e., radial flow exists at late times during the transient period, and when θr > 0, the power-law response is present at these times.In Figures 8 and 9 the fractal dimension, dfz, and the connectivity index, θz, are varied in the vertical direction, respectively, keeping the other parameters constant, including dfr = 2, and θr = 0.The influence of dfz and θz is observed only in the period before the radial flow.In these cases, when dfz = 1.0 and θz = 0, the traditional Euclidean case is obtained again, with the presence of spherical flow before the radial period.
In Figures 8 and 9 the fractal dimension, d fz , and the connectivity index, θ z , are varied in the vertical direction, respectively, keeping the other parameters constant, including d fr = 2, and θ r = 0.The influence of d fz and θ z is observed only in the period before the radial flow.In these cases, when d fz = 1.0 and θ z = 0, the traditional Euclidean case is obtained again, with the presence of spherical flow before the radial period.In Figures 8 and 9 the fractal dimension, dfz, and the connectivity index, θz, are varied in the vertical direction, respectively, keeping the other parameters constant, including dfr = 2, and θr = 0.The influence of dfz and θz is observed only in the period before the radial flow.In these cases, when dfz = 1.0 and θz = 0, the traditional Euclidean case is obtained again, with the presence of spherical flow before the radial period.In Figures 10 and 11, the influence of hpD and the mechanical skin is shown, respectively, keeping the other parameters constant, including the fractal parameters.At large times within the transient period, the power-law behavior can be detected.In fact, in Figure 11, the presence of two power-law periods is observed.In Figures 10 and 11, the influence of h pD and the mechanical skin is shown, respectively, keeping the other parameters constant, including the fractal parameters.At large times within the transient Fractal Fract.2019, 3, 23 9 of 17 period, the power-law behavior can be detected.In fact, in Figure 11, the presence of two power-law periods is observed.In Figures 10 and 11, the influence of hpD and the mechanical skin is shown, respectively, keeping the other parameters constant, including the fractal parameters.At large times within the transient period, the power-law behavior can be detected.In fact, in Figure 11, the presence of two power-law periods is observed.

Figures 12 and 13
show the influence of the fractal parameters in the vertical direction, considering a fractal condition in the radial direction.In Figure 12, it is observed that the effect of the fractal dimension, dfz, is not very strong; however, it can be expected that with the arrival of undesirable fluids to the producing well, this parameter could play an important role.In both figures, the presence of two power-law periods is observed.Figure 13 shows that when the connectivity of fractures or pores in the vertical direction decreases, or even becomes null (i.e., θz = 1), the late power-law period is delayed, which is an expected behavior.
Figures 12 and 13 show the influence of the fractal parameters in the vertical direction, considering a fractal condition in the radial direction.In Figure 12, it is observed that the effect of the fractal dimension, d fz , is not very strong; however, it can be expected that with the arrival of undesirable fluids to the producing well, this parameter could play an important role.In both figures, the presence of two power-law periods is observed.Figure 13 shows that when the connectivity of fractures or pores in the vertical direction decreases, or even becomes null (i.e., θ z = 1), the late power-law period is delayed, which is an expected behavior. = 1.0,  = 0.2.
Figures 12 and 13 show the influence of the fractal parameters in the vertical direction, considering a fractal condition in the radial direction.In Figure 12, it is observed that the effect of the fractal dimension, dfz, is not very strong; however, it can be expected that with the arrival of undesirable fluids to the producing well, this parameter could play an important role.In both figures, the presence of two power-law periods is observed.Figure 13 shows that when the connectivity of fractures or pores in the vertical direction decreases, or even becomes null (i.e., θz = 1), the late power-law period is delayed, which is an expected behavior.Considering the above results, it can be deduced that the new proposed analytical solution may provide useful information for the proper development of a reservoir.However, it can be intuited that to determine all the parameters involved in the proposed analytical solution, it is necessary to use a robust optimizer, since a visual adjustment is expected to be very difficult to apply for a complex model such as the one proposed in this work.

Discussion
Taking into account the above results, and those presented by Posadas and Camacho [14], and the fact that there are many unknown parameters (S, CD, ,  ,  ,  ,  , kr) to fully characterize Considering the above results, it can be deduced that the new proposed analytical solution may provide useful information for the proper development of a reservoir.However, it can be intuited that to determine all the parameters involved in the proposed analytical solution, it is necessary to use a robust optimizer, since a visual adjustment is expected to be very difficult to apply for a complex model such as the one proposed in this work.

Discussion
Taking into account the above results, and those presented by Posadas and Camacho [14], and the fact that there are many unknown parameters (S, C D , ε, d f r , θ r , d f z , θ z , k r ) to fully characterize this system, it is necessary to use robust optimization software in the type-curve matching process of both the pressure and its semi-logarithmic pressure derivative in order to obtain all of these parameters from well test data.

Conclusions
The novel analytical solution presented in this paper considers for the first time the application of fractal geometry to the problem of partial penetration.This is relevant because it allows the consideration of the variation of petrophysical properties with the scale or it takes into account the tortuosity of the flow lines in a cylindrical system.The solution was deduced by applying the methods of the Laplace's transform, separation of variables, and Newman's product using instantaneous source functions.Considering the results presented in this article, we can conclude the following: 1.
The new fractal analytical solution for a constant rate describes the pressure-transient behavior for partially penetrating wells in a single-porosity naturally fractured reservoir and includes the traditional Euclidean solution as a special case.

2.
The proposed fractal solution generates a power-law response at late times during the transient period after the wellbore storage, mechanical skin, and partial penetration effects have ended.This behavior occurs when the radial fractal parameters are different from the Euclidean values, i.e., d fr < 2 and θ r > 0.

3.
A different behavior to the power-law response occurs when d fz < 1 and θ z > 0. The effect of these parameters is shown only during the partial penetration period, and after this period, the traditional radial behavior (if d fr = 2 and θ r = 0) or a power-law behavior (when d fr < 2 and/or θ r > 0) can be present.

4.
The typical spherical flow regime due to partial penetration is only present when the fractal parameters in the radial direction have the Euclidean values, i.e., d fr = 2 and θ r = 0. 5.
An expression is provided to evaluate the pseudo-skin due to the partial penetration effects that consider fractal behavior in both the radial and vertical directions.6.
To determine the pseudo-damage due to restricted penetration, horizontal permeability, vertical to horizontal permeability ratio, mechanical skin, and the four fractal parameters, it is necessary to resort to a type-curve matching process of the pressure data and its semi-logarithmic derivative using a robust optimizer that minimizes the difference between the real data and the analytical solution.Connectivity index in the radial direction (0 ≤ θ r ≤ 1) θ z Connectivity index in the vertical direction (0 ≤ θ z ≤ 1) ξ(s) Pseudo-skin due to partial penetration considering fractal behavior

Appendix A. Solution in the Radial Direction
The flow in the radial direction is obtained from Equation (19) as follows: where Using the Laplace transform: Applying the Levedev [24] technique, we are able to obtain the solution in terms of the modified Bessel functions: where: Using the following boundary conditions in the radial direction: It is found that the solution in the radial direction for total penetration in the Laplace space is given by:

Appendix B. Solution in the Vertical Direction
From Equation ( 19), the continuity equation in the vertical direction is given by: z where Applying separation of variables in Equation (A8): Thus, the solution for u(t D ) is given by: and the solution for w(z D ) is obtained from: Applying the Levedev [24] technique, we are able to obtain the solution for w (z D ), in terms of the Bessel functions: where: Substituting Equations (A10) and (A12) into Equation (A9), the general solution for the problem in the vertical direction is given by: Considering the following boundary conditions in the vertical direction: We obtain from Equation (A15), D = 0. From Equation (A16) we obtain: From the roots of Equation (A17) the characteristic values λ are obtained.
Applying the superposition principle with Equation (A17) we obtain the following expression: Considering an instantaneous source plate with its center in z wmD , which agrees with the midpoint of the producing interval, the following expression is obtained: Multiplying Equation (A19) by x = z 2+θz 2

D
and then applying the orthogonality property, we obtain: To evaluate the term of the integral on the right-hand side of Equation (A20), Abramowitz and Stegun [25] is used, obtaining the following expression: where . To evaluate the term of the integral on the left-hand side of Equation (A20), we use Gradshteyn and Ryzhik [26], which is expressed as: Obtaining the following: Substituting Equations (A21) and (A23) into Equation (A20), we obtain: Substituting Equation (A24) into Equation (A18), the following solution is obtained: According to Razminia et al. [10] the instantaneous source function for a partial penetration is defined as a function of the instantaneous source function for total penetration, such as: Using the method of Newman's product, the instantaneous source function for partial penetration can be obtained as follows: S(r D , z D , t D ) = S r (r D , t D ) • S z (z D , t D ). the solution in the Laplace space is given by: Substituting Equation (A7) into Equation (A30), the final solution is given as follows: √ s+λn]+Kv r +1 [( 2 2+θr ) √ s+λn] a n J vz −1 (an)− vz an J vz (an) (A32) Finally, the following expression is obtained for the wellbore pressure drop: If Equation (A33) is evaluated for h pD = 1, z wD = 0 and h wD = 1, we obtain the fully penetrated well solution, given by: Thus, the pseudo-skin due to partial penetration considering fractal behavior is given by: (A36)

2 .
Equations (A25) and (A26) into Equation (A27):S(r D , z D , t D ) = 2 h pD ∂p D f (r D , t D ) v z −1 (a n ) − v z a n J v z (a n ) (r D , z D , t D ) = t D 0 S(r D , z D , τ)dτ,(A29) λ n are the roots of Equation (A17).The evaluation at the wellbore is obtained by evaluating Equation (A31) for r D = 1 in the producing interval, using: (r D = 1, z D , s)dz D .