Optimization Algorithm of Effective Stress Coefﬁcient for Permeability

: The effective stress coefﬁcient for permeability is a signiﬁcant index for characterizing the variation in permeability with effective stress. The realization of its accuracy is essential for studying the stress sensitivity of oil and gas reservoirs. The determination of the effective stress coefﬁcient for permeability can be mainly evaluated using the cross-plotting or response surface method. Both methods preprocess experimental data and preset a speciﬁc function relation, resulting in deviation in the calculation results. To improve the calculation accuracy of the effective stress coefﬁcient for permeability, a 3D surface ﬁtting calculation method was proposed according to the linear effective stress law and continuity hypothesis. The statistical parameters of the aforementioned three methods were compared, and the results showed that the three-dimensional (3D) surface ﬁtting method had the advantages of a high correlation coefﬁcient, low root mean square error, and low residual error. The principal of using the 3D surface ﬁtting method to calculate the effective stress coefﬁcient of permeability was to evaluate the inﬂuence of two independent variables on a dependent variable by means of a 3D nonlinear regression. Therefore, the method could be applied to studying the relationship between other physical properties and effective stress.


Introduction
Effective stress has always been a notable research topic in pore elastomechanics theory, which has substantial guiding significance for oil and gas development and is widely used in reservoir numerical simulations, oil and gas productivity prediction, and reservoir reconstruction.
Terzaghi [1] first proposed the expression of effective stress as follows: σ T e f f = σ − P, Equation (1) is suitable for studying the settlement of loose soil and mechanical failure of most rocks. However, the cementation mode, pore structure, and compressibility of rock in oil and gas reservoirs differ from those of soil, and thus Equation (1) is not suitable to express effective stress.
According to the theory of Biot [2,3] and the principle of stress-strain superposition, Nur [4] derived an accurate expression of the effective stress applicable to the elastic deformation of pores: where α is the Biot coefficient, σ eff is the effective stress (MPa), σ is the overlying stress or confining pressure (MPa), and P is the pore fluid pressure (MPa). Many theoretical and experimental studies have shown that the predicted value of the Biot coefficient is less than 1 [3,5,6]. For rocks in general, the Biot coefficient is between 0.12 and 0.91. The Biot coefficient of effective stress mainly focuses on measuring the pore elastic parameters of rock and represents the degree of the relative influence of pore fluid pressure and overburdened pressure on the pore elastic parameters.
However, the effective stress coefficient is less than 1, which is not applicable to all types of rock. Ghabezloo [7], Meng [8], Wang [9], Li Min and Xiao Wenlian et al. [10] found that the effective stress coefficient was greater than 1 when they tested the permeability of rocks containing clay minerals, which indicated that the influence of pore fluid pressure on permeability was greater than confining pressure.
According to various scholars including Berryman [11], Robin [12], and Bernabe [13,14], there is no uniform law of effective stress for different physical properties, and the effective stress coefficients differ significantly. Therefore, different symbols should be used to represent the effective stress coefficients for different physical properties. When using permeability as the research object, such coefficient should be called the effective stress coefficient for permeability (ESCP, expressed as α k to distinguish it from the Biot coefficient).
Bernabe [13,14] first proposed the concept of effective stress for permeability and then expressed the relationship among permeability and effective stress, confining pressure, and pore pressure as follows: where α k is ESCP, which reflects the relative magnitude of the influence of pore fluid pressure and confining pressure on permeability. The calculation of the ESCP is helpful for studying the deformation law and seepage characteristics of reservoir rock. This calculation is the basis of the evaluation of reservoir stress sensitivity and is of great significance for oil and gas productivity evaluation and efficient development [15,16].

Calculation Methods for the ESCP
The determination of the ESCP depends on the test data obtained from the stresssensitive experiment. The ESCP can be obtained after the data are processed using the cross-plotting method [17] or the response surface method [18].

Cross-Plotting Method
Walsh [17] used the cross-plotting method to calculate the effective stress coefficient while studying the permeability of fractured rock that changed with the confining and pore fluid pressures. The main steps of this method are as follows.
(1) Linear fitting was performed for determining the relationship between k − 1 3 and pore fluid pressure under different confining pressures.
(2) In the linear diagram of the relationship between k − 1 3 and pore fluid pressure, isolines of several groups of permeability were drawn, and the confining pressure and pore fluid pressure under the same permeability could be obtained by cross-plotting.
(3) The relationship between the confining pressure and pore fluid pressure under the same permeability was linearly fitted, and the slope of the linear relationship between the confining pressure and pore fluid pressure was determined as the effective stress coefficient. The effective stress coefficient of the rock sample was obtained by averaging several groups of the ESCP.
The cross-plotting method implies two prerequisites: first, the effective stress should be linear, that is, Equation (2), and second, the relationship of permeability with effective stress was preset as follows [17]: where k 0 and a 0 are the permeability and half-width of the crack at a reference stress σ e f f 0 , respectively, and h is the root mean square of the height distribution of the crack surface. Jones [19] demonstrated Equation (4) in an experimental study of fractured carbonate rocks. Equation (4) was derived by Walsh [17,20] and is based on the plate fracture model, considering fracture roughness. Therefore, Equation (4) is mainly applicable to fractured rocks.

Response Surface Method
Xiao [18] and Li [21,22] used the response surface method to investigate the stress sensitivity and effective stress coefficient of tight sandstone. The steps of this method are as follows.
(1) Transform the permeability or volumetric strain data into a simpler form and appropriately weight the variance.
where k is the original experimental data of permeability, k is the converted permeability data, and λ demonstrates power, which is typically between −3 and +3. A transformation can simplify the surface or standardize the known variances or errors in the data; however, this must be performed carefully to avoid skewing the statistical analysis.
(2) Fit the converted data to the quadric surface for both σ and P: where x i coefficients are determined from the fit.
(3) Statistical parameters and quadric surface figures were used to infer the physical and mechanical effects, permeability evolution, and effective stress coefficients of rock.
Equation (6) could be simplified as follows: where y i coefficients are determined from the fit. Equation (6) can be derived by expanding Equation (7). By introducing an intermediate variable (i.e., effective stress σ k ), Equations (6) and (7) can be written as follows: Equations (8) and (9) demonstrate that the response surface method implies two conditions for use: first, the effective stress should be linear, that is, Equation (2), and second, the evolution of permeability with effective stress is assumed to be a quadratic polynomial function. Therefore, the essence of the response surface method is to use a bivariate quadratic polynomial to fit the experimental data and then obtain statistical parameters and effective stress coefficients.
The aforementioned analysis showed that the calculation of the effective stress coefficient by the cross-plotting method and response surface method had the following defects: (1) The specific functional relationship between permeability and effective stress could not satisfy the experimental data for all rock types.
(2) Preprocessing of the original experimental data may lead to fitting errors and even distortion of the analysis results.

3D Surface Fitting Method
On the basis of defects of the cross-plotting method and response surface method, this study proposes a method of fitting experimental data with a 3D surface to obtain the effective stress coefficient. The functional relationship between rock permeability and effective stress could be expressed as follows: Equation (10) mainly contains three types: (1) Power law model [23,24] (2) Exponential model [25] (3) Quadratic polynomial model [26] where a, b, c, and d are fitting parameters, and a is the ESCP of each model. An assumption was that the effective stress linearly correlated with the confining and pore fluid pressures, satisfying Equation (2). Assuming that Equation (10) is continuously differentiable within the range of the stress variation, The ESCP represents the relative influence of pore fluid and confining pressures on permeability, which was used to evaluate the contribution of two independent variables (σ, P) to a dependent variable (k). The ESCP can be expressed as the ratio of the partial derivative coefficients of the two independent variables: The calculation of the ESCP was transformed into the problem of fitting the 3D surface of experimental data on the confining pressure, pore fluid pressure, and permeability. The Curve Fitting toolbox in MATLAB was used to fit the experimental data through three models, and three fitting parameters, a, were obtained. Finally, the a in the optimal fitting scheme was selected as the ESCP of the rock sample.
The 3D surface fitting method is proposed according to linear effective stress law and continuity hypothesis. If the effective stress cannot be expressed by the linear relationship between confining pressure and pore pressure, the 3D surface fitting method is no longer applicable.
Two tight sandstones from British Columbia, Canada, were Jurassic sandstones. The samples were dark gray and mainly consisted of quartz, black silica, and feldspar clasts. The porosity of the sandstone ranged from 1.65% to 3.23%. The rock samples were cylindrical specimens with a diameter of 38 mm and a height of 51 mm. The experimental data for rock samples 1 and 2 are listed in Table 1. The 3D surface fitting method was used to fit the experimental data of the two rock samples (Table 1). Table A1 in the Appendix A lists the precision and calculation results of different fitting functions, and Figure 1 shows the corresponding fitting curves.  The 3D surface fitting method was used to fit the experimental data of the two rock samples (Table 1). Table A1 in the Appendix A lists the precision and calculation results of different fitting functions, and Figure 1 shows the corresponding fitting curves.  Two lithic sandstones were taken from the Permian and buried at a depth of approximately 2700 m. The porosity of D13−4 was 6.48%, and that of D15−2 was 6.88%. Experimental data and calculation results of the lithic sandstones are shown in Tables A2 and A3 in the Appendix A. The 3D surface fitting method was used to fit the experimental data of the two rock samples, D13−4 and D15−2. Table A4 in the Appendix A lists the precision and calculation results of the different fitting functions. The ESCP of rock samples D13−4 and D15−2 was set as 1, three functional relationships were used to fit the experimental data; the fitting correlation coefficients were all lower than that of the 3D surface fitting method. Two lithic sandstones were taken from the Permian and buried at a depth of approximately 2700 m. The porosity of D13−4 was 6.48%, and that of D15−2 was 6.88%. Experimental data and calculation results of the lithic sandstones are shown in Tables A2 and A3 in the Appendix A.
The 3D surface fitting method was used to fit the experimental data of the two rock samples, D13−4 and D15−2. Table A4 in the Appendix A lists the precision and calculation results of the different fitting functions. The ESCP of rock samples D13−4 and D15−2 was set as 1, three functional relationships were used to fit the experimental data; the fitting correlation coefficients were all lower than that of the 3D surface fitting method.  Figure 2 shows that the fitting correlation coefficient is only 0.7096 when the ESCP is set as 1. It cannot accurately describe the variation law of sandstone permeability with confining and pore fluid pressures if the effective stress coefficient is set as 1. Therefore, the first step in the evaluation of stress sensitivity of reservoir permeability involves obtaining the ESCP based on the experimental data, and the second step involves the evaluation of reservoir stress sensitivity.   Figure 3a,b shows that when the permeability is between 0.03 and 0.05 mD, the 3D surface fitting method has a better effect than the response surface method. When the effective stress coefficient was set to 0.6479, most data points were not on the fitting surface. The 3D surface fitting method was used to fit the experimental data of the two rock samples, D13−4 and D15−2. Table A4 in the Appendix A lists the precision and calculation results of the different fitting functions. The ESCP of rock samples D13−4 and D15−2 was set as 1, three functional relationships were used to fit the experimental data; the fitting correlation coefficients were all lower than that of the 3D surface fitting method.  Figure 2 shows that the fitting correlation coefficient is only 0.7096 when the ESCP is set as 1. It cannot accurately describe the variation law of sandstone permeability with confining and pore fluid pressures if the effective stress coefficient is set as 1. Therefore, the first step in the evaluation of stress sensitivity of reservoir permeability involves obtaining the ESCP based on the experimental data, and the second step involves the evaluation of reservoir stress sensitivity. Figure 3a,b shows that when the permeability is between 0.03 and 0.05 mD, the 3D surface fitting method has a better effect than the response surface method. When the effective stress coefficient was set to 0.6479, most data points were not on the fitting surface.  Figure 2 shows that the fitting correlation coefficient is only 0.7096 when the ESCP is set as 1. It cannot accurately describe the variation law of sandstone permeability with confining and pore fluid pressures if the effective stress coefficient is set as 1. Therefore, the first step in the evaluation of stress sensitivity of reservoir permeability involves obtaining the ESCP based on the experimental data, and the second step involves the evaluation of reservoir stress sensitivity. Figure 3a,b shows that when the permeability is between 0.03 and 0.05 mD, the 3D surface fitting method has a better effect than the response surface method. When the effective stress coefficient was set to 0.6479, most data points were not on the fitting surface. Figures 1-3 demonstrate that the fitting effect of the 3D surface fitting method is better than the cross-plotting method and response surface method. The following section presents a quantitative analysis of the various methods from the perspective of data statistics.

Discussion
The accuracy and reliability of the 3D surface fitting method, cross-plotting method, and response surface method were quantitatively analyzed by selecting three common statistical parameters: correlation coefficient, root mean square error (RMSE), and residual error.

Comparison of the 3D Surface Fitting Method and Cross-Plotting Method
The ESCP values of rock samples 1 and 2 obtained by the 3D surface fitting method are 1.551 and 1.539, respectively. The results obtained by the different fitting types are relatively stable and consistent (Table A1 in the Appendix A). The ESCPs of rock samples 1 and 2 obtained by the cross-plotting method were 0.509 and 0.612, respectively. An ESCP of less than 1 indicates that the pore fluid pressure has less effect on permeability than the confining pressure; more than 1 means the opposite. The results obtained by the two methods are the opposite. The original experimental data must be analyzed to determine which method is the most reliable.
Three groups of experimental data were selected ( Table 2). Using the permeability under a confining pressure of 30 MPa and pore fluid pressure of 5 MPa as a reference, the permeability change values were obtained when the confining pressure was reduced by 20 MPa, and the pore fluid pressure was increased by 20 MPa, respectively. The changing values of the confining and pore fluid pressures were the same. If the change in permeability during the period of change in the pore fluid pressure is greater than that during confining pressure change, the pore fluid pressure has a greater influence on permeability, and the ESCP should be greater than 1. The calculation results in Table 2 show that when the pore fluid pressure changes by 20 MPa, the permeability change values of rock samples 1 and 2 are 0.02583 and 0.0289 mD, respectively. When the confining pressure changes by 20 MPa, the permeability change values of rock samples 1 and 2 are 0.00643 and 0.0118 mD, respectively. The effect of pore fluid pressure on permeability is significantly greater than that of confining pressure; that is, the ESCPs of rock samples 1 and 2 are greater than 1.
The ESCP calculated using the cross-plotting method is inconsistent with the experimental data. The main reason for this deviation is that if the original data points cannot coincide with the cross plot, the confining and pore fluid pressures of the cross point must be calculated using the interpolation method. Bernabe [13] pointed out that the stress-permeability data were too scattered and that interpolation could lead to statistical errors; thus, interpolation was generally applied to stress-strain analysis.

Comparison of Three Fitting Functions
Three functions were used to fit the experimental data of rock samples 1 and 2. The quadratic polynomial type had the highest correlation coefficient, followed by the exponential type, and the power law type had a negative correlation coefficient. The relationship between permeability and effective stress in rock samples 1 and 2 did not agree with the power law function.
The quadratic polynomial model that Yin and Wang [26] proposed was based on the thick-walled tube theory, mainly used to describe the quadratic parabolic relationship between permeability and effective stress of capillary rock. The matrix rock and fractured rock can be described by exponential and power law models, respectively. Zhao [28] and Dong have posited that the logarithmic, quadratic polynomial, and power law models are the best for characterizing fractured rock, porous tight sandstone, and fine-grained sandstone, respectively. The power law model, introduced by Shi [24], was using to fit experimental data and achieved good results. However, Xiao [18] and Dong [23] have pointed out that the power law relationship is an empirical formula based on the fitting of experimental data, without clear physical meaning and theoretical support. The fitting results in this study also showed that the power law relation was invalid and not universally applicable.
David et al. [25] discussed the exponential model in detail and posited that the exponential relation is an empirical formula based on experimental data. When the experimental pressure was less than the critical pressure, the equation was in good agreement with the experimental data. Notably, the exponential model can be simplified using Equation (18).
McKee et al. [29] deduced the following formulas for the functional relationships between the effective stress and two factors, porosity and permeability: where φ 0 and k 0 are the porosity and permeability when the effective stress is 0, c p is the average pore compressibility, and ∆σ is the change in the effective stress. Because the initial porosity of tight sandstone is small, Equation (18) can be simplified as follows: Equation (19) is an exponential model that is simple in mathematical form, easy to use, and widely used to fit stress-sensitive experimental data. Table A5 in the Appendix A lists the experimental research information on stress sensitivity. Many experimental results have demonstrated that this simplification was acceptable for rocks with small initial porosity.

Correlation Coefficient and RMSE
The ESCP obtained by the cross-plotting and 3D surface fitting methods was used to fit the experimental data. The correlation coefficient (R 2 ) and RMSE are listed in Table 3. In Table 3, the correlation coefficients of the 3D surface fitting method are all greater than 0.9, and the correlation coefficients of the cross-plotting method are less than 0.3. The larger the correlation coefficients, the closer the statistical results are to the actual situation. The RMSE of the 3D surface fitting method is approximately 0.003, and that of the crossplotting method is approximately 0.009, which is three times that of the former. The bigger the RMSE, the bigger the difference between the statistical results and the actual situation is. Comparing the correlation coefficient and RMSE of the two methods demonstrates the advantages of the 3D surface fitting method.

Residual Error
To directly reflect the difference between the cross-plotting and 3D surface fitting methods, distribution maps of the residual error for the two methods under different fitting relations were drawn (Figures 4 and 5). Figure 4a,b shows that the residual error of the cross-plotting method ranges from −0.012 to 0.018, and the maximum residual error is 0.01784. Figure 4c,d shows that the distribution interval of the residual error of the 3D surface fitting method is from −0.004 to 0.0068, and the maximum value is 0.00672. The scope and maximum value of the residual error of the 3D surface fitting method are smaller than those of the cross-plotting method by nearly three times. Figure 5a,b shows that the residual error of the cross-plotting method ranges from −0.012 to 0.018, and the maximum residual error is 0.01736. Figure 5c,d shows that the distribution interval of the residual error of the 3D surface fitting method is from −0.004 to 0.0068, and the maximum value is 0.00672. The scope and maximum value of the residual error of the 3D surface fitting method are nearly three times smaller than those of the cross-plotting method, indicating that the accuracy and reliability of the 3D surface fitting method are much higher than those of the cross-plotting method.
Comparative analysis of the three statistical parameters shows that the ESCP obtained by the cross-plotting method cannot accurately fit the experimental data of rock, and the correlation coefficients are less than 0.3. The correlation coefficients of the 3D surface fitting method were all greater than 0.9, demonstrating the method's benefit. In Table 3, the correlation coefficients of the 3D surface fitting method are all greater than 0.9, and the correlation coefficients of the cross-plotting method are less than 0.3. The larger the correlation coefficients, the closer the statistical results are to the actual situation. The RMSE of the 3D surface fitting method is approximately 0.003, and that of the cross-plotting method is approximately 0.009, which is three times that of the former. The bigger the RMSE, the bigger the difference between the statistical results and the actual situation is. Comparing the correlation coefficient and RMSE of the two methods demonstrates the advantages of the 3D surface fitting method.

Residual Error
To directly reflect the difference between the cross-plotting and 3D surface fitting methods, distribution maps of the residual error for the two methods under different fitting relations were drawn (Figures 4 and 5).    Figure 4c,d shows that the distribution interval of the residual error of the 3D surface fitting method is from −0.004 to 0.0068, and the maximum value is 0.00672. The scope and maximum value of the residual error of the 3D surface fitting method are smaller than those of the cross-plotting method by nearly three times. Figure 5a,b shows that the residual error of the cross-plotting method ranges from −0.012 to 0.018, and the maximum residual error is 0.01736. Figure 5c,d shows that the distribution interval of the residual error of the 3D surface fitting method is from −0.004 to 0.0068, and the maximum value is 0.00672. The scope and maximum value of the residual error of the 3D surface fitting method are nearly three times smaller than those of the cross-plotting method, indicating that the accuracy and reliability of the 3D surface fitting method are much higher than those of the cross-plotting method.
Comparative analysis of the three statistical parameters shows that the ESCP obtained by the cross-plotting method cannot accurately fit the experimental data of rock, and the correlation coefficients are less than 0.3. The correlation coefficients of the 3D surface fitting method were all greater than 0.9, demonstrating the method's benefit. Table A4 in the Appendix A lists the precision and calculation results of the 3D surface fitting method using different fitting functions. The fitting accuracy of the power law type is the highest, but the calculation results of the effective stress coefficient significantly differ from those of the other two types. The exponential and quadratic polynomial models have the theoretical basis support, the difference in calculation result is small. The power law relation had no theoretical basis and failed to fit rock samples 1 and 2. Therefore, the exponential and quadratic polynomial models are more reliable than the power  Table A4 in the Appendix A lists the precision and calculation results of the 3D surface fitting method using different fitting functions. The fitting accuracy of the power law type is the highest, but the calculation results of the effective stress coefficient significantly differ from those of the other two types. The exponential and quadratic polynomial models have the theoretical basis support, the difference in calculation result is small. The power law relation had no theoretical basis and failed to fit rock samples 1 and 2. Therefore, the exponential and quadratic polynomial models are more reliable than the power law relation. In this study, the power law relation was not considered in the calculation of ESCP.

Comparison of the 3D Surface Fitting Method and Response Surface Method
The ESCP values of rock samples D13−4 and D15−2 obtained by the 3D surface fitting method were 0.3537 and 0.7767, respectively. The ESCP values of rock samples D13−4 and D15−2 obtained by the response surface method were 0.2592 and 0.6497 [18], respectively.

Correlation Coefficient and RMSE
The ESCPs obtained by the response surface and 3D surface fitting method were used to fit the test data of samples D13−4 and D15−2. The correlation coefficient and RMSE are listed in Table 4.
For rock sample D13−4, whether a quadratic polynomial model or exponential model is used for fitting, the correlation coefficient of the 3D surface fitting method is larger than that of the response surface method, while the RMSE is smaller than that of the response surface method. The same situation is observed in rock sample D15−2.
The comparison of the correlation coefficient and RMSE of the two methods showed a slightly improved data processing accuracy when using the 3D surface fitting method.

Residual Error
Experimental data of rock samples D13−4 and D15−2 were processed to draw residual error distribution maps of the response surface and 3D surface fitting methods by using different fitting functions (Figures 6 and 7).
In Figure 6a,b, the distribution range of the residual error of the response surface method is from −0.04 to 0.07, and the maximum value of the residual error is 0.06804. Figure 6c,d shows that the distribution interval of the residual error of the 3D surface fitting method is from −0.04 to 0.06, and the maximum value of the residual error is 0.05842. The distribution interval and maximum value of the residual error of the latter are approximately 0.01 smaller than those of the former, indicating that the fitting accuracy of the latter is slightly higher than that of the former.
In Figure 7a,b, the distribution range of residual error of the response surface method is from −0.014 to 0.01, and the maximum value of the residual error is 0.01391. Figure 7c,d shows that the distribution interval of the residual error of the 3D surface fitting method is from −0.0107 to 0.01, and the maximum value of the residual error is 0.0107. The distribution interval and maximum value of the residual error of the latter are approximately 0.003 smaller than those of the former, which also indicates that the fitting accuracy of the latter is higher than that of the former.  For rock sample D13−4, whether a quadratic polynomial model or exponential model is used for fitting, the correlation coefficient of the 3D surface fitting method is larger than that of the response surface method, while the RMSE is smaller than that of the response surface method. The same situation is observed in rock sample D15−2.
The comparison of the correlation coefficient and RMSE of the two methods showed a slightly improved data processing accuracy when using the 3D surface fitting method.

Residual Error
Experimental data of rock samples D13−4 and D15−2 were processed to draw residual error distribution maps of the response surface and 3D surface fitting methods by using different fitting functions (Figures 6 and 7).  In Figure 6a,b, the distribution range of the residual error of the response surface method is from −0.04 to 0.07, and the maximum value of the residual error is 0.06804. Figure 6c,d shows that the distribution interval of the residual error of the 3D surface fitting method is from −0.04 to 0.06, and the maximum value of the residual error is 0.05842. The distribution interval and maximum value of the residual error of the latter are approximately 0.01 smaller than those of the former, indicating that the fitting accuracy of the latter is slightly higher than that of the former.
In Figure 7a,b, the distribution range of residual error of the response surface method is from −0.014 to 0.01, and the maximum value of the residual error is 0.01391. Figure 7c,d shows that the distribution interval of the residual error of the 3D surface fitting method is from −0.0107 to 0.01, and the maximum value of the residual error is 0.0107. The distribution interval and maximum value of the residual error of the latter are approximately 0.003 smaller than those of the former, which also indicates that the fitting accuracy of the latter is higher than that of the former.
The 3D surface fitting method performed better than the response surface method in three statistical parameters; although the advantages were not obvious, the superiority was assured. The 3D surface fitting method is of great significance for improving the processing accuracy of experimental data, and lays a good foundation for the research of effective stress law, permeability evolution law, and evaluation of reservoir stress sensitivity.

Conclusions
This study analyzed the essence of the cross-plotting and response surface methods. The analysis demonstrated that the preprocessing of the original experimental data and the preset of the function relation might lead to the calculation deviation of ESCP. To avoid calculation deviation and improve the fitting accuracy, a 3D surface fitting method was proposed to calculate the ESCP. The 3D surface fitting method was compared with the cross-plotting and response surface methods by quantitative analysis of the statistical parameters. The 3D surface fitting method performed better than the response surface method in three statistical parameters; although the advantages were not obvious, the superiority was assured. The 3D surface fitting method is of great significance for improving the processing accuracy of experimental data, and lays a good foundation for the research of effective stress law, permeability evolution law, and evaluation of reservoir stress sensitivity.

Conclusions
This study analyzed the essence of the cross-plotting and response surface methods. The analysis demonstrated that the preprocessing of the original experimental data and the preset of the function relation might lead to the calculation deviation of ESCP. To avoid calculation deviation and improve the fitting accuracy, a 3D surface fitting method was proposed to calculate the ESCP. The 3D surface fitting method was compared with the cross-plotting and response surface methods by quantitative analysis of the statistical parameters.
(1) Compared with the cross-plotting method, the 3D surface fitting method had obvious advantages: the correlation coefficient was increased by 0.6, and the RMSE and maximum residual error were reduced by nearly three times. The preprocessing of the original experimental data led to the distortion of the calculation results by the crossplotting method.
(2) The three statistical parameters of the 3D surface fitting method were better than those of the response surface method, and the calculation results were more accurate and reliable.
(3) The power law relationship model between permeability and stress failed to fit the stress-sensitive experimental data of some rock samples, and no basic theoretical support was observed. Therefore, a power law model should be used with caution.
(4) The exponential model of permeability change with stress can be simplified by a theoretical formula, and its mathematical form is simple and widely used. The key to investigating the effective stress and fluid-structure interaction of porous media lied in the mechanism and functional relationship of the confining pressure and pore fluid pressure on the flow characteristics [30]. When the mechanism of action and function relation were not determined, only the effective stress was insufficient, especially when the effective stress could not be expressed by a simple linear relation.
Author Contributions: All authors have contributed to this article. Conceptualization, methodology, data curation, writing-original draft, X.Z.; data curation, J.S.; and writing-review and editing, J.L. All authors have read and agreed to the published version of the manuscript.