A Semi-Analytical Method for Simulating Two-Phase Flow Performance of Horizontal Volatile Oil Wells in Fractured Carbonate Reservoirs

Two-phase flow behavior in fractured carbonate reservoirs was investigated due to the importance for geothermal and petroleum resource recovery, such as in this study the phase change of volatile oil. This paper presents a semi-analytical method for accurately modeling two-phase flow behavior and quickly predicting the production performance. The fractured carbonate reservoir was modeled with a dual-porosity model, and the phase change and two-phase flow were modeled using the black oil model. The production of the oil phase was obtained through linearizing and solving the mathematical model. The gas phase production was forecast using the producing gas-oil ratio (GOR), calculated using flowing material balance equations. By comparing the semi-analytical solution to the solution of the commercial numerical simulator and applying it to a field case, the accuracy and practicability of the proposed semi-analytical method could be validated. Based on the semi-analytical model, the influences of several critical parameters on production performance were also analyzed. The proposed model was shown to be efficient in evaluating two-phase production performance of horizontal volatile oil wells. Furthermore, the new technique is able to serve as a useful tool for analyzing two-phase production data and making forecasts for volatile oil wells in fractured carbonate reservoirs.


Introduction
With the degree of deepening exploration and adoption of advanced technology, the exploitation of fractured carbonate volatile oil reservoirs has been attracting considerable attention in recent years [1][2][3][4].Meanwhile, the application of horizontal wells not only increases the drainage area but also enhances production of a single well.Analyzing the performance of horizontal wells has become more and more popular among engineers [5][6][7][8][9][10][11][12].The two-phase flow of oil and gas is quite common and hard to simulate in volatile oil reservoirs.The fractured carbonate volatile oil reservoir is rather complicated because of its complex reservoir architecture and fluid properties.A big challenge is that the phase change, pressure dependent properties, and two-phase flow behavior and stress-dependent fracture permeability should be all considered when developing the mathematical model of fractured carbonated volatile oil reservoirs.The current method to establish horizontal volatile oil wells is mainly depletion.It is crucial to predict the production and analyze two-phase production data from this kind of well, and much research has been done on theoretical models of horizontal wells in fractured carbonate reservoirs [13][14][15][16][17].However, most of them are only applicable to the single-phase flow problem.Moreover, little to no research has been performed on the production performance of horizontal volatile oil wells in fractured carbonate reservoirs.
Recently, researchers became interested in the study of two-phase flow modeling of horizontal wells in volatile oil reservoirs [1][2][3][4]7].Empirical, semi-analytical and numerical methods have been widely used for two-phase production data analysis.Empirical methods [18][19][20][21][22][23][24] are commonly applied in the actual oil field to conduct history matching and make forecasts for production wells because they are easy to use.However, these methods are only suitable for specific conditions and generally inappropriate for the two-phase flow problem.Conventional methods have used numerical simulations [25][26][27][28][29] for modeling two-phase flow behavior of horizontal volatile oil wells, but gridding and simulating are time-consuming.
Modeling the two phase flow characteristics is a key part in the numerical simulation of petroleum development.There are mainly two semi-analytical methods to capture the two-phase flow behavior.The first method is based on Perrine's [30] and Martin's [31] assumption, in which two-phase flow could be reduced to single-phase flow by defining a two-phase pseudo-pressure.According to Perrine's approach [30], the relationship between saturation and pressure is constant and should be known in advance.In addition, laboratory PVT tests are needed to evaluation this relationship.Fetkovich [32], Raghavan et al. [33][34][35][36], as well as Ng and Aguilera [37] proposed a procedure to obtain the relationship from producing the gas-oil ratio.Actually, the flow equations are still mathematically nonlinear; this is because the relative permeability is a function of saturation.Also, oil density and the formation volume factor are functions of reservoir pressure.However, both the saturation and pressure in the reservoir are unknown parameters.The second method is by using the similarity method [38][39][40][41][42][43][44][45], which uses the Boltzmann transformation to turn the partial governing equations (PDEs) into ordinary differential equations (ODEs).O'Sullivan [38] firstly introduced the similarity method in the study of the two-phase flow problem.Bøe et al. [39] and Kissling et al. [40] extended O'Sullivan work by using the Boltzmann transformation to turn PDEs into ODEs.Ayala and Kouassi [41] and Zhang and Ayala [42] expanded the Bøe's approach to cover the transient linear flow condition.The advantages of using the similarity method lie in two aspects, the first is that the relationship between saturation and pressure does not need to be known in advance, and the second is that the solutions of the equations could be solved accurately if the PDEs could be transformed into ODEs [43][44][45].However, the similarity method is only suitable for transient linear flow regime, rather than dual porosity models.Recently, Clarkson and Qanbari [46,47] proposed an approximate semi-analytical method to analyze two-phase production data for fractured tight reservoirs, in which the producing gas-oil ratio (GOR) is used to forecast gas production.Generally, current research on the semi-analytical methods mostly focuses on the fractured tight reservoirs.Moreover, relatively less research on horizontal volatile wells in fractured carbonate reservoirs has been found in public literature.Therefore, proposing an effective semi-analytical method to capture the two-phase flow behavior of horizontal volatile oil wells in fractured carbonate reservoirs is of significant importance.
In this study, a novel semi-analytical method is presented to simulate the two-phase flow behavior of horizontal volatile oil wells in fractured carbonate reservoirs.A dual-porosity, black oil model considering phase change and two-phase flow is applied to model the fractured carbonate reservoirs.The model takes into account the complexities of phase change, pressure dependent PVT properties, two-phase flow behavior, and permeability stress sensitivity characteristics in the reservoir.The governing equations are linearized using a successive iteration method.With Laplace transform and separation of variable, one can obtain the analytical solution of the mathematical model.After getting the production of the oil phase from the previous time step, the average pressure and saturation within the investigated area are calculated with the flowing material balance equation.Calculations are then repeated for this time step to update the model parameters and the producing GOR with the Newton-Raphson scheme, followed by a novel proposed procedure for history matching of field production data and production forecasts.The semi-analytical method is benchmarked against the commercial simulator Eclipse.Furthermore, the effects of some critical parameters on production Energies 2018, 11, 2700 3 of 21 performance are also analyzed with the proposed model.Finally, a field example from the eastern Pre-Caspian basin is used to demonstrate the practicability of the method.

Model Assumption
As shown in Figure 1, the fractured carbonate reservoir is modeled with a classical dual-porosity model [15], which is composed of fracture and matrix systems.A radial cylindrical dual-porosity medium reservoir is considered in which a single horizontal well is located at the center, which completely penetrates the formation.The dual-porosity model assumes the fracture system is the flow pathway directly connected to the wellbore, fluids in the fracture system first flow into the horizontal wellbore, followed by the matrix-fracture inter-porosity flow.The volatile oil consists of two components, which are dead oil and dry gas.During depletion development, the process of gas dissolving and separating from the oil phase is considered, but oil dissolved in the gas phase is neglected.Some simplifying physical model assumptions for the derivation of the governing equation are listed as: The horizontal well is produced with constant bottom-hole pressure in the fractured carbonate reservoir, and the external boundaries of the top, bottom, and side are all assumed to be closed.The spherical matrix blocks are considered in the model to describe the matrix-fracture inter-porosity flow [48].Fluid flow follows the law of Darcy seepage, and stress-dependent fracture permeability is considered.Also, capillary and gravity forces are neglected to simplify the model.parameters on production performance are also analyzed with the proposed model.Finally, a field example from the eastern Pre-Caspian basin is used to demonstrate the practicability of the method.

Model Assumption
As shown in Figure 1, the fractured carbonate reservoir is modeled with a classical dual-porosity model [15], which is composed of fracture and matrix systems.A radial cylindrical dual-porosity medium reservoir is considered in which a single horizontal well is located at the center, which completely penetrates the formation.The dual-porosity model assumes the fracture system is the flow pathway directly connected to the wellbore, fluids in the fracture system first flow into the horizontal wellbore, followed by the matrix-fracture inter-porosity flow.The volatile oil consists of two components, which are dead oil and dry gas.During depletion development, the process of gas dissolving and separating from the oil phase is considered, but oil dissolved in the gas phase is neglected.Some simplifying physical model assumptions for the derivation of the governing equation are listed as: The horizontal well is produced with constant bottom-hole pressure in the fractured carbonate reservoir, and the external boundaries of the top, bottom, and side are all assumed to be closed.The spherical matrix blocks are considered in the model to describe the matrixfracture inter-porosity flow [48].Fluid flow follows the law of Darcy seepage, and stress-dependent fracture permeability is considered.Also, capillary and gravity forces are neglected to simplify the model.

Mathematical Model
In this paper, a semi-analytical method is proposed to simulate two-phase flow behavior and make forecasts of horizontal volatile oil wells in fractured carbonate reservoirs.The forecasted producing GOR is used to reduce the equations of the two-phase flow model, so only the flow equations of the oil phase are used to develop the mathematical model.The production rate of the oil phase is obtained with the solution of the mathematical model, and the production rate of the gas phase is calculated by using the producing GOR.In the following, a simplified black oil model is used to analyze two-phase performance from this kind of well.The process of gas dissolving and separating from the oil phase is considered, while oil dissolved in the gas phase is neglected.After some reasonable assumptions and simplification, the mathematical model was established in a radial cylindrical system.
The governing differential equation of the fracture system is expressed as follows [10]: where r is the radial distance in reservoir formation, m; z is the perpendicular distance from bottom, m; μo represents the oil viscosity, mPa·s; Bo represents the oil formation volume factor, fraction; pf represents the fracture-system pressure, MPa; ϕf represents the fracture-system porosity, fraction; kfh represents the formation permeability of fracture-system, mD; kfv is the perpendicular permeability

Mathematical Model
In this paper, a semi-analytical method is proposed to simulate two-phase flow behavior and make forecasts of horizontal volatile oil wells in fractured carbonate reservoirs.The forecasted producing GOR is used to reduce the equations of the two-phase flow model, so only the flow equations of the oil phase are used to develop the mathematical model.The production rate of the oil phase is obtained with the solution of the mathematical model, and the production rate of the gas phase is calculated by using the producing GOR.In the following, a simplified black oil model is used to analyze two-phase performance from this kind of well.The process of gas dissolving and separating from the oil phase is considered, while oil dissolved in the gas phase is neglected.After some reasonable assumptions and simplification, the mathematical model was established in a radial cylindrical system.
The governing differential equation of the fracture system is expressed as follows [10]: where r is the radial distance in reservoir formation, m; z is the perpendicular distance from bottom, m; µ o represents the oil viscosity, mPa•s; B o represents the oil formation volume factor, fraction; p f represents the fracture-system pressure, MPa; ϕ f represents the fracture-system porosity, fraction; k fh represents the formation permeability of fracture-system, mD; k fv is the perpendicular permeability of the fracture-system, mD; k fro represents the oil relative permeability of the fracture-system, fraction; q oMf represents the inter-porosity flow rate per unit matrix volume; α c represents the conversion factor, fraction; S o represents the oil saturation, fraction; t is the time, h.The inter-porosity flow rate, q oMf can be written as [48] where r 1 is the radius of a spherical matrix block, m; k Mro represents the oil relative permeability of matrix-system, fraction; k M represents the matrix-system permeability, mD; p M represents the matrix-system pressure, MPa; r M is the radial distance in a spherical matrix block, m.Taking Equation (2) into Equation ( 1), the flowing equation can be rewritten as where p i is the initial formation pressure, MPa.The external boundaries, including top, bottom, and side are assumed to be closed and given by ) ) Inner boundary condition where ε represents a variable in the z direction, m; z w represents the perpendicular distance of the horizontal well from the bottom, m; q represents the production rate from a point source dominated by a point source production, m 3 /d.The governing differential equation of the matrix system is expressed by where ϕ M represents the matrix-system porosity, fraction; k Mro represents the oil relative permeability of the matrix-system, fraction.External boundary condition

Solution to Mathematical Model
As mentioned above, the production of the oil phase is obtained by solving the mathematical model, and the output of the gas phase is forecast with the producing GOR.According to the established Energies 2018, 11, 2700 5 of 21 mathematical model, one can find that the flow equations of the oil phase are still nonlinear.This is because the relative permeability is a function of oil saturation.Meanwhile, oil density and oil formation volume factor are functions of reservoir pressure.However, both the saturation and pressure in the reservoir are unknown parameters.
To linearize the mathematical model, the period is divided into several steps, and the production is predicted step by step.At each time step, the saturation related parameters are treated explicitly and updated with the average pressure within the investigation area.Meanwhile, pressure dependent parameters are handled using pseudo-pressure.Thus, the model can be linearized and solved.Then the Laplace transformation is adopted to address the linearized model, and Stehfest [49] numerical inversion is used to calculate the production rate in real space.After obtaining the production rate of the oil phase from the previous time step, the average pressure and fluid saturations within the investigation area are obtained with the solution of coupled material balance equations.Calculations are then repeated for this time step to update the model parameters with the Newton-Raphson scheme.In the following, the linearization of flow equations and solution of the mathematical model are presented.

Linearization of Flow Equation
As the PVT properties such as fluid viscosity and volume factor of oil and gas phases are quite sensitive to formation pressure for volatile oil reservoirs, a transformed two-phase pseudo-pressure equation is adopted to capture pressure sensitive PVT properties.
For convenience, the pseudo-pressure drop is defined to make the equations homogeneous.
where m i represents the initial pseudo pressure, MPa 2 /(mPa•s); m f represents the pseudo pressure of the fracture-system, MPa 2 /(mPa•s); m M represents the pseudo pressure of the matrix-system, MPa 2 /(mPa•s); ∆m f is the difference between the initial pseudo pressure and the pseudo pressure of the fracture-system, MPa 2 /(mPa•s); ∆m M is the difference between the initial pseudo pressure and the pseudo pressure of the matrix-system, MPa 2 /(mPa•s).Note that saturation is handled explicitly, which means that average saturation is assumed to be constant in each time step.The values of saturation dependent parameters are updated in each step using the average saturation in the previous time step.
New pressure transmitting coefficients of the two systems are also introduced to simplify the form of the governing equation.1 where η represents the transmitting coefficient, m 2 /(ks); β c represents the conversion factor, fraction; p represents average pressure, MPa; S p represents the average saturation, fraction.
As with the analogous equations for the fracture and matrix system, it can be seen from Equation ( 14) that the pressure and saturation dependent parameters are also updated with the average pressure and saturation within the investigated area.In this way, η j could be regarded as a constant variable over a short period of several days and weeks.Consequently, the equations could be linearized and solved analytically.

Solution of Equations
Laplace transform is used to solve the established model and is defined as (15) where s represents the Laplace transform variable, fraction; ∆ m represents the variable in Laplace space.Therefore, for the fracture system, the governing differential equation is rewritten as External boundary condition Inner boundary condition For the matrix system, the governing differential equation is given by External boundary condition Inner boundary condition Combined with the separation variable method [9,10], the solution of the point source rate in the Laplace domain can be written as in which Energies 2018, 11, 2700 where q represents the point source rate in the Laplace domain, m 3 /d; h is the formation thickness, m; r e is the radial distance of the side external boundary, m; I 0 represents the modified Bessel function of the first kind, zero order; ∆m wf represents the difference between the initial pseudo-pressure and the pseudo-pressure of the bottom-hole, MPa 2 /(mPa•s); K 0 represents the modified Bessel function of the second kind, zero order; I 1 represents the modified Bessel function of the first kind, first order; K 1 represents the modified Bessel function of the second kind, first order.According to the superposition principle, the production rate of the oil phase can be obtained.
where q represents the production rate of the oil phase in the Laplace domain, m 3 /d; L is the horizontal wellbore length, m.One can apply the Stehfest [49] numerical inversion to obtain the production rate in real space.

Produced GOR Calculation
The solution derived in the previous section is only for oil phase of horizontal volatile oil wells in fractured carbonate reservoirs, the production of the gas phase is forecast by prediction of the producing GOR, and the method for prediction of the producing GOR is presented in this part.
As stated above, for volatile oil reservoirs, only gas dissolved in the oil phase is considered in this study with ignorance of the oil dissolved in the gas phase.The conventional method used to calculate the producing GOR is based on the assumption that there is a constant relationship between saturation and pressure [32][33][34].The equation for the producing GOR is expressed by where R s represents the solution gas-oil ratio; q gsc represents the surface gas production rate, 10 4 m 3 /d; q osc represents the surface oil production rate, m 3 /d.The producing GOR in Equation ( 28) is a function of pressure and saturation at the wellbore.The bottom-hole flowing pressure is easy to obtain, however, the saturation could not be calculated with analytical models.Therefore, in this approach, an approximate method was proposed to deal with this problem.The saturation dependent parameters in Equation ( 28) are updated with the average pressure calculated using flowing material balance equations within the investigated area.
The equation for the producing GOR can be rewritten as where GOR represents the producing gas-oil ratio, m 3 /m 3 ; k fro S o represents the relative permeability to oil at the average oil saturation; R s (p) represents the solution gas at average pressure.
To calculate the average saturation and pressure in Equation ( 29), the nonlinear equations must be solved iteratively for the oil and gas phases respectively.In the following, a coupled flowing material balance technique is developed to forecast oil and gas production.

Material-Balance Equations
In this section, the detailed derivation of material balance equations for a two-phase system is provided.The basic material balance equation for the oil phase is where IOIP is initial oil in place, m 3 ROIP is remaining oil in place, m 3 ; N p is cumulative oil production, m 3 .Initial oil in place Remaining oil in place Cumulative oil production Substituting Equations ( 31) through (33) into Equation (30) gives In this paper, the investigated radius for radial flow is calculated by use of [50] where r inv represents the investigated radius, m; k fi represents the fracture permeability of dual-porosity media, MPa.
Similarly, the basic material balance equation for gas is where IGIP is initial gas in place, 10 4 m 3 ; RGIP is remaining gas in place, 10 4 m 3 ; G p is cumulative gas production, 10 4 m 3 .Initial gas in place Remaining gas in place Cumulative gas production Energies 2018, 11, 2700 9 of 21 Substituting Equations ( 37) through (39) into Equation (36), one can obtain The saturation relation equation is expressed by Combining Equation (34) with Equation (40) and Equation ( 41) yields that In which where S oi and S gi are the saturations of oil and gas at initial pressure respectively; B oi and B gi are the formation volume factors of oil and gas at initial pressure respectively; R si is the solution gas oil ratio at initial pressure; S g and S o are the saturations of oil and gas at average pressure, respectively; B o and B g are the formation volume factors of oil and gas at average pressure, respectively.According to Equation (42), it is easy to deduce that Then average pressures using the Newton-Raphson scheme for Equation ( 45) is where k is the previous time step, and k+1 is the present time step; ω is the iterative coefficient.

History Matching and Production Prediction
This section provides a flow chart for two-phase history matching and production forecasting with the proposed semi-analytical method.The procedure for estimation of oil and gas production rate and history matching is shown in Figure 2. To conduct history matching and production prediction with the new method successfully, the reliability of model inputs such as reservoir parameters, relative permeability curves, fluid properties, and field data are very significant and should be checked carefully before being used in the model.The production is calculated with the semi-analytical model step by step, and the time increment is given by one day in accordance with the production data.Firstly, the initial average pressure and average saturation are used to calculate the production of the oil phase with Equation (27).The output of the gas phase is forecast with the producing GOR (Equation ( 29)).The average pressure and saturation within the investigated area are calculated with the flowing material balance equation.Calculations are then repeated for this time step using the Newton-Raphson scheme (Equation ( 46)).Thus, nonlinear parameters and GOR are updated with the average pressure and saturation at each time step.After that, the production rate curves of oil and gas phases using the proposed model are employed to match the field data.In this paper, we adopted the nonlinear least-square curve fitting method to conduct field data matching by using the semi-analytical method.

Model Verification
The proposed semi-analytical method is benchmarked against the commercial simulator Eclipse.A numerical model is built to simulate the two-phase flow of the horizontal volatile oil well.The input parameters for the semi-analytical and numerical model are identical, in which, the horizontal well has a wellbore length of 500 m and is located in the center of the formation.The reservoir diameter and thickness are 880 m and 38 m respectively.The initial reservoir pressure is 24.48 MPa and the bottom-hole pressure is fixed at 6 MPa during the production periods.The basic reservoir, horizontal well, and fluid properties are listed in Table 1.The initial pressure of the formation is higher than the bubble point of fluids, so the primary phase in the reservoir is the oil phase.In this case, the oil is assumed to be live oil, and the gas is assumed to be dead gas.The fluid properties and relative permeability curves used as input for the semi-analytical and numerical model are shown in Figures 3 and 4 respectively.Calculate q o , q g using Eq. ( 27) and Eq. ( 29) respectively .
Calculate ⎯p and⎯S o using the material balance method and Newton-Raphson scheme.
Input field data and plot q osc , q gsc versus t .
Match the field data using the model type curves.
Choose the appropriate parameters to calculate reservoir properties.
Average pressure, saturation Flowchart of history matching procedure using the proposed semi-analytical method.

Model Verification
The proposed semi-analytical method is benchmarked against the commercial simulator Eclipse.A numerical model is built to simulate the two-phase flow of the horizontal volatile oil well.The input parameters for the semi-analytical and numerical model are identical, in which, the horizontal well has a wellbore length of 500 m and is located in the center of the formation.The reservoir diameter and thickness are 880 m and 38 m respectively.The initial reservoir pressure is 24.48 MPa and the bottom-hole pressure is fixed at 6 MPa during the production periods.The basic reservoir, horizontal well, and fluid properties are listed in Table 1.The initial pressure of the formation is higher than the bubble point of fluids, so the primary phase in the reservoir is the oil phase.In this case, the oil is assumed to be live oil, and the gas is assumed to be dead gas.The fluid properties and relative permeability curves used as input for the semi-analytical and numerical model are shown in Figures 3 and 4 respectively.Figures 5 and 6 show the comparison of production rate and cumulative production between the semi-analytical method and the numerical simulation.As shown, a good match between the semianalytical model and the numerical simulation in both flow rate and cumulative production is obtained.The matching error between these two methods for oil production rate and cumulative oil production is less than 2%, which can be seen in Figures 5a and 6a.For the gas production rate and cumulative gas production in Figures 5b and 6b, the matching errors are 3% and 5% respectively.The    Figures 5 and 6 show the comparison of production rate and cumulative production between the semi-analytical method and the numerical simulation.As shown, a good match between the semianalytical model and the numerical simulation in both flow rate and cumulative production is obtained.The matching error between these two methods for oil production rate and cumulative oil production is less than 2%, which can be seen in Figures 5a and 6a.For the gas production rate and cumulative gas production in Figures 5b and 6b, the matching errors are 3% and 5% respectively.The matching error of the gas production rate and cumulative gas production is higher than the oil phase.Figures 5 and 6 show the comparison of production rate and cumulative production between the semi-analytical method and the numerical simulation.As shown, a good match between the semi-analytical model and the numerical simulation in both flow rate and cumulative production is obtained.The matching error between these two methods for oil production rate and cumulative oil production is less than 2%, which can be seen in Figures 5a and 6a.For the gas production rate and cumulative gas production in Figures 5b and 6b, the matching errors are 3% and 5% respectively.The matching error of the gas production rate and cumulative gas production is higher than the oil phase.This is because the PVT properties of the gas phase are sensitive to the pressure, but less sensitive for the oil phase.Meanwhile, to characterize the pressure sensitive PVT properties, the pseudo-pressure concept is applied in this paper.Although there is a slight mismatch in gas production rate and cumulative gas production, the error is appropriate for engineering purposes (less than 10%).Therefore, the proposed semi-analytical method can accurately simulate the two-phase flow behavior of horizontal volatile oil wells in fractured carbonate reservoirs.

Sensitivity Analysis
After validation, the influences of some critical parameters on gas and oil production performance are analyzed with the semi-analytical model.Table 2 shows the input parameters used for the sensitivity analysis.

Sensitivity Analysis
After validation, the influences of some critical parameters on gas and oil production performance are analyzed with the semi-analytical model.Table 2 shows the input parameters used for the sensitivity analysis.

Sensitivity Analysis
After validation, the influences of some critical parameters on gas and oil production performance are analyzed with the semi-analytical model.Table 2 shows the input parameters used for the sensitivity analysis.The permeability stress sensitivity of the fracture is expected to have a significant influence on the production performance in fractured carbonate reservoirs.Therefore, the permeability is evaluated with the average pressure calculated by using flowing material balance equations within the investigated area, which is expressed by where γ represents the permeability modulus, MPa −1 .
Figure 7 shows the influence of permeability stress sensitivity on production performance.One can see that permeability stress sensitivity affects all production periods, especially in the late period.As the permeability modulus increases, the production rate of both the oil and gas phases decreases greatly during the first year.However, those cases with stronger stress sensitivity have a higher production rate int the late period.This is because the fluid flow will be difficult and more oil and gas are left in the reservoir with the increase of the permeability modulus.The stress-sensitivity reflects the damage of the permeability and a stronger stress-sensitivity will increase the damage of the permeability.Consequently, the permeability stress sensitivity of the fracture decreases the cumulative production.
Energies 2018, 10, x FOR PEER REVIEW 13 of 21 The permeability stress sensitivity of the fracture is expected to have a significant influence on the production performance in fractured carbonate reservoirs.Therefore, the permeability is evaluated with the average pressure calculated by using flowing material balance equations within the investigated area, which is expressed by where γ represents the permeability modulus, MPa −1 .
Figure 7 shows the influence of permeability stress sensitivity on production performance.One can see that permeability stress sensitivity affects all production periods, especially in the late period.As the permeability modulus increases, the production rate of both the oil gas phases decreases greatly during the first year.However, those cases with stronger stress sensitivity have a higher production rate int the late period.This is because the fluid flow will be difficult and more oil and gas are left in the reservoir with the increase of the permeability modulus.The stress-sensitivity reflects the damage of the permeability and a stronger stress-sensitivity will increase the damage of the permeability.Consequently, the permeability stress sensitivity of the fracture decreases the cumulative production.

Fracture Porosity
Figure 8 shows the influence of the fracture porosity on the oil and gas production rate calculated by the semi-analytical model.It can be seen from Figure 8 that the production rate increases with the increase of the fracture porosity during the early and intermediate flow periods.However, the fracture porosity has little effect on the late production period.This is because the fracture porosity reflects the relative capacity of the fluid stored in the fracture system, a bigger fracture porosity is the response of the relative abundant reserves in the fracture system.The production drop should decrease to maintain the constant bottom-hole pressure of the production well when increasing the

Fracture Porosity
Figure 8 shows the influence of the fracture porosity on the oil and gas production rate calculated by the semi-analytical model.It can be seen from Figure 8 that the production rate increases with the Energies 2018, 11, 2700 14 of 21 increase of the fracture porosity during the early and intermediate flow periods.However, the fracture porosity has little effect on the late production period.This is because the fracture porosity reflects the relative capacity of the fluid stored in the fracture system, a bigger fracture porosity is the response of the relative abundant reserves in the fracture system.The production drop should decrease to maintain the constant bottom-hole pressure of the production well when increasing the fracture porosity.Meanwhile, compared with the fracture storing capacity of the affected region of the pressure wave, the increment of fracture storing capacity could be neglected during the late flow period.That is the reason why the fracture porosity has little effect on the late production period.

Radial Distance of External Boundary
Figure 9 shows the effect of the radial distance of the external boundary on the production rate.The value of the radial distance of the external boundary reflects the size of the drainage area.The results in Figure 9 indicate that the smaller the drainage area, the faster the production rate declined.As the radial distance of the external boundary increases, the production rates of both oil and gas phases increase greatly during the early and intermediate flow periods.This is because when the radial distance of the external boundary increases, the time of the external boundary response will be delayed, so the production drop should decrease to maintain the constant bottom-hole pressure of the well.In addition, the radial distance of the external boundary has little effect on the production rate before 40 days as shown in Figure 9.This is because the pressure wave is just propagated to the top and bottom boundary.As shown in Figure 10, the horizontal wellbore length affects the production rate during the whole process.Increasing the horizontal wellbore length will lead to an increase in production rate, particularly in the intermediate production period.It can be seen from Figure 11 that a greater

Radial Distance of External Boundary
Figure 9 shows the effect of the radial distance of the external boundary on the production rate.The value of the radial distance of the external boundary reflects the size of the drainage area.The results in Figure 9 indicate that the smaller the drainage area, the faster the production rate declined.As the radial distance of the external boundary increases, the production rates of both oil and gas phases increase greatly during the early and intermediate flow periods.This is because when the radial distance of the external boundary increases, the time of the external boundary response will be delayed, so the production drop should decrease to maintain the constant bottom-hole pressure of the well.In addition, the radial distance of the external boundary has little effect on the production rate before 40 days as shown in Figure 9.This is because the pressure wave is just propagated to the top and bottom boundary.

Radial Distance of External Boundary
Figure 9 shows the effect of the radial distance of the external boundary on the production rate.The value of the radial distance of the external boundary reflects the size of the drainage area.The results in Figure 9 indicate that the smaller the drainage area, the faster the production rate declined.As the radial distance of the external boundary increases, the production rates of both oil and gas phases increase greatly during the early and intermediate flow periods.This is because when the radial distance of the external boundary increases, the time of the external boundary response will be delayed, so the production drop should decrease to maintain the constant bottom-hole pressure of the well.In addition, the radial distance of the external boundary has little effect on the production rate before 40 days as shown in Figure 9.This is because the pressure wave is just propagated to the top and bottom boundary.As shown in Figure 10, the horizontal wellbore length affects the production rate during the whole process.Increasing the horizontal wellbore length will lead to an increase in production rate, particularly in the intermediate production period.It can be seen from Figure 11 that a greater reservoir thickness will generate a higher production rate during the early production periods.This is because the pressure drop will increase and the pressure wave needs more time to propagate to the closed boundary when increasing the value of the reservoir thickness.

Ratio of Horizontal Permeability to Perpendicular Permeability
Figure 12 shows the effect of the ratio of horizontal permeability to perpendicular permeability on the oil and gas production rate.It can be observed that a larger kh/kv will result in a higher production rate during the early production periods.This is because the fluid flow will be more difficult in the vertical direction and the pressure drop will increase when increasing kh/kv, so the larger the kh/kv, the later the pressure wave propagates to the top and bottom boundary, and the higher the production rate will be.

Ratio of Horizontal Permeability to Perpendicular Permeability
Figure 12 shows the effect of the ratio of horizontal permeability to perpendicular permeability on the oil and gas production rate.It can be observed that a larger kh/kv will result in a higher production rate during the early production periods.This is because the fluid flow will be more difficult in the vertical direction and the pressure drop will increase when increasing kh/kv, so the larger the kh/kv, the later the pressure wave propagates to the top and bottom boundary, and the higher the production rate will be.

Ratio of Horizontal Permeability to Perpendicular Permeability
Figure 12 shows the effect of the ratio of horizontal permeability to perpendicular permeability on the oil and gas production rate.It can be observed that a larger k h /k v will result in a higher production rate during the early production periods.This is because the fluid flow will be more difficult in the vertical direction and the pressure drop will increase when increasing k h /k v , so the larger the k h /k v , the later the pressure wave propagates to the top and bottom boundary, and the higher the production rate will be.

Figure 1 .
Figure 1.Description of the dual-porosity model in this paper: (a) horizontal well-scheme in formation; (b) dual-porosity model flow scheme; (c) dual-porosity model with spherical matrix.

Figure 1 .
Figure 1.Description of the dual-porosity model in this paper: (a) horizontal well-scheme in formation; (b) dual-porosity model flow scheme; (c) dual-porosity model with spherical matrix.

Figure 2 .
Figure 2. Flowchart of history matching procedure using the proposed semi-analytical method.

Figure 3 .
Figure 3. PVT properties used as input for the semi-analytical and numerical model: (a) oil volume factor and viscosity curves; (b) solution gas-oil ratio curve.

Figure 4 .
Figure 4. Relative permeability curve used as input for the semi-analytical and numerical model.

Figure 3 .
Figure 3. PVT properties used as input for the semi-analytical and numerical model: (a) oil volume factor and viscosity curves; (b) solution gas-oil ratio curve.

Figure 3 .
Figure 3. PVT properties used as input for the semi-analytical and numerical model: (a) oil volume factor and viscosity curves; (b) solution gas-oil ratio curve.

Figure 4 .
Figure 4. Relative permeability curve used as input for the semi-analytical and numerical model.

Figure 4 .
Figure 4. Relative permeability curve used as input for the semi-analytical and numerical model.

Figure 5 .Figure 6 .
Figure 5. Results comparison of production rate between numerical simulation and the semianalytical method for validation case: (a) log-log plot of oil production rate; (b) log-log plot of gas production rate.

Figure 5 .Figure 5 .Figure 6 .
Figure 5. Results comparison of production rate between numerical simulation and the semi-analytical method for validation case: (a) log-log plot of oil production rate; (b) log-log plot of gas production rate.

Figure 6 .
Figure 6.Results comparison of cumulative production between numerical simulation and the semi-analytical method for validation case: (a) Cartesian plot of cumulative oil production; (b) Cartesian plot of cumulative gas production.

Figure 7 .
Figure 7. Influence of permeability stress sensitivity of fracture on production rate: (a) influence of permeability stress sensitivity on oil production rate; (b) influence of permeability stress sensitivity on gas production rate.

Figure 7 .
Figure 7. Influence of permeability stress sensitivity of fracture on production rate: (a) influence of permeability stress sensitivity on oil production rate; (b) influence of permeability stress sensitivity on gas production rate.

Figure 8 .
Figure 8.Effect of fracture porosity on production rate: (a) effect of fracture porosity on oil production rate; (b) effect of fracture porosity on gas production rate.

Figure 9 .
Figure 9.Effect of radial distance of external boundary on production rate: (a) effect of radial distance of external boundary on oil production rate; (b) effect of radial distance of external boundary on gas production rate.3.2.4.Horizontal Wellbore Length and Reservoir ThicknessFigures10 and 11show the effects of horizontal wellbore length and reservoir thickness on the production dynamics of a horizontal well in fractured carbonate volatile oil reservoirs, respectively.As shown in Figure10, the horizontal wellbore length affects the production rate during the whole process.Increasing the horizontal wellbore length will lead to an increase in production rate, particularly in the intermediate production period.It can be seen from Figure11that a greater

Figure 8 .
Figure 8.Effect of fracture porosity on production rate: (a) effect of fracture porosity on oil production rate; (b) effect of fracture porosity on gas production rate.

Figure 8 .
Figure 8.Effect of fracture porosity on production rate: (a) effect of fracture porosity on oil production rate; (b) effect of fracture porosity on gas production rate.

Figure 9 .Figure 9 .
Figure 9.Effect of radial distance of external boundary on production rate: (a) effect of radial distance of external boundary on oil production rate; (b) effect of radial distance of external boundary on gas production rate.3.2.4.Horizontal Wellbore Length and Reservoir Thickness Figures 10 and 11 show the effects of horizontal wellbore length and reservoir thickness on the production dynamics of a horizontal well in fractured carbonate volatile oil reservoirs, respectively.As shown in Figure 10, the horizontal wellbore length affects the production rate during the whole

Energies 2018 ,
10, x FOR PEER REVIEW 15 of 21is because the pressure drop will increase and the pressure wave needs more time to propagate to the closed boundary when increasing the value of the reservoir thickness.

Figure 10 .Figure 11 .
Figure 10.Effect of length of horizontal well section on production rate: (a) effect of length of horizontal well section on oil production rate; (b) effect of length of horizontal well section on gas production rate.

Figure 10 .
Figure 10.Effect of length of horizontal well section on production rate: (a) effect of length of horizontal well section on oil production rate; (b) effect of length of horizontal well section on gas production rate.

Energies 2018 ,
10, x FOR PEER REVIEW 15 of 21is because the pressure drop will increase and the pressure wave needs more time to propagate to the closed boundary when increasing the value of the reservoir thickness.

Figure 10 .Figure 11 .
Figure 10.Effect of length of horizontal well section on production rate: (a) effect of length of horizontal well section on oil production rate; (b) effect of length of horizontal well section on gas production rate.

Figure 11 .
Figure 11.Effect of reservoir thickness on production rate: (a) effect of reservoir thickness on oil production rate; (b) effect of reservoir thickness on gas production rate.

Table 1 .
Parameters used for validation case.Set⎯p=p i ,⎯S o =S oi at initial time step and determine the time step t N .

Table 1 .
Parameters used for validation case.

Table 2 .
Basic parameters used for the sensitivity analysis.

Table 2 .
Basic parameters used for the sensitivity analysis.

Table 2 .
Basic parameters used for the sensitivity analysis.