An Analytical Solution of the Pseudosteady State Productivity Index for the Fracture Geometry Optimization of Fractured Wells

The pseudosteady state productivity index is very important for evaluating the production from oil and gas wells. It is usually used as an objective function for the optimization of fractured wells. However, there is no analytical solution for it, especially when the proppant number of the fractured well is greater than 0.1. This paper extends the established fitting solution for proppant numbers less than 0.1 by introducing an explicit expression of the shape factor. It also proposes a new asymptotic solution based on the trilinear-flow model for proppant numbers greater than 0.1. The two solutions are combined to evaluate the pseudosteady state productivity index. The evaluation results are verified by the numerical method. The new solution can be directly used for fracture geometry optimization. The optimization results are consistent with those given by the unified fracture design (UFD) method. Using the analytical solution for the pseudosteady state productivity index, optimization results can be obtained for rectangular drainage areas with arbitrary aspect ratios without requiring any interpolation or extrapolation. Moreover, the new solution provides more rigorous optimization results than the UFD method, especially for fractured horizontal wells.


Introduction
Hydraulic fracturing is an important method currently used for low-permeability oil and gas reservoir stimulation and development.It can be used for vertical and horizontal wells.This paper focuses on well production within closed rectangular drainage areas.In this case, the well works under the transient and pseudosteady state.
Gringarten and Ramey [1] proposed a list of instantaneous Green's and source functions that can be used with Newman's product method to generate solutions for a wide variety of reservoir flow problems.These source functions have been used to study the unsteady-state pressure distributions created by wells with a single vertical fracture [2].By dividing the fracture length into multiple segments, the non-uniform flux along the fracture can be considered.However, the pressure drop along the fracture has been neglected in previous studies and the conductivity of the fracture has been assumed to be infinite.Pressure-transient behavior and inflow performance of horizontal wells with multiple fractures have been studied using source functions [3,4].However, the flux along the fracture was assumed to be uniform.To consider the pressure drop and non-uniform flux in the fracture, the fracture was considered as discrete infinite sink points and the flow in the fracture was assumed to be a steady Darcy linear flow [5][6][7].Lord Kelvin's point-source solution and superposition principle were Energies 2019, 12, 176 2 of 22 used to construct flux equations for all sink points.Cinco-Ley et al. [8] divided the fracture length into multiple segments and coupled a reservoir flow equation with a fracture flow equation to solve the pressure-transient behavior in finite-conductivity fractures.
Gringarten and Ramey provided instantaneous Green's and source functions.However, if the source is continuous and the flow rate is variable, then the source function must be integrated with respect to time, which is not easy to handle in the time domain.These problems can be solved by performing the Laplace transform of the integration [9][10][11].Chen et al. [12] extended Newman's product method into the Laplace domain.However, it was only valid for a limited range of conditions.Laplace transform and Fourier transform have been used directly to solve reservoir and fracture flow equations [13][14][15][16][17][18][19][20].However, these equations must be solved differently for different reservoir and fracture configurations.
Ozkan [21] proposed an extensive list of solutions for point, line, and slab sources in terms of the Laplace transform variable.These solutions were used to obtain pressure distributions and well responses for a wide variety of wellbore configurations and reservoir conditions.Without the use of integration with respect to time, it becomes easier to consider variable rate conditions such as constant pressure production [22,23].Laplace inversion can be done through the Stehfest algorithm [24].The line-source solution has been widely used for complicated applications of multiple fractured horizontal wells [25][26][27][28].Slab sources [29][30][31][32][33][34][35] and distributed volumetric sources (DVS) [36][37][38][39][40] have recently been proposed for pressure-transient analysis and productivity prediction.More generally, numerical methods have also been used for well performance evaluation, such as the finite element method (FEM) [41], the extended finite element method (XFEM) [42], the finite difference method (FDM) [43], and the boundary element method (BEM) [44].
However, in the above methods, the fracture must be divided into points or small segments.Thus, calculation is time consuming and it is difficult to obtain analytical or asymptotic solutions.Another generic method is to assume linear-flow in both the reservoir and fracture.Flow equations can be simplified and solved by Laplace transforms.first proposed bilinear-flow theory in an infinite reservoir.Lee and Brockenbrough [46] extended it to the trilinear-flow model.Ozkan et al. [47,48] obtained trilinear-flow solutions for closed reservoirs.Brown et al. [49] used the model in unconventional shale reservoirs and introduced many asymptotic solutions for pressure and its derivative.Wang et al. [50] incorporated fractal theory into the trilinear-flow model.Wang et al. [51] considered asymmetric configurations in multi-fractured horizontal wells.The trilinear-flow model has been further extended to the five-linear-flow model for more complex reservoirs [52][53][54][55].
For closed reservoirs, wells work under transient and pseudosteady states.During the pseudosteady state, the productivity index remains constant.So, it is usually used as the objective function to be maximized for the optimization of fractured wells [56][57][58].However, it is difficult to obtain an analytical formula for the productivity index under the pseudosteady state.Most approximate analytical formulas [59][60][61] neglect flow in the fracture.Guo et al. [62] divided the flow system into four regions and coupled the radial flow in the non-fractured region of the reservoir, linear flow toward the fractures in the fractured region, linear flow in the fracture, and radial flow in the fracture toward the horizontal wellbore.Zhang et al. [63] extended the above work to consider non-Darcy flow and fracture heterogeneity.However, the pressure condition was not rigorously continuous between adjacent regions and the flux distribution along the boundary was assumed to be uniform.
Economides et al. [64] proposed a fitting formula of the productivity index based on Cinco-Ley and Samaniego's work [45].However, the formula is only valid for square drainage areas and proppant numbers less than 0.1.Since the trilinear-flow model can provide analytical solutions in the Laplace space, asymptotic solutions can be obtained for short or long time periods, which is the basis of this paper.Solutions in the time space can then be readily derived by analytical Laplace inversions [49].Nevertheless, the flow in the reservoir-fracture system becomes pseudo-radial rather than trilinear with proppant numbers less than 0.1 [65].In other words, the trilinear-flow model is no longer valid with small proppant numbers.To this end, this paper extends the established fitting solution for proppant numbers less than 0.1 by introducing an explicit expression of the shape factor [66] and proposes a new asymptotic solution based on the trilinear-flow model [49] for proppant numbers greater than 0.1.The two solutions are combined to evaluate the pseudosteady state productivity index.
Ozkan [21] proposed an asymptotic solution for line sources in the Laplace space under the pseudosteady state.This solution has been used to obtain the productivity index of fractured horizontal wells [67].However, with this method, the fracture must be dispersed and large linear equations must be solved, which is time consuming.Nevertheless, the results of this method have been plotted as curves.Fitting formulas of the optimal dimensionless fracture conductivity and maximum dimensionless productivity index have been obtained for a given proppant number and reservoir aspect ratio [68][69][70].This optimization method is referred to as unified fracture design (UFD) [64].
Since no simple analytical formula of the productivity index exists for proppant numbers greater than 0.1, the fitting formulas of the optimal dimensionless fracture conductivity and maximum dimensionless productivity index given by the UFD method are rather complex.In particular, when the reservoir aspect ratio does not equal 1, several parameters in the fitting formulas must be extrapolated [65] and become incorrect if they are far away from the existing parameter values.Using the analytical solution of the pseudosteady state productivity index proposed in this work, optimization results can be obtained for rectangular drainage areas with arbitrary aspect ratios without any interpolation or extrapolation.
For fractured horizontal wells, the flow near the wellbore in the fracture becomes radial and an additional skin factor can be incorporated into the productivity index formula for corresponding vertical wells [71].Wei and Economides [72] used this skin factor for the fracture geometry optimization of horizontal wells in the UFD method.However, in the above work, this skin factor was assumed to be constant, yet it actually varies with the permeability of the reservoir and fracture, the fracture height and width, and the wellbore radius.Therefore, the existing fitting formulas for the optimal dimensionless fracture conductivity and maximum dimensionless productivity index given by the UFD method are no longer accurate for fractured horizontal wells.With the analytical solution of the pseudosteady state productivity index proposed in this work, more rigorous optimization results for fractured horizontal wells can be derived.

The Model
In this section, we propose an analytical solution of the dimensionless pseudosteady state productivity index for rectangular drainage areas with arbitrary aspect ratios, given a certain proppant number and dimensionless fracture conductivity.The proppant number can be determined once the proppant volume injected into the well is known [64].The proppant volume can be determined by the optimization method proposed in our previous paper [57].Thus, the proppant number is considered as a known constant in this work.Using the function of the dimensionless pseudosteady state productivity index, the optimal dimensionless fracture conductivity can be determined to maximize the dimensionless productivity index.
For proppant numbers less than 0.1, we extend the established fitting equation of the productivity index for rectangular drainage areas by introducing the equivalent proppant number and analytical expression of the shape factor, which are detailed in Appendix A.
For proppant numbers greater than 0.1, we derive a new asymptotic solution based on the trilinear-flow model.We first define the related dimensionless variables, and then start the derivation from the analytical solution of the dimensionless wellbore pressure in the Laplace space.Using some approximate relations, we obtain the asymptotic solution of dimensionless wellbore pressure for long time periods (i.e., the Laplace-transform parameter with respect to dimensionless time is small).Next, applying the analytical inversion of Laplace transform, we can obtain the expression for dimensionless wellbore pressure in the time space.According to the definition of the dimensionless productivity index, we can obtain its analytical form under the pseudosteady state.
In order to derive the relationships between the dimensionless productivity index, proppant number, and dimensionless fracture conductivity, we extend the definition of proppant number for rectangular drainage areas, which is given in Equation A40.Using the relationship between proppant number and dimensionless drainage size, we can obtain the final asymptotic solution of the dimensionless productivity index under the pseudosteady state for rectangular drainage areas with arbitrary aspect ratios, given a certain proppant number and dimensionless fracture conductivity value.The detail derivation is described in Appendix B.
The analytical solution of the productivity index can be summarized as follows.
(1) When N prop ≤ 0.1: where N prop is the proppant number and C f D is the dimensionless fracture conductivity.The fitting function f is given in Equation (A7) and the analytical expression of the shape factor C A is given in Equation (A9).(2) When N prop > 0.1: where k y is the aspect ratio of the rectangular drainage area defined in Appendix A.
In the preceding derivations, we have assumed 1D (linear) flow within the hydraulic fracture; that is, we have ignored the radial convergence of flow toward the wellbore within the hydraulic fracture for horizontal wells.Mukherjee and Economides [71] provided the following equation to compute the skin factor caused by radial flow choking within the fracture: where k f is the permeability of the propped fracture (md); k is the permeability of the reservoir (md); h is the thickness of the reservoir (m); w f is the average width of the propped fracture (m); and r w is the radius of the wellbore (m).Thus, the dimensionless productivity index for fractured horizontal wells can be described as follows:

Verification
Since Equation ( 1) is an extension of established theories, it does not need to be verified in this section.The optimization results obtained from both Equations ( 1) and (2) will be verified by comparing them with the UFD method in the next section.
Through numerical Laplace inversion [24], we obtain the exact solution in the time space for dimensionless wellbore pressure using Equation (A22).Combining this solution with Equations (A37) and (A38), we can obtain the dimensionless productivity index under both transient and pseudosteady states.The productivity index under the pseudosteady state is used for the verification of the new asymptotic solution in Equation ( 2).
We take the data in Table 1 as an example.The variable values in dimensionless form, SI units, and field units are shown in the table.Note that the values are different when using different units.Usually, we prefer the dimensionless form to avoid unit transform issues without affecting the analysis result.The dimensionless variables are defined in Appendices A and B. We take four aspect ratios, namely, 0.1, 0.4, 0.7, and 1.0, for the verification.The curves of the productivity index for the four aspect ratios under both transient and pseudosteady states are plotted in Figure 1.The constant pseudosteady state productivity indexes calculated by the new asymptotic solution are also plotted as dashed horizontal lines in Figure 1.The horizontal lines agree well with the straight segments of the corresponding curves.
Energies 2019, 12, x FOR PEER REVIEW 5 of 23 and (A38), we can obtain the dimensionless productivity index under both transient and pseudosteady states.The productivity index under the pseudosteady state is used for the verification of the new asymptotic solution in Equation ( 2).We take the data in Table 1 as an example.The variable values in dimensionless form, SI units, and field units are shown in the table.Note that the values are different when using different units.Usually, we prefer the dimensionless form to avoid unit transform issues without affecting the analysis result.The dimensionless variables are defined in Appendices A and B. We take four aspect ratios, namely, 0.1, 0.4, 0.7, and 1.0, for the verification.The curves of the productivity index for the four aspect ratios under both transient and pseudosteady states are plotted in Figure 1.The constant pseudosteady state productivity indexes calculated by the new asymptotic solution are also plotted as dashed horizontal lines in Figure 1.The horizontal lines agree well with the straight segments of the corresponding curves.Observing the horizontal lines representing the pseudosteady state productivity indexes for each aspect ratio, we find that the productivity indexes do not increase or decrease monotonously with aspect ratio.In Figure 1, the productivity index for aspect ratio 0.4 is the greatest, while that for aspect ratio 0.7 is the second-largest.However, the productivity index for aspect ratios 1.0 and 0.1 are smaller than those of the others.More interestingly, if we plot the curve of the productivity index with the aspect ratio in Figure 2, we find that the productivity index reaches a maximum value when the aspect ratio is about 0.3.The relationship between the productivity index and the aspect ratio is not the focus of this paper.Nevertheless, it could be studied in future work.Observing the horizontal lines representing the pseudosteady state productivity indexes for each aspect ratio, we find that the productivity indexes do not increase or decrease monotonously with aspect ratio.In Figure 1, the productivity index for aspect ratio 0.4 is the greatest, while that for aspect ratio 0.7 is the second-largest.However, the productivity index for aspect ratios 1.0 and 0.1 are smaller than those of the others.More interestingly, if we plot the curve of the productivity index with the aspect ratio in Figure 2, we find that the productivity index reaches a maximum value when the aspect ratio is about 0.3.The relationship between the productivity index and the aspect ratio is not the focus of this paper.Nevertheless, it could be studied in future work.

Comparison with the Unified Fracture Design Method
Given a certain proppant number, we can easily obtain the optimal dimensionless fracture conductivity to maximize the dimensionless productivity index for different aspect ratios of rectangular drainage areas using proposed Equations ( 1) and (2).
When the optimal dimensionless fracture conductivity is determined, the optimal fracture halflength and width can be obtained accordingly [57]: where  is the permeability of the propped fracture (md);  is the permeability of the reservoir (md); ℎ is the thickness of the reservoir which, in this model, equals the fracture height (m);  is the optimal dimensionless fracture conductivity; is the single wing volume of the propped fracture (m 3 ); and  is the total volume of the propped fracture (m 3 ).
We compare the optimization result for the optimal dimensionless fracture conductivity and maximum dimensionless productivity index obtained from the proposed analytical solutions with the result of the UFD method without considering the radial flow skin factor of fractured horizontal wells [65].The influence of the radial flow skin factor on the optimization result will be discussed in the next section.Both results are compared with the numerical result by the direct boundary element method [66,67].
We consider the following proppant numbers: 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, and 100.0.For aspect ratios, we divide the values into two groups.The first group comprises 0.1, 0.5, and 1.0, whereas the second comprises 0.05, 0.35, and 2.0.For the aspect ratios in the first group, all parameter values used in the UFD method (see Appendix C) can be found in related tables and no interpolation or extrapolation is needed.For aspect ratios in the second group, some of the parameter values cannot be found in the given tables and interpolation or extrapolation is needed.Interpolation or extrapolation will cause errors, which will be discussed in this section.However, fracture geometry optimization using the proposed analytical solutions does not need any interpolation or extrapolation and, hence, will not introduce extra errors.
Figure 3 shows the optimal dimensionless fracture conductivity (  ) and maximum dimensionless productivity index ( ) obtained by three methods in the drainage area with aspect

Comparison with the Unified Fracture Design Method
Given a certain proppant number, we can easily obtain the optimal dimensionless fracture conductivity to maximize the dimensionless productivity index for different aspect ratios of rectangular drainage areas using proposed Equations ( 1) and (2).
When the optimal dimensionless fracture conductivity is determined, the optimal fracture half-length and width can be obtained accordingly [57]: where k f is the permeability of the propped fracture (md); k is the permeability of the reservoir (md); h is the thickness of the reservoir which, in this model, equals the fracture height (m); C f Dopt is the optimal dimensionless fracture conductivity; 2 is the single wing volume of the propped fracture (m 3 ); and V p is the total volume of the propped fracture (m 3 ).
We compare the optimization result for the optimal dimensionless fracture conductivity and maximum dimensionless productivity index obtained from the proposed analytical solutions with the result of the UFD method without considering the radial flow skin factor of fractured horizontal wells [65].The influence of the radial flow skin factor on the optimization result will be discussed in the next section.Both results are compared with the numerical result by the direct boundary element method [66,67].
We consider the following proppant numbers: 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, and 100.0.For aspect ratios, we divide the values into two groups.The first group comprises 0.1, 0.5, and 1.0, whereas the second comprises 0.05, 0.35, and 2.0.For the aspect ratios in the first group, all parameter values used in the UFD method (see Appendix C) can be found in related tables and no interpolation or extrapolation is needed.For aspect ratios in the second group, some of the parameter values cannot be found in the given tables and interpolation or extrapolation is needed.Interpolation or extrapolation will cause errors, which will be discussed in this section.However, fracture geometry optimization using the proposed analytical solutions does not need any interpolation or extrapolation and, hence, will not introduce extra errors.
Figure 3 shows the optimal dimensionless fracture conductivity (C f Dopt ) and maximum dimensionless productivity index (J Dmax ) obtained by three methods in the drainage area with aspect ratio (k y ) equal to 1.0 for various proppant numbers (N prop ).Both the values of C f Dopt and J Dmax from the proposed method and the UFD method are consistent with those obtained by the numerical method.
proppant numbers.The maximum error of  with the UFD method is 6.67% and occurs when  = 1.0, whereas the maximum error of  with the proposed method is 7.14% and occurs when  = 10.0.In this case, the errors from both methods are not high.Table 3 lists the errors in  obtained from the UFD and proposed methods for several proppant numbers.The maximum error of  with the UFD method is 0.49% when  = 10.0, whereas that for the proposed method is 11.49% when  = 1.0.Note that proposed Equation ( 1) is valid for  less than 0.1 and Equation (2) applies for  larger than 0.1.With the increase of  , the reservoir-fracture flow system transitions from pseudo-radial flow to trilinear flow.Within the transition region, error occurs for the optimization using the proposed method.4 shows the curves of  and  with  = 0.5 for various values of  , while Figure 5 shows these for  = 0.1.Note that with the decrease of  , the  values obtained by both the proposed and UFD methods deviate from the true value given by the numerical method.Table 2 lists the errors in C f Dopt obtained from the UFD and proposed methods for several proppant numbers.The maximum error of C f Dopt with the UFD method is 6.67% and occurs when N prop = 1.0, whereas the maximum error of C f Dopt with the proposed method is 7.14% and occurs when N prop = 10.0.In this case, the errors from both methods are not high.Table 3 lists the errors in J Dmax obtained from the UFD and proposed methods for several proppant numbers.The maximum error of J Dmax with the UFD method is 0.49% when N prop = 10.0, whereas that for the proposed method is 11.49% when N prop = 1.0.Note that proposed Equation ( 1) is valid for N prop less than 0.1 and Equation (2) applies for N prop larger than 0.1.With the increase of N prop , the reservoir-fracture flow system transitions from pseudo-radial flow to trilinear flow.Within the transition region, error occurs for the optimization using the proposed method.Figure 4 shows the curves of C f Dopt and J Dmax with k y = 0.5 for various values of N prop , while Figure 5 shows these for k y = 0.1.Note that with the decrease of k y , the C f Dopt values obtained by both the proposed and UFD methods deviate from the true value given by the numerical method.For k y = 0.1, the error becomes obvious, especially within the transition region.However, the error in J Dmax from the proposed method disappears within the transition region when k y becomes small.For  = 0.1, the error becomes obvious, especially within the transition region.However, the error in  from the proposed method disappears within the transition region when  becomes small.All aspect ratios considered above are from the first group, and all parameter values used in the UFD method (see Appendix C) can be found in related tables.However, for the aspect ratios from the second group such as  = 0.05 or 2.0, the value of the shape factor ( ) must be extrapolated from Table A1, and the value of parameters , , , and  must be extrapolated from Table A2.For  = 0.35, those parameters must be interpolated from the related tables.Since we have no special explanation for Tables A1 and A2, we use linear interpolation and extrapolation in this work.
Figure 6 shows the curves of  and  with  = 0.05 for various values of  .Since  is extremely small, the values of  obtained from both the proposed and UFD methods deviate from the result given by the numerical method, which is similar to the case presented in Figure 5.As is shown in Table 4, the maximum error in  with the UFD method is For  = 0.1, the error becomes obvious, especially within the transition region.However, the error in  from the proposed method disappears within the transition region when  becomes small.All aspect ratios considered above are from the first group, and all parameter values used in the UFD method (see Appendix C) can be found in related tables.However, for the aspect ratios from the second group such as  = 0.05 or 2.0, the value of the shape factor ( ) must be extrapolated from Table A1, and the value of parameters , , , and  must be extrapolated from Table A2.For  = 0.35, those parameters must be interpolated from the related tables.Since we have no special explanation for Tables A1 and A2, we use linear interpolation and extrapolation in this work.
Figure 6 shows the curves of  and  with  = 0.05 for various values of  .Since  is extremely small, the values of  obtained from both the proposed and UFD methods deviate from the result given by the numerical method, which is similar to the case presented in Figure 5.As is shown in Table 4, the maximum error in  with the UFD method is All aspect ratios considered above are from the first group, and all parameter values used in the UFD method (see Appendix C) can be found in related tables.However, for the aspect ratios from the second group such as k y = 0.05 or 2.0, the value of the shape factor (C A ) must be extrapolated from Table A1, and the value of parameters a, b, c, and d must be extrapolated from Table A2.For k y = 0.35, those parameters must be interpolated from the related tables.Since we have no special explanation for Tables A1 and A2, we use linear interpolation and extrapolation in this work.
Figure 6 shows the curves of C f Dopt and J Dmax with k y = 0.05 for various values of N prop .Since k y is extremely small, the values of C f Dopt obtained from both the proposed and UFD methods deviate from the result given by the numerical method, which is similar to the case presented in Figure 5.As is shown in Table 4, the maximum error in C f Dopt with the UFD method is 153.96%,whereas it is Energies 2019, 12, 176 9 of 22 160.31% with the proposed method, both of which occur within the transition region.Nevertheless, the extrapolation used in the UFD method has no significant effect on C f Dopt .
However, because of the extrapolation of the shape factor (C A ), J Dmax becomes invalid when N prop ≤ 0.1, as listed in Table 5.Some values of J Dmax are less than 0, which is also invalid when N prop is 0.1-10.0.Even for valid J Dmax values when N prop > 10.0, the values deviate seriously from the true values given by the numerical method.The data deviation leads to a maximum J Dmax error of 401.18% when N prop = 100.0.For the proposed method, although the value of C f Dopt deviates from the true value seriously, the J Dmax error is low and reaches a maximum of 15.52% when N prop = 10.0.In this sense, J Dmax calculated by the proposed method is closer to the true value than that from the UFD method, which suffers from extrapolation error.
Energies 2019, 12, x FOR PEER REVIEW 9 of 23 153.96%, whereas it is 160.31% with the proposed method, both of which occur within the transition region.Nevertheless, the extrapolation used in the UFD method has no significant effect on  .However, because of the extrapolation of the shape factor ( ),  becomes invalid when  ≤ 0.1, as listed in Table 5.Some values of  are less than 0, which is also invalid when  is 0.1-10.0.Even for valid  values when  > 10.0, the values deviate seriously from the true values given by the numerical method.The data deviation leads to a maximum  error of 401.18% when  = 100.0.For the proposed method, although the value of  deviates from the true value seriously, the  error is low and reaches a maximum of 15.52% when  = 10.0.In this sense,  calculated by the proposed method is closer to the true value than that from the UFD method, which suffers from extrapolation error.8 shows these with  = 2.0.Note that the interpolation for  = 0.35 and the extrapolation for  = 2.0 used in the UFD method have no significant effects on  , whereas they do have serious effects on  , especially when  is large.Nevertheless, interpolation induces less error than extrapolation in the UFD method.However, the proposed method does not experience any interpolation or extrapolation problem.The only error associated with the proposed  Figure 7 shows the curves of C f Dopt and J Dmax with k y = 0.35 for various values of N prop , while Figure 8 shows these with k y = 2.0.Note that the interpolation for k y = 0.35 and the extrapolation for k y = 2.0 used in the UFD method have no significant effects on C f Dopt , whereas they do have serious effects on J Dmax , especially when N prop is large.Nevertheless, interpolation induces less error than extrapolation in the UFD method.However, the proposed method does not experience any interpolation or extrapolation problem.The only error associated with the proposed method for the aspect ratio of the drainage area becomes too small or too large, flow in the reservoir-fracture system deviates from the trilinear-flow model and new models should be applied together with the current model in the future study.

Considering the Radial Flow Skin Factor
Note that from Equations ( 3) and ( 4), the dimensionless productivity index of fractured horizontal wells depends not only on the dimensionless fracture conductivity, proppant number, and aspect ratio, but also on the permeability of the fracture and the reservoir, the thickness of the reservoir, the radius of the wellbore, and the average width of the fracture.Even given a certain proppant number and aspect ratio, the dimensionless productivity index cannot be expressed as a single-valued function of dimensionless fracture conductivity.Thus, the formulas for optimal dimensionless fracture conductivity given by the UFD method in Equations (A45) and (A46) are no longer accurate, and no optimization solution exists, especially for fractured horizontal wells.
However, thanks to the analytical solutions of the dimensionless productivity index proposed in Equations ( 1)-( 4), it is possible to obtain the optimization result for fractured horizontal wells.The proppant number for rectangular reservoirs is defined as follows [65]: where x e is the length of the drainage area parallel to the direction of the fracture (m) and y e is the length of the drainage area vertical to the direction of the fracture (m), as illustrated in Figure A1.
Combining Equations ( 6) and ( 7), we can obtain: Substituting Equation (8) into Equation (3), the skin factor becomes: Substituting Equation ( 9) into Equation (4), we can obtain the analytical formula of the dimensionless productivity index for fractured horizontal wells expressed in terms of dimensionless fracture conductivity and proppant number.Thus, given a certain proppant number and the related reservoir-well parameters, the optimal dimensionless fracture conductivity and maximum dimensionless productivity index can be obtained.
For example, setting the following parameters: h = 20 m, r w = 0.1 m, and x e = y e = 1200 m, we can plot the relationship between optimal dimensionless fracture conductivity and proppant number, as is shown in Figure 9.The optimal dimensionless fracture conductivity is compared with the results from the UFD and numerical methods.Note that since the skin factor is taken as a constant in the UFD method, the optimal results from the UFD method deviate seriously from the true values of the numerical method, especially when N prop is less than 0.1, whereas the results from the proposed method agree well with the true values.

Conclusions
This paper extended the established fitting solution for the pseudosteady state productivity index of fractured wells when the proppant number is less than 0.1 by introducing an explicit expression of the shape factor.It also proposed a new asymptotic solution based on the trilinear-flow model for proppant numbers greater than 0.1.The two solutions were combined together in order to evaluate the pseudosteady state productivity index.The new solution can be directly used for fracture geometry optimization.Based on the analyses contained in this paper, the following conclusions can be drawn: (1) The proposed analytical formulas provide solutions of the dimensionless productivity index for any dimensionless fracture conductivity value and proppant number, whereas the UFD method only provides fitting solutions of the maximum productivity index corresponding to the optimal fracture conductivity.(2) Using the proposed analytical solution of the pseudosteady state productivity index, optimized fracture geometry dimensions can be obtained for arbitrary aspect ratios of rectangular drainage areas, whereas the UFD method suffers from interpolation or extrapolation error if the value of the aspect ratio is not given in related tables.(3) By explicitly considering the skin factor caused by the radial flow choking of fractured horizontal wells in dimensionless fracture conductivity optimization, the result is more accurate than that obtained by the UFD method.
Although most of the proposed method's optimization results agree well with the true values given by the numerical method, there are errors associated with some ranges of rectangular drainage area aspect ratios.Determining more accurate and adaptable analytical solutions that overcome such shortcomings is the focus of future work.
Figure A1a represents a fractured vertical well in a square drainage area.Figure A1b is a fractured vertical well in a rectangular drainage area with  < 1. Figure A1c is a fractured vertical well in a rectangular drainage area with  > 1. Figure A1d is a horizontal well intersected with a single fracture in a square drainage area.Figure A1e is a horizontal well intersected with multiple fractures in a rectangular drainage area.Each fracture is in an individual drainage area with  ≤ 1 or  > 1.We first consider the square drainage area and define the penetration ratio as [64,68]: and dimensionless fracture conductivity as: where  is the permeability of the propped fracture (md);  is the permeability of the reservoir (md); and  is the average width of the propped fracture (m).We define the proppant number as [64]: where  is the volume of the propped fracture (m 3 ) and  is the volume of the drainage area (m 3 ).If the volume of the propped fracture is constant, then the proppant number is also constant.The We first consider the square drainage area and define the penetration ratio as [64,68]: and dimensionless fracture conductivity as: where k f is the permeability of the propped fracture (md); k is the permeability of the reservoir (md); and w f is the average width of the propped fracture (m).We define the proppant number as [64]: where V p is the volume of the propped fracture (m 3 ) and V res is the volume of the drainage area (m 3 ).If the volume of the propped fracture is constant, then the proppant number is also constant.The volume of the propped fracture is related to the injected proppant volume, which can be determined by the optimization method proposed in our previous paper [57].This paper studies the relationship between the dimensionless fracture conductivity and dimensionless productivity index for a certain proppant number.From Equations (A2) and (A3), the fracture half-length can be described as: When the proppant number is less than 0.1, the dimensionless productivity index for the square drainage area can be described as [64]: Substituting Equation (A4) into (A5), the dimensionless productivity index becomes: For the rectangular drainage area, the proppant number can be replaced by the equivalent proppant number in Equation (A8).
We further extend the analytical form of the dimensionless productivity index by introducing an explicit expression of the shape factor C A [66]  where the Euler constant is γ = 0.57721566 and ε is a suitably small positive number (e.g., 10 −6 ).
Taking the lower left corner of the drainage area as the origin of the coordinate system, x w and y w are the coordinates of the wellbore, which is located in the center of the drainage area while calculating the shape factor.a is the influence function given by [66,68]: function of the maximum dimensionless productivity index or the dimensionless fracture conductivity with the proppant number can be fit out.We implemented the numerical method and illustrated the curves in Figure A2 for aspect ratio = 0.5 as an example.Without loss of generality, we showed the curves for proppant number = 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, and 100.0 [65].The maximum dimensionless productivity index and the optimal dimensionless fracture conductivity in rectangular reservoirs can be described by the following equations [57,65].The maximum dimensionless productivity index and the optimal dimensionless fracture conductivity in rectangular reservoirs can be described by the following equations [57,65].The Dietz shape factors are taken from Table A1.The related constants are taken from Table A2.

Figure 1 .
Figure 1.Comparison of true values and asymptotic solutions.

Figure 1 .
Figure 1.Comparison of true values and asymptotic solutions.

Figure 2 .
Figure 2. Curve of the pseudosteady state dimensionless productivity index vs.aspect ratio.

Figure A1 .
Figure A1.Geometric relationships of reservoirs and fractures.(a) A fractured vertical well with  = 1; (b) A fractured vertical well with  < 1; (c) A fractured vertical well with  > 1; (d) A horizontal well intersected with a single fracture with  = 1; (e) A horizontal well intersected with multiple fractures.

Figure A1 .
Figure A1.Geometric relationships of reservoirs and fractures.(a) A fractured vertical well with k y = 1; (b) A fractured vertical well with k y < 1; (c) A fractured vertical well with k y > 1; (d) A horizontal well intersected with a single fracture with k y = 1; (e) A horizontal well intersected with multiple fractures.

Figure A2 .
Figure A2.Relationships between dimensionless productivity index and dimensionless fracture conductivity.
Figure A2.Relationships between dimensionless productivity index and dimensionless fracture conductivity.

Table 1 .
Variable values used for verification.

Table 1 .
Variable values used for verification.
Curve of the pseudosteady state dimensionless productivity index vs.aspect ratio.

Table A2 .
Constants in F-function.