On the numerical simulation of wellbore pressure in gas reservoirs incorporating the phenomena of slippage, formation damage and wellbore storage

— In this work, we consider the effects of gas slippage and wellbore storage in shale gas reservoirs. We use the Finite Difference Method for the discretization of the nonlinear governing equation, and the iterative Gauss-Seidel method is applied to obtain the solution of the algebraic system. We also perform pressure tests in the producing well through the use of numerical simulation using cylindrical coordinates. The results, obtained in the well testing analysis context, show the relevance of the introduction of the slip, formation damage and wellbore storage on the ﬂow simulation in shale gas reservoirs.


I. INTRODUCTION
Archaeological records show that ancient civilizations in different parts of the world used oil. Incas and Babylonians used the oil to lay bricks, to pave roads, and to waterproof ceramic artifacts [6,34]. However, it was only in the 20th century that oil began to stand out as economic activity on the world stage, with the emergence of the large automotive industries and oil companies, with the latter beginning the search for more oil reserves. On the other hand, natural gas has been in the background for a long time, mainly due to its specificities concerning transport and storage. However, this situation has changed, mostly in the transition between the 20th and 21st centuries.
The market for natural gas of fossil origin has grown worldwide for several reasons, such as the discovery of new deposits, the lower degree of pollution that its burning causes, when compared to oil, and its diversity of applications. The economic relevance of natural gas was evident in the years 1973 and 1979 when the Organization of Petroleum Exporting Countries (OPEC) raised oil prices to a level that led to changes in hydrocarbon exploration and production and the general dynamics of generation and con-sumption of energy in the world. As a consequence, activities associated with natural gas increased [4]. In this context, Fig. 1 presents the Brazilian energy matrix in the 1970s and 2010. In this figure, we can see that the consumption of natural gas increases from 3 to 9%.
In Brazil, in 1941, natural gas began to appear in the economic scenario due to the discovery of the first commercial oil field in Candeias, Bahia, the industrial sector being its first destination. Due to the adverse effects of the two oil shocks, there was an increase in the production of natural gas due to the exploration and development of the Campos Basin in Rio de Janeiro [9]. At the end of the 1990s, with the completion of the construction of the Bolivia-Brazil gas pipeline (GASBOL) [7], the import of natural gas from Bolivia began, and it consolidated itself as a representative energy source in the Brazilian energy matrix, leaving aside the oil by-product label.
Currently, with the discoveries of pre-salt reserves in the Santos Basin [32], interest in natural gas reserves in Brazil has been growing even more rapidly in recent years. Figure 2 shows the growth of gas production over the years and the dependence on Bo-livian gas imports, which is necessary to supply national consumption. However, the Energy Research Company (EPE) disclosed, through its National Energy Plan [14], that the national natural gas supply will be adequate to meet the demand of all the commitments assumed until 2030, thus making self-sufficient Brazil in the sector.  Fig. 1: Brazilian energy matrix in the 1970s and 2010 [15].

Obtaining natural gas
We define a reservoir of natural gas as a deposit that presents a mixture of hydrocarbons in the gaseous state [17]. As stated by Ezekwe [17], we can classify gas reservoirs as a reservoir of dry gas, wet gas, or retrograde condensation gas, according to the type and behavior of the fluid. When there are no molecules heavy enough to form liquid hydrocarbons after surface separation processes, we say that the reservoir is to be dry gas. We call the deposit a wet gas reservoir when it produces a certain amount of liquid in the well and/or in the surface facilities. In turn, as we produce gas from a retrograde condensation gas reservoir, the pressure decreases, and the temperature remains approximately constant. Thus, there is a liquid phase formation in the reservoir through the condensation of the gas.
Another classification that can be applied is that of conventional and unconventional reservoirs. In conventional gas reservoirs, it is possible to produce from vertical wells while it is required specific production techniques for unconventional reservoirs. For example, we can cite the horizontal wells and hydraulic fracturing since there is higher resistance to the gas flow due to the lower absolute permeability of the reservoir. However, with the advancement of technological development, the gas of unconventional origin began to gain prominence in the world economic scenario. India [21], Brazil [27], and USA [40]  nations with the largest reservoirs of shale gas in the world, a type of unconventional reservoir. The trend is that gas of unconventional origin will be essential in the world economy in the coming decades. Unconventional gas reservoirs can be of different types: deep gas reservoirs (located beyond 4,500 meters in depth); low permeability (tight gas); shale gas; gas adhered to coal veins (coalbed methane); gas from pressurized zones (at very high pressure when compared to other reservoirs with the same depth); and underwater and Arctic hydrates (methane hydrates) [2]. The focus of this work, however, is the shale gas reservoirs.

Shale gas reservoirs
The shale gas reservoir sedimentary rock has a fine granulometry and very low permeability [22]. Also, the shales have quite variable mineralogical formation, with a predominance of brittle minerals such as quartz, carbonates, and feldspars [20]. Although it has characteristics of unconventional reservoirs, its exploitation has increased over the years.
Hasan et al. [21] has shown that two-thirds of the world's hydrocarbon reserves are unconventional and that this fact is directly related to the growing importance of shale gas in the world energy matrix. On the other hand, Gomes [19] claims that technological advances in the sector are allowing its production to become attractively economical. It has been drilled more than 50,000 wells over two years in gas reservoirs in the United States [28], making the country's hydrocarbon production almost double.
Countries such as the USA and Brazil have a large number of gas reservoirs, which is very important in economic terms and of influence in the international socio-political framework. As a result of the growth in non-conventional gas exploration, these countries' dependence on the world's largest oil and gas producers, such as Venezuela and Russia, can end.
Despite the environmental concern with fossil fuels and the rapid growth in the use of fuels from renewable sources, Jia et al. [22] believe that fossil energy should still account for 78% of global energy consumption in the year 2040. Even though natural gas of fossil origin is not renewable, technological advances in the area and the growth of discoveries of shale gas reservoirs guarantee its use for still many years. Knudsen et al. [24] show that in recent years the United States has been reducing the use of coalfired power plants and replacing them with natural gas in the generation of electricity. In part, we can explain this by the desire to reduce emissions of polluting carbon dioxide, which harms the environment so much.
Wellbore Testing Analysis studies the pressure and flow changes as a function of production time, through measurements at the bottom of the well and the flow at the surface. From the measured pressure response, it is possible to determine the reservoir properties useful for production planning [41]. In a well test, a transient pressure response occurs due to a production/injection flow. Depending on the objectives of the test, we record the response of the well over a relatively short period when compared to the productive life of the reservoir. In this work, in addition to a sensitivity analysis, well pressure tests using numerical simulation is also done.

II. POROUS MEDIA GAS FLOW
Petroleum is a mixture of hydrocarbons, which has its physical state-properties determined, in general, by its composition, temperature (T ), and pressure (p). According to Ezekwe [17], oil is the part that remains in the liquid state when a mixture of hydrocarbons is brought from the reservoir conditions to the surface conditions, while the natural gas presents the gaseous state in the surface conditions. Under reservoir conditions, natural gas can be present in gaseous form or dissolved in oil.
Regarding the composition of natural gas, the amount of each component can vary depending on the type of reservoir and its characteristics, for example, the location (land or sea), the type of soil, and the geological formation process of the basin, among other factors [30]. However, as can be seen in Table 1 which presents the typical composition of a natural gas reservoir, it is evident that the primary component is methane, which may represent an amount of 70 to 98% of the total natural gas, and in smaller quantities, considered as impurities, carbon dioxide, hydrogen sulfide, and nitrogen. Besides temperature and pressure conditions, another fundamental parameter for calculating the properties of the gas is its relative density (or specific gravity), γ, which is the ratio between the molecular mass of the gas, M , and the molecular mass of air, M air . In this work, it is considered a reservoir of dry gas, produced without the appearance of liquid at any time of production. up to 0.5% For gas mixtures of hydrocarbons, we use the pseudo-critical pressure and temperature coordinates, p pc and T pc , respectively, to determine the socalled pseudo-reduced coordinates [26], and we use them to calculate the physical properties of natural gas.
Sutton [38] presents, depending on the gas density, the correlations that we apply here to obtain the pseudo-critical pressure and temperature, which are fundamental for the reservoir simulations [16]. We use them in determining, for example, the compressibility factor (Z), volume formation factor (B) and viscosity (µ) [37]. From Z and the universal gas constant, R, it is possible to determine the density, ρ, from the equation of state for a real gas, ρ = pM/ZRT .
The gas volume formation factor is the relationship between the volumes it occupies under reservoir conditions (V ) and standard conditions (V sc ) (pressure, p sc , and temperature, T sc , in standard conditions) [16].
On the other hand, we calculate the viscosity of natural gas using the correlation suggested by Lee et al. [25], widely used in reservoir simulation.
Here, the effective porosity (φ) varies depending on the pressure [16]: where φ 0 and p 0 are, respectively, the porosity and pressure in the reference conditions. c φ is the coefficient of compressibility of the rock, and we assume that the compressibility of the rock is small and constant.
In addition to porosity, the economic viability of a reservoir also depends on the permeability of the rock. This property is a measure of a porous material's ability to allow fluids to pass through its pores. We usually represent the absolute permeability by the tensor k.

Slip in gas flow
Studies and predictions about the flow of gas in porous media are more difficult to carry out than those of liquid because the gas properties generally depend more strongly on pressure and also due to the different mass transport mechanisms that can be present [28,31]. Therefore, in some cases, the classic Darcy's law does not adequately describe the flow physics, and experimental data suggest corrections for the calculation of permeability and, thus, we introduce a modified Darcy's law [28] where k a is the apparent permeability tensor, v is the surface velocity of the fluid, g is the acceleration of gravity and D is the depth. Specifically for gases, the slip flow regime occurs when the average free path of the gas molecules has a scale comparable to the pore size [18]. So, both the reservoir and fluid properties influence the determination of apparent permeability. We can mention, among the non-Darcy effects, that we can incorporate in the apparent permeability, high flow rates (inertial and turbulent effects), non-Newtonian fluid flow (for liquids), and slip flow, which occurs only for gas under certain reservoir conditions [5].
When the fluid is a gas, the Klinkenberg effect shows that the permeability measurements made in the laboratory result in values higher than the absolute values, due to the slipping of the gas on the walls of the porous medium. This slip results in a higher flow and leads to a correction of the apparent permeability [23], where b is the Klinkenberg parameter and k the absolute permeability tensor. In reality, the mass is transported in the porous medium by a variety of mechanisms, one of which is the so-called Knudsen diffusion. The Knudsen number measures the relationship between the mean free path of the molecules, λ, and the character- 2τ k/φ, and τ is the tortuosity of the porous medium.
The slip flow regime occurs for 10 −3 < Kn < 0.1 and we can also introduce a different model for determining apparent permeability [28] k a = 1 + 4Kn

Governing equation
In obtaining the partial differential equation (PDE) that governs the isothermal flow of a gas in a porous medium, we employ the mass conservation equation and the modified Darcy's law. We also take into account the effects of slippage, wellbore storage and formation damage, disregarding the phenomenon of gas adsorption, the gravitational force and non-Darcy behaviors related to inertial effects.
As it takes time to the hydrocarbons in the reservoir to reach the surface, in the first moments, we produce the fluid initially stored in the well. This effect is called wellbore storage. The formation damage concerns the reduction of permeability in the region close to the well, caused by wellbore fluids used during drilling and completion.
According to Li et al. [29], mass conservation for the flow of gas in porous media can be described, excluding the adsorption effects and source terms [3], by ∂ ∂t Then, replacing Eq. (2) in Eq. (5) and ignoring the effects of gravity due to the low specific gravity of the gas and the thickness of the reservoir, We can rewrite the term ∂(φ/B)/∂t if we take into account the fluid and rock properties [13] ∂ ∂t where we also employed Eq. (1).
To study the flow dynamics in the region close to the producing well, we assume a two-dimensional flow in cylindrical geometry in the rz-plan and a diagonal permeability tensor, so that: where k ar and k az are the apparent permeabilities in the rand zdirections, respectively. As we are considering the effects of storage in the well, the production flow in it is given by [29] Q sc = q sc + C sc dp wf dt , where q sc is the flow from the porous medium, C sc is the storage coefficient (that already incorporates B) and p wf the pressure in the well. The flow rate q sc is calculated by [35] where J w is the productivity index. Finally, as an initial condition we impose where p ini represents the initial pressure before the reservoir undergoes any changes due to fluid production/injection. On the other hand, the external boundary conditions are of null flow at the external borders, where L z is the thickness of the reservoir and r e is the outer radius of the reservoir. In turn, for the internal boundary condition, where h is the thickness of the production region considered and r w is the radius of the well.

III. NUMERICAL RESOLUTION METHODOLOGY
We employ a computational mesh of centered blocks [1,5,10,16] and the cylindrical coordinate system (r and z), Fig. 3, and we also assume the angular symmetry of the flow. We obtain the numerical solution in the nodes of the computational mesh, located in the centers of the cells, with n r and n z being the numbers of cells in the rand zdirections, respectively. We also use fractional indexes i ± 1/2 and k ± 1/2 to indicate cell interfaces of the computational mesh. So, for the cell i, k and the time n + 1 we can write that where the governing equation was multiplied by the volume of the cell, V i,k , and T stands for the transmissibility.
Following the techniques traditionally applied in reservoir simulation [5,16] and employing centered differences and Similarly, for the discretization of pressure derivatives in the r-direction, and ∂p ∂r where ∆r i±1/2,k is the distance between nodes of cells i, k and i ± 1, k. We can obtain approximations for the derivatives in the z-direction analogously.
In the case of the accumulation term, we employ a conservative expansion [16] Γ n+1 where ∆t is the time step. Next, we obtain a totally implicit formulation for Eq. (14) using a backward Euler approximation, 3.1. Grid refinement For an accurate determination of the pressure gradient, we use to construct the non-uniform mesh in the r-direction (where we suppress k for simplicity of notation) [16], such as: 1. we space pressure calculation points using where i = 1, 2, ..., n r − 1; 2. we define cell boundaries through where i = 1, 2, ..., n r − 1; and 3. we calculate cell volumes employing In the case of the z-direction, we use ∆z = L z /n z . We now introduce the definitions of transmissibilities: and where, considering that the total angle is 2π and the apparent permeabilities [16], and and we apply arithmetic mean to calculate the fluid properties [16].

Numerical approximation for storage
When we consider the wellbore storage effect, the total flow rate is determined as proposed by Li et al. [29] and Tavares [39]: with where 3.3. The solution of the system of equations After the process of discretizing the governing equation, as the resulting algebraic equations are non-linear, we must apply a linearization technique, and we chose the method of Picard [33]. So this results in where the index v refers to the previous iterative level when obtaining the pressure. As we can see, the properties and coefficients of the discretized equation are determined at the iterative level v and used after to calculate the new pressure values at the iterative level v + 1, n + 1.
In the case of cells in direct contact with the well, we determine its pressure (p wf ) through Eq. (10) for a prescribed flow rate.
In solving the linearized system, to obtain the pressures in the porous medium and the well, it was chosen to use the Gauss-Seidel iterative method [16,36].

IV. NUMERICAL RESULTS
In all numerical simulations, we adopted an arrangement consisting of a vertical producer well, of length L wf and centered on the rθ-plane, a maximum production time (t max ), and an initial time step (∆t ini ). The time step can vary depending on the growth rate (δ ∆t ). This methodology allows a growing time step, and we use it until the final time step (∆t max ), preestablished, is reached. It is of general use in reservoir simulation [1,16]. We aim to enhance the accuracy of the calculated well pressure in the initial moments when the pressure drop is more accentuated.
We obtain the results using a standard set of data, based on the non-Darcy model discussed in Li et al. [29], in which the authors incorporated the slippage effect into the apparent permeability. Table 2 shows the parameters for the standard case, and we choose the physical properties as stated by Li et al. [29] (shale gas), and by de Souza [12], which simulated the flow of natural gas, using cylindrical coordinates, to assess

www.ijaers.com
Page | 230 the pressure of a producer well (vertical). As a simplification, we took C n+1 sc as constant and equal to C sc (we intend to modify this in future work).

Numerical verification
We employed different meshes in the study of numerical convergence: Meshes 1, 2, 3, and 4 with n r = 10, 20, 40, and 80 cells, respectively. We kept, in all simulations, n z constant and equal to 3. Figures 4 and 5 bring the results for the four meshes mentioned in the last paragraph. The results show the pressure in the well as a function of time, obtained considering a production time equal to 375 days. For this purpose, we use specialized and diagnostic graphics, respectively. By the way, in this work, on the diagnostic plots, continuous lines represent pressure drop, and dashed lines represent Bourdet derivative [8].
From the figures, we can see that we achieve numerical convergence as the number of cells n r increases, with the consequent overlap of pressure curves. Then, as a consequence of the results ob-tained, the 40 cells mesh was adopted as the standard (lowest computational cost). From the pressure curves in Fig. 4, we realize that we capture two typical regimes: a nearly horizontal line that corresponds to the storage in the well and an inclined straight line related to the transient flow regime. However, the effects of external borders are still absent, reflecting the fact that the pressure at the reservoir frontier (r = r e ) remains equal to the pressure p ini . From the Bourdet derivative (Fig. 5), we can more clearly distinguish the two distinct flow regimes. First, we observe the inclined line related to the wellbore storage and after a straight line corresponding to the transient flow regime, without boundary effects.   We also carried out tests to observe the behavior of the solution concerning the variation of the time step ∆t n+1 = δ ∆t ∆t n . We consider two situations for prescribed ∆t ini and ∆t max : numerical simulations

www.ijaers.com
Page | 231     Fig. 9: Well pressure variation due to the growth of ∆t, diagnostic plot.
As we can see from these figures, the results did not show significant variations concerning changes in the well pressure and pressure derivative when we employ these ∆t ini and δt values. Based on this verification, we established the standard values for the initial time step, the growth rate, and the maximum time step of Table 2. As already mentioned, the use of a small initial time step contributes to enabling the analysis of the pressure variation in the initial moments of production (typical in Well Test Analysis using numerical simulation).
We also studied the effect of the tolerance used to stipulate convergence in the internal (Gauss-Seidel method) and external (Picard method) iterations. Figures 10 and 11 present the results, which corroborate the choice of tolerance for the standard case since we did not detect significant differences in the results.  Fig. 11: Wellbore pressure variation as a function of tolerance, diagnostic plot.
As an additional check, Figs. 12 (specialized plot) and 13 (diagnostic plot) show our results obtained with higher values for porosity and permeability, 0.2 and 1.0 10 −3 Darcy respectively, and without taking into account the phenomena of slippage, formation damage, and wellbore storage. This simulation is the same as those performed by de Souza [12] for the flow considering only the classic Darcy law, whose results were verified through a direct confrontation with those obtained with the commercial simulator IMEX/CMG [11]. IMEX is a commercial simulator widely known and used in reservoir simulation. The values here numerically calculated accurately reproduce those determined by de Souza [12]. Therefore, indirectly, it was possible to validate our simulator in the case of flow governed by the classic Darcy law, without the effects of slippage, formation damage, and wellbore storage.

Sensitivity analysis
We begin by the sensitivity analysis of the effects on the wellbore pressure when we change the reservoir's permeability. In Figs. 14 and 15, we can observe the specialized and diagnostic plots for wellbore pressure variation. The higher the permeability, the lower the pressure drop, following the modified Darcy's law. It is worth mentioning that permeability has a significant impact on flow in porous media, and is often responsible for defining the economic viability of the reser-

www.ijaers.com
Page | 233 voir, including decision making on a hydraulic fracturing operation. In the case of reservoirs, where the slip flow regime occurs, the permeability value also influences Kn, and the Knudsen number increases if we decrease it, with a consequent increase in the apparent permeability value. We also noticed that the higher the permeability, the shorter the duration of the transition for the transient regime. In this test, and for all curves, we note that boundary effects are absent. We can verify this by inspecting the end of the transient flow regime curve, in the specialized plot, and the curve corresponding to the transient flow regime in the diagnostic plot (Bourdet derivative). Besides that, we can highlight that wellbore storage physically depends on permeability in a way that for lower permeability, wellbore storage effect is more prolonged.
Just as we did with permeability, we used different porosity values to check its influence in the wellbore pressure variation. We show the results in Figs. 16 and 17. However, for the porosity values that we retained, we only note small variations in the well pressure curves (p wf ). From the figures, we observe that higher porosity values lead to lower pressure drop for maximum production time, and this is due to the higher volume of gas that we can produce. Therefore, in this situation, the well pressure drops less for a fixed production flow rate. Porosity also influences the value of the Knudsen number but acting in the opposite direction. Regarding the boundary effects, the smaller the porosity, the shorter the transient flow regime. Moreover, as expected, more significant well-bore storage effects occur for lower porosity values, although we only capture slight differences for the values used in simulations.  On the other hand, changing the tortuosity leads only to slight changes in the pressure values for the range that we considered (see Figs. 18 and 19).
We note from the results in Fig. 18 a higher pressure drop as the value of tortuosity increases. This is consistent with the fact that the higher the tortuosity of the medium, the lower the apparent permeability must be and, therefore, the higher the flow resistance. From the results shown in these figures, we realize that we did not detect significant differences in the Bourdet derivative for this input simulation data.

Well Test Analysis
We now focus on the specific cases of well test analysis. The initial transient response, during a well test, comes from the pressure difference applied to the volume of fluid initially present in the well. Figure 20 shows the storage effect on pressure behavior in the well. The higher the storage effect, the smaller the representative range of the transient flow regime. This means that a well test must be long enough for the storage effects to disappear so that the flow recorded in the well test identifies the transient regime of reservoir production. In practice, it is from the well test that we obtain the data for the characterization of the well-reservoir system. The diagnostic plot (Fig. 21) shows that the wellbore storage strongly influences the results in the initial times, in both the pressure drop and pressure derivative. The higher the wellbore storage, the smaller the pressure drop and the later the transient flow regime will be.   We also notice significant results when comparing the pressure variation for the flow without (k a = k) and with the slippage effects (k a = f (Kn)k). The same simulation conditions and parameters were maintained for both models, varying only the property of interest.
Initially, only the responses resulting from the slip flow effect are studied, Fig. 22. Despite presenting a similar initial behavior (storage and transition to the transient regime), the model taking into account the phenomenon of slippage presents a lower pressure drop as production time increases. We can explain  Based on the results in Fig. 23, we note that the use of the classical Darcy's law (lower apparent permeability) leads to wellbore storage more significant in magnitude and duration.  We present now the results corresponding to three different flow rates in Fig. 24. We only change the flow rate values: -5,000, -10,000, and -20,000 scf/day. We call attention to the fact that we show the results in the form of the pressure variation (∆p wf ) divided by the production flow rate (Q sc ), a graph applied in the area of well test analysis. We perceive non-linear behavior through the variation of pressure curves as time progress. Figures 25 and 26 display pressure variation curves in a region damaged by contact with fluids external to the reservoir during drilling and completion, and, as a consequence, we have a decreased permeability in this region (it was applied the Hawkins' formation damage model [39]). We assumed that the damaged region measures 2.77 ft in all tests run.
After the beginning of the transient regime, we can verify a significant variation in pressure drop as the damage increases (the relationship between the permeability of the damaged region, k s , and the nondamaged region, k). Furthermore, we observe that the effect of storage has a longer duration and magnitude for higher formation damage values (as a result of reduced permeability around the wellbore).  Finally, to show the appearance of the border effects, even for a shale gas reservoir, we performed a numerical simulation using the data presented in Table 2, except for the following three values: t max =4,280 days, Q sc =-5,000 scf/day and L r =625 ft.
Therefore, at the end of 12 years of production and after the transient regime, we realize that the pressure curve changes its inclination (downward curvature), Fig. 27, as a result of the emergence of border effects. As we can see, the border effects can (in some cases) appear only after a long period of time. This fact shows the practical difficulty of carrying out well pressure tests, which can take from a few days, in general, to months when long-lasting. Finally, we can also detect the same effect in Fig. 28. It corresponds to the change in the curve slope in the final simulation times, after the horizontal line for the Bourdet derivative.  In all tests performed, the Knudsen number kept its value at Kn ≤ 10, outside the range of molecular flow.

V. CONCLUSION
We carried out a study to better understand singlephase flow in petroleum reservoirs, aiming to maximize hydrocarbon recovery. The numerical simulation allows testing different production scenarios in less time so that we can choose an optimized production plan that leads to economically viable exploration.
The results of this work showed the importance of applying a complete model (including the effects of slippage, formation damage and wellbore storage) to study single-phase flow in shale gas reservoirs. The use of the classic Darcy's law can lead to well pressure results that do not correspond to reality due to the non-consideration of these effects. This is essential information that must be taken into account if we are to have reservoirs producing and generating profits.
From pressure tests in vertical wells, we could obtain the reservoir properties by solving an inverse problem using the physical model of this work. Therefore, depending on the results, the use of horizontal wells or hydraulic fracturing could make production viable for some reservoirs.
As expected, the nonlinear behavior of the results was adequately detected, especially in cases where we varied the permeability values, which appear explicitly in the governing partial differential equation. Also, we captured the wellbore storage and formation damage effects, expanding the scope of this study, and we verified the influence of these phenomena on