Analysis of nonlinear effects in fluid flows through porous media

This article examines the nonlinear effects of fluid flow in a porous medium, governed by a new semi-analytical equation, from three aspects: equation derivation, experimental verification, and macroscale simulation modelling. The rigorous derivation of the new equation is presented with a semi-analytical approach in which the gas slippage effect and inertial forces are described. The latter effect is controlled by Fochheimer number, which is defined as a product of tortuosity and Reynolds number. The new equation successfully predicts the deviations from Darcy’s law in low-permeability media when the gas slippage effect occurs. The Klinkenberg gas slippage factor is obtained as a function of porous media’s structural parameter (porosity and intrinsic permeability) and gas property (mean free path of gas molecules). The equation validations are performed by core flow experiments for a wide range of reservoir properties, which yield good matching relationship between modelled and observed values. In addition, the proposed semi-analytical equation is used to simulate gas flow in the radial model.


Introduction
The continuing decline in conventional hydrocarbon resources accelerates the development of unconventional formations, such as tight sands, shales, coal beds, etc. Porous media permeability is an essential indicator in unconventional reservoirs and is frequently determined through laboratory experiments. It becomes standard to use gas instead of water as a fluid for permeability measurements because of the advantages of inert properties and low sensitivity to ambient temperature (Freeman and Bush 1983;Rodwell and Nash 1992;Tanikawa et al. 2009). Theoretically, the permeability value does not depend on the type of porous medium. However, numerous studies have shown that for gas flows in porous media with low permeability such as tight sandstones, coal-bed and shale gas, the permeability measured during the experiments (gas permeability) is greater than the intrinsic (liquid) permeability (k) and increases with decreasing average gas pressure p (Heid et al. 1950;Jones 1972;Jones and Owens 1980;Faulkner and Rutter 2000;Tanikawa et al. 2009).  conducted the first laboratory tests of oil-bearing rock samples to detect this phenomenon, the results of which were further reviewed by Muscat (1937). Overall, it was found that for high-permeable media, the differences between liquid and gas permeability were small. In contrast, these differences were significant for low and ultra-low (10 -19 -10 -22 m 2 ) permeable porous media. To distinguish it from the intrinsic permeability, the term of apparent gas permeability k a has been introduced, which is expressed by Klinkenberg (1941) as inversely proportional to the pore pressure: (1) k a = k 1 + b p Klinkenberg's research has shown that gas permeability is a function of the mean free path( ) , which in its term depends on pressure, temperature, and the nature of the gas. According to rigorous studies (Heid et al. 1950;Jones and Owens 1980;Sampath and Keighin 1982), the capillary channel diameter becomes comparable to the gas mean free path, the frequency of collisions of gas molecules with the walls of pore channels increases. This phenomenon can be described by assigning a non-zero velocity of the gas at the pore wall u wall (gas slippage effect) or by considering the pore channel with increased effective pore radius (r) (Zolotukhin and Ursin 2000) (Fig. 1).
Equation 1 is convenient and easy to use, but the required Klinkenberg coefficient (b) value approximated by Eq. 2 is difficult to determine for various porous systems. Consequently, researchers mainly focused on developing a better model for calculating the b value. The best-known and cited works proposing the Klinkenberg slip coefficient are summarized and shown in Table 1.
The Knudsen number characterizes the degree of influence of gas slippage in porous media (Kn) . Kn is a dimensionless parameter defined as the ratio of the mean free path of gas molecules ( ) to the average pore diameter (d): Here is commonly calculated by the following relationship: Gas flow in porous media can be classified into four regimes according to the classification by Schaaf and Chambre (1961): continuum fluid flow ( Kn ≤ 0.001 ); slip flow ( 0.001 < Kn ≤ 0.1 ); transition flow ( 0.1 < Kn ≤ 10 ) and free 2M molecular flow ( Kn > 10 ). Figure 2 schematically shows the transition of flow regimes based on the Knudsen number. Several analytical solutions (Wu and Pruess 1998;Innocentini and Pandolfelli 2001;Zhu et al. 2007;Hu et al. 2009;Civan 2010;Hayek 2015;Chen 2016) and numerical (Zhu et al. 2007;Li et al. 2016;Tao et al. 2016;Li and Sultan 2017) have been developed for characterizing the gas flow in tight porous media considering the Klinkenberg effect. While the developed numerical models can simulate gas flow under complex conditions accurately, analytical solutions provide a simple tool for describing gas flow properties in porous media. All the studies discussed above confirm that the Klinkenberg effect exists in a wide range of porous media and transport processes.
Like gas slippage, liquid slippage has also been confirmed to exist between liquid molecules and pore inner walls by experimental methods and molecular dynamics simulations (Afsharpoor and Javadpour 2016). Studies about liquid flow in nanotubes suggest that the slip boundary condition is related to wettability and significantly different from the no-slip boundary condition (Barrat and Bocquet 1999;Majumder et al. 2005;Wang et al. 2016;Afsharpoor and Javadpour 2016).
In addition to Klinkenberg effects, fluid flow through porous media may be affected by inertial forces (Tek et al. 1962;Dranchuk et al. 1969;Firoozabadi and Katz 1979;Skjetne 1995). At the same time, significant inertial flow usually occurs in formations with high permeability. Forchheimer (1901) discovered that Darcy's law underestimates the pressure drop in high-velocity gas flow in porous media. In this case, Darcy's law is generally corrected by a square term in seepage velocity (Forchheimer 1901) as: Different names have been used for the term in the literature. In this paper, we use the term "non-Darcy coefficient" when referring to as it appears in the Forchheimer equation. It is an empirical value representing the inertial resistance in a porous medium and depending on the pore   Wu et al. (1998) Zolotukhin and Ursin (2000) 9 b = (0.15 ± 0.06)k (−0.37±0.038) Tanikawa and Shimamoto (2006) Moghadam and Chalaturnyk (2014) Hooman et al. (2014)  Many previous studies have focused on non-Darcy coefficient determination, which was discussed in detail in our previous work (Zolotukhin and Gayubov 2021). Reynolds and Forchheimer numbers have been used by researchers to determine the onset of non-Darcy flow. Reynolds number is usually written as: As shown in Table 2, researchers used different parameters for characteristic length (d) based on preference and convenience. Besides, Table 2 contains information on critical Reynolds numbers, in the range of 0.10-1000, corresponding to the boundary between laminar and inertial flow regimes.
Forchhimeier number is another criterion for flow regime determination, initially defined by Green and Duwez (1951) as: This dimensionless number has the advantages of a clear definition and easy applicability. All involved parameters in Eq. 7 can be determined experimentally. Therefore, in our studies, we consider the Forchheimer number as a flow (6) Re = ud rule criterion instead of Reynolds number since the former includes both Re and tortuosity. This paper presents a set of semi-analytical solutions developed to analyze steady-state gas flow through porous media with Klinkenberg and Forchheimer effects. To demonstrate the application of the proposed technique, a new semi-analytical model is compared to existing experimental data in the literature to observe the excellent agreement. We also report on Klinkenberg and non-Darcy coefficients and compare our predictions with available experimental results in the literature.
The current paper is presented in three parts. In the first part, we briefly recall Darcy's, Forchheimer's and Klinkenberg's laws and consider detecting the transition between the ranges of their applicability. The initial theoretical equations are estimated, and empirical dependencies are proposed based on the study of published data. The second section of the paper discusses developing a new semi-analytical equation, which explains the deviations from Darcy's law. This new equation provides enhanced accuracy in interpreting the experimental data on permeability and provides a basis for improved reservoir simulation models. In the last part, comparisons between the predicted values, experimental values, and field measurement data were conducted to verify the proposed model's accuracy. The presented results advance our understanding of fluid flow in low and high-permeability media, where pre and post-Darcy are the dominant flow regimes.  Green and Duwez (1951) 7 Ma and Ruth (1993) 9 Re = ru Re c = 0.11 Thauvin and Mohanty (1998) 10 Zolotukhin and Gayubov (2021) Theoretical modelling For each of linear and radial flows in porous media, first, Darcy's law is presented. Then it is combined with pre-Darcy and post-Darcy fluxes to obtain the universal semianalytical equation. The appropriate boundary conditions are used to solve the differential equations and obtain pressure distribution and flux as a distance function.

Linear Darcy flow
The first mention of flow measurements through a porous medium was in 1856, when Henry Darcy published his water flow studies through a sand filter (Darcy 1856). The results of the experiment were generalized into empirical law that bears his name today. The law relates the flow rate to the hydraulic gradient through a linear proportionality coefficient, which Darcy called "perméabilité". Darcy claimed that this coefficient depends on the permeability of the sand layer. In 1933, Wyckoff, Botset, Muscat, and Reed published an article entitled "The Measurement of the Permeability of Porous Medium for Homogeneous Fluids", which outlined first time the definition of permeability and its physical concept in petroleum engineering (Wyckoff et al. 1933): No theoretical justification was given for the above equation and no references were made to earlier studies. In the 1950s, Hubbert provided a complete theoretical foundation to Darcy's empirical expression, deriving it from the general Navier-Stokes equation (Hubbert 1956). For the case of compressible steady-state flow at room temperature, Eq. 8 can be written as: Taking into account the constant mass flow rate and the perfect gas law Eq. 9 can be rewritten in the following form: Equation 10 allows determining the mass flow rate u D,m if the inlet and outlet pressures are specified (Dirichlet boundary conditions). Note that the values 0 and p 0 are referred to an arbitrarily selected point with index 0 (it could be referred to standard conditions or any other point).
The value of the mass flow rate u D,m enables us to determine the distribution of pressure, fluid density and the volumetric flow rate through the core sample. By integrating Eq. 9 from the inlet point x 1 to an arbitrary point x the pressure distribution can be obtained by the following expression:

Linear post-Darcy flow
For many years, Darcy's law was considered the fundamental equation governing the fluid flow in a porous medium. However, Darcy's law is valid for a certain range of fluxes and some basic assumptions (laminar single-phase flow, an isothermal condition, constant fluid viscosity, and no rockfluid interaction) for fluid flow in porous media (Hubbert 1956;Scheidegger 1960;Wu et al. 1998). Indeed, when the flow rate increases, the pressure drop is no longer proportional to the seepage velocity. To describe the nonlinear effect, a quadratic term was included in Eq. 8 by Dupuit (1863) andForchheimer (1901) to generalize the flow equation, i.e.: Based on the Forchheimer equation, a plotting method was developed to determine absolute permeability and inertial coefficient through linear regression of experimental data. In gas wells, the coefficient is usually determined by the analysis of multi-rate pressure test results. Unfortunately, such data measurements are expensive and not readily available in many cases. Thus, it is a regular practice to use theoretical and empirical correlations that can be derived by exploiting experimental values of the inertial coefficient .
Using the dimensional analysis and machine learning methods, discussed in detail in our previous studies (Zolotukhin and Gayubov 2019, 2021), it was possible to obtain the following modified Forchheimer equation for the compressible linear steady-state flow: where is the tortuosity, a parameter that characterizes a porous medium at a mesoscale and is determined based on laboratory studies conducted on cores.
Here, B 1 and are constants equal to 9.67 ⋅ 10 −5 and 0.47 , respectively. To obtain a semi-analytical expression for a flow rate, we provide the following expression for the non-Darcy coefficient, which is advantageous in that both hydraulic properties of the rock and the fluid properties are considered simultaneously: Figure 3 depicts a comparison of experimentally measured (Cooper et al. 1999, Muljadi et al. 2016) and predicted slippage factors using various correlations (Jones 1987;Liu et al. 1995;Coles and Hartman 1998;Friedel and Voigt 2006). The non-Darcy coefficients values ( ) are normalized by dividing by the experimentally measured values exp . The greater the distance between the normalized slippage factor and the dashed black line, the less accurate the calculation. For all samples, the non-Darcy coefficient was underestimated by Jones, Liu et al., and Coles and Hartman correlations. Furthermore, for permeabilities less than 215 mD, Friedel and Voigt's correlation underestimated the non-Darcy coefficient. The graph shows that the values obtained from Eq. 15 are much closer to unity than other correlations.
Considering constant mass flow rate and the perfect gas law, Eq. 13 can be rewritten as mass velocity: Rewriting Eq. 16 as quadratic equation and solving it for u F,m we obtain: As shown in our previous work (Zolotukhin and Gayubov 2019), the use of the dimensional analysis method can significantly reduce the number of variables involved in the analysis, improve the quality of the resulting solutions and at the same time simplify the obtained equations and increase the reliability of the analysis.
Since u F,m is now determined (Eqs. 18 and 19), the pressure distribution can also be found:

Linear post-Darcy and pre-Darcy flow
A common problem in technology and science is the derivation of simple equations describing complex phenomena. Complex control equations are often known, but they are too difficult to analyse. This section will describe a universal semi-analytical model that couples the modified-Forchheimer model (mentioned in Sect. 2.2) with the gas slippage effect. This model can be used both for fluid flow with low-pressure gradients when the effect of gas slippage effects occurs and for high-pressure gradients when inertial effects appear. It should be noted that analytical solutions are Fig. 3 Comparison of the experimental and the predicted non-Darcy coefficient by various correlations obtained under steady-state flow conditions, and the effects of adsorption and diffusion are not considered in them. Since gas is more compressible than solid and liquid, their permeability characteristics are different from one another. Therefore, the difference between the gas (apparent) and liquid permeability of most porous media can be accounted for by the Klinkenberg effect mentioned in the introduction (Sect. 1). Assuming that the apparent reservoir permeability is a function of pressure k a (p) , we can rewrite Eq. 13 in the following form: The tortuosity of the pore channels, calculated by Eq. 14, with consideration of apparent permeability k a (p) , can be expressed in the following form: Let us now define a formula for the parameter b . As earlier proposed (Zolotukhin and Ursin 2000), the equation for the Klinkenberg gas slippage factor can be written in the following form: The further estimates of the apparent permeability based on Eq. 23, made by the authors with the assistance of machine learning methods and dimension analysis give the following best-fit estimate: Equation 24 includes a parameter that considers the type of gas when the gas slip effect occurs, which has been ignored by many researchers. Rushing et al. (2004) report the steady-state and unsteady-state permeability measurement results for both helium and nitrogen. The results show that although gas slippage effect is generally more pronounced with helium compared to the nitrogen. The sensitivity of the permeability values to the flowing gas type is shown in different sections of Salahshoor and Fahes (2017), demonstrating the need to capture this effect in the calculations through some modification or corrections. Experimentally measured gas permeabilities can be misleading if the effect of gas type is not considered.
Comparison of the experimentally measured (Tessem 1980;Torsaeter et al. 1981;Wu et al. 1998) and the predicted slippage factors by various correlations (Tanikawa and Shimamoto 2006;Heid et al. 1950;Jones and Owens 1980) is shown in Fig. 4. The slippage factors values (b) are normalized by dividing by the experimentally measured values b exp . The farther normalized slippage factor is from the dashed black line (unity), the less accurate the calculation. Heid et al. and Tanikawa and Shimamoto correlations overestimated the gas slippage factor for permeabilities less than 1 µD. Also, Jones and Owens (1980)

Fig. 4
Comparison of the experimental and the predicted gas slippage factor by various correlations correlation overestimated the gas slippage factor for the entire experimental data range. According to the graph, the values obtained using Eq. 24 (Zolotukhin and Gayubov) are much closer to unity than other correlations. However, a more extensive data set will be more helpful to investigate this relationship further.
Substituting the obtained relation (Eq. 22) into the flow equation (Eq. 21), we obtain: Following the logic of the previous section, integration of the left-hand side of Eq. 25 yields: Rewriting Eq. 26 as quadratic equation and solving it for u KF,m we obtain: The sample's pressure distribution can be found from Eq. 26 using the following rearrangement of its terms:

Fluid flows at the macroscale
In addition to the mesoscale model discussed in Sect. 2, a macroscopic model (i.e., the reservoir scale model) for radial flow has also been developed. The problem concerns gas injection into a well in a large horizontal, uniform, and isothermal formation. A constant gas mass injection rate is imposed at the well, and the initial pressure is uniform throughout the formation. Herein, we estimate the influence of the main reservoir (reservoir porosity, permeability, and tortuosity), fluid (oil and gas viscosity, and density), and technological parameters (pressure drawdown and well spacing) on the efficiency of field development. Let us rewrite the equation of fluid flow (Eq. 21) in the radial geometry: Using the obvious relation for the radial flow u KFm = q KF,m S = q KF,m 2 hr and substituting it into Eq. 32, we obtain:

Separation of the variables and integration yields:
Integration of the left-hand side term of Eq. 34 yields: Integration of the right-hand side part of Eq. 34 results in: Approximate integration of the second term of Eq. 36 yields: we can rewrite the integral Eq. 37 as follows: and the result of the integration of Eq. 33 in the following form: Equation 40 is the basis for determining the fluid mass flow rate q KF,m . Rewriting it as a quadratic equation: and solving it for q KF,m we obtain: The pressure distribution within the drainage zone can be found from Eq. 40 using the following rearrangement of its terms. Denoting y = p(r) , we obtain: Solution of Eq. 44 can be written as follows:

Results and discussion
In Sect. 3, the newly derived semi-analytical model for gas slippage and inertia factors for low and high permeability porous media has been derived. In order to examine the validity of such an approach, the steady gas flow in porous media will be compared with the available data in the literature. Firstly, we will conduct this study at the mesoscale-a scale at which most of the experiments occur and where the principles of continuum mechanics are fully applicable. Then, the well production rates will be estimated at the macroscale level (10 2 -10 4 m).

Experimental data
This section aims to provide an overview of the laboratory study conducted by various researchers on low permeable (Wu et al. 1998;Amao 2007) and high-permeable cores (Tessem 1980;Torsaeter et al. 1981). Permeability was measured by the steady-state flow method when the differential pore pressure p 1 − p 2 is kept constant and outflow from samples is monitored. The upstream side of pore pressure p 1 was controlled by the gas regulator for gas permeability measurement, and fluid that flows out from the downstream of a specimen was released to atmospheric pressure p 0 . Nitrogen (M = 28.014 g/mol) at isothermal condition was injected in all experiments with different pressure gradients. Core number, intrinsic permeability, porosity and length of the cores, and their derivatives are shown in Table 3. Tortuosity of pore channels ( ) , non-Darcy coefficient ( ) and slip factor (b) are calculated using Eqs. 14, 15 and 24, respectively. In the next section, these data will be interpreted according to the traditional Darcy's law (Eq. 9) and the new semianalytical model (Eq. 21).

Boundary conditions
Usually, the following types of boundary conditions are used in solving differential equations: • Conditions of the first kind, where the inlet and outlet pressure values are specified at the flow region's boundaries (Dirichlet boundary conditions).
• Conditions of the second kind, where the pressure gradient or flow rate is set at the flow region's boundaries (Neumann boundary conditions). • Mixed conditions, when the first type's condition is set on one of the boundaries (inlet), and the second type's condition is set on the other (outlet).
In our simulation analysis, Dirichlet and mixed boundary conditions are the most appropriate for the considered tasks.

Dirichlet boundary conditions
The numerical model of one-dimensional gas flow problem is carried out for a rock column, the parameters of which are listed in Table 3. Figure 5 presents the comparison of flow velocity along the rock column from our numerical result and Darcy's law, which indicates that the semi-analytical solution proposed by this study is in good agreement with the experimental data. Due to the lack of space, we place only part of the numerical results in Fig. 5. The markers are the experimental data and corresponding dashed green lines are the fitting results using Eq. 21. At the lower permeability (high Klinkenberg factor b ), slippage effect on the gas velocity is significant. This effect disappears as the permeability at very high gas pressure is increased to 10 -12 m 2 , where the inertial factor has a significant effect. Published laboratory experiments have focused on the clastic rocks, though it is questionable whether other rocks (carbonate rocks, granite, fractured rock, etc.) show the same features.
The experimental data shown in Fig. 5 are divided into high-permeable (9.2; 2.3; 3.2; 5.1) and low-permeable (1.1; 1; 6; 9b; 9a; 36) samples, which will be analysed in detail below. In samples 9.2, 2.3 and 3.2, the inertial factor prevails over the Klinkenberg effect, and the model correctly describes the experimental data over the entire pressure gradient range. In the case of cores # 5.1 and 1.1, in which both nonlinear effects are equally manifested, our model gives an underestimated value of the flowrate ratios at the pressure gradient below 12 MPa/m. This might be caused by the assumption that all pore sizes in the specimen are the same, far from natural clastic rocks. Samples # 1, 6, 9b and 9a mainly exhibit the gas slip effect, which gives an overestimated value of the velocity compared to the Darcy equation. Our model accurately describes this nonlinear effect in low-permeability samples. For sample # 36 with the lowest permeability, the model does not completely reproduce the fluid flow rate with increasing pressure gradient. This discrepancy might be associated with such effects as adsorption and diffusion that are not included in the model. These phenomena can significantly impact ultra-low permeable porous  (Wu et al. 1998) with extremely low permeabilities (Amao 2007), the gas flow velocity is an order of magnitude higher than that calculated by Darcy's law. These results imply that the difference in the permeability measurement using nitrogen and liquid is more considerable for less permeable specimens (Klinkenberg effect). Figure 6 presents the comparison of pressure distribution ratios p KF ∕p D along the rock column calculated by using the new Eq. 30 p KF and Darcy's law p D . At the lower permeability (high Klinkenberg factor b ), slippage effect on the gas pressure reaches 6.5% (Fig. 6b). This effect disappears as the permeability at very high gas pressure is increased to 10 -12 m 2 (Fig. 6a). The ratio p KF ∕p D = 1 suggests that Darcy's law can be used to describe gas transport in porous media, and slippage gas effect can be ignored.
To investigate the influence of nonlinear effects on the gas flow, we plotted the ratio of flow rates with and without inertial and viscous forces u KF ∕u D along the rock column (x∕l) at the pressure drop along with the sample of 0.1 MPa. Figure 7b illustrates that nonlinear effect, namely the gas slip effect, is very salient for low and ultra-low permeable samples (high value of b ). For example, for the core with permeability 0.000017 mD, this difference reaches 12.3 times. For samples with high permeability, where inertial forces prevail, the difference reaches 12%.

Mixed boundary conditions
In this section, the gas flow velocity along the core samples is set as the Neumann boundary condition at the core outlet. Setting u KF,m = u D,m enables us to define the pressure distributions for our semi-analytical model and Darcy's law for different permeability. This condition makes it possible to correctly control the pressure drop along the rock column and prevent the flow rate overestimation, neglecting nonlinear effects. Figure 8 represents a plot of pressure distribution along the sample calculated using Eq. 30 and Darcy's law. The pressure distribution ratio for cores with high permeability (12.14-1132 mD) varies within 0.75-1.2 times, and for samples with low permeability (0.000017-5.05 mD) this ratio increases by the factor of 2. Interestingly, the sample's pressure distribution ratio curve with a permeability of 1132 mD is directed in the opposite direction compared to the less permeable samples.
Provided that the mass velocities are equal u KF,m = u D,m , the volumetric flow rates in the sample may differ significantly. This effect can be seen in Fig. 9, which shows the volumetric velocity ratios u KF ∕u D for various samples under mixed boundary conditions. The difference for samples with 12.14-1132 mD is in the range 0.8-1.3 times, and for samples with 0.000017-5.05 mD the difference reaches 0.45 times. Figure 9a shows a similar effect as in Fig. 8a, where the trend of the velocity ratio curve for the most highly permeable sample is quite different compared to the less permeable samples.

Dirichlet boundary conditions
Production from unconventional gas reservoirs usually involves complex mechanisms. According to different simulations, inaccurately estimation of the matrix permeability in unconventional gas reservoirs will lead to overestimating the initial production rate and more importantly, an overestimation of the long-term production (Bustin et al. 2008). Using the Darcy equation for such reservoirs can lead to   Table 4 an incorrect estimate in long-term production and financial planning. These overestimations can be significantly reduced by using the approaches proposed in this section.
In the above section, the gas flow model's validation and its implementation based on theoretical solutions are given when gas flow with Klinkenberg and Forchheimer effects is considered. In this section, the coupled model that has been numerically implemented in Sect. 3 is used to simulate gas migration in the radial model. This example concerns gas injection into a well in a large vertical, uniform, and isothermal formation, which follows Eq. 40. A constant gas mass injection rate is applied at the well, and the initial pressure is uniform throughout the formation. Herein, we estimate the influence of the main reservoir (reservoir porosity, permeability, and tortuosity), fluid (oil and gas viscosity, and density), and technological parameters (pressure drawdown and well spacing) on the efficiency of field development. The well spacing is defined as the area of the field per production well and expressed in acres. Some of the reservoir properties (permeability, porosity, tortuosity) are taken from Table 3. The rest of the reservoir and fluid characteristics, along with technological parameters, are summarized in Table 4. Figure 10 presents our numerical result, which compares the flow rate ratio distribution along with the sample from our semi-analytical model and Darcy equation. In the figure, the dimensionless flow rate is the ratio of the rate calculated using Eq. 42 q KF,m to the rate given by Darcy equation p D . Interestingly to note, that at the certain permeability value q KF,m ∕q D,m = 1 , which indicates that both nonlinear effects are mutually compensated, and Darcy's law can be used to describe gas transport in porous media. This phenomenon in a porous medium can be seen in the example of a synthetic rock sample with a permeability of 0.8 mD, in which q KF,m ∕q D,m = 1 . The ratio q KF,m ∕q D,m > 1 indicates that fluid flow is dominated by slippage effect (Fig. 10b). However, the inertial effect becomes significant when the permeability at very high gas pressure increases to 10 -12 m 2 (Fig. 10a). Figure 11 demonstrates the effect of the nonlinear flow on the pressure-distance distribution for radial pre-and post-Darcy flows in a homogeneous porous medium using Eq. 45 with Δp = 5MPa and parameters given in Table 4. Because the basic features of fluid flow are manifested in the nearwellbore zone, logarithmic grid coordinates are selected that allow the identification of those features and are calculated as follows: where n the number of intervals into which the reservoir is divided. The results shown in Fig. 11b reveal that the pressure distribution ratio p KF ∕p D for pre-Darcy flow is smaller than the one for Darcy flow under ultra-low permeabilities (0.000017-0.002 mD). On the other hand, the results in Fig. 11a indicate that more deviations from Darcy flow happen at high-permeabilities. Note that the most significant changes in calculated parameters are observed in the interval from 25 cm to 15 m. This interval is much larger than the usual wellbore zone (10-50 cm). Figure 12 shows the nonlinear effect on the radial preand post-Darcy flux q KF ∕q D versus distance ln r∕r w in a homogeneous porous medium using parameters given in Table 4. The results demonstrated in Fig. 12a reveal that the flux at the well for post-Darcy flow is smaller than the one for Darcy law. These results suggest that the use of Darcy's (62) ln r r w = ln r r w ⋅ i n i = 0, 1, … , n Fig. 11 Comparison of the pressure ratios distribution along the reservoir from our numerical result p KF and Darcy equation p D for highpermeable (a) and low-permeable (b) rocks law to predict the well production rate from a high-permeable formation (177.68-1132 mD) leads to significant overestimation of the rate for the entire life of a well.

Mixed boundary conditions
In this section, the gas flow rate is set as the Neumann wellbore boundary condition. Setting q KF,m = q D,m enables us to define the pressure distributions for semi-analytical model and Darcy's law for different permeability values. Figure 13 represents a plot of pressure distribution along the sample calculated using Eq. 45 and Darcy's law. The pressure distribution ratio for high permeability (12.14-1132 mD) varies up to 0.75 times, and for porous media with low permeability (0.000017-5.05 mD), this coefficient varies in the range of 1-2%. Figure 14 represents the volumetric flow rate ratios q KF ∕q D for various porous media under mixed boundary conditions. The difference for samples with 12.14-1132 mD could reach up to 35%, and for samples with 0.000017-5.05 mD is below 1.5%.