An Analytical Model for Capturing the Decline of Fracture Conductivity in the Tuscaloosa Marine Shale Trend from Production Data

Fracture conductivity decline is a concern in the Tuscaloosa Marine Shale (TMS) wells due to the high content of clay in the shale. An analytical well productivity model was developed in this study considering the pressure-dependent conductivity of hydraulic fractures. The log-log diagnostic approach was used to identify the boundary-dominated flow regime rather than the linear flow regime. Case studies of seven TMS wells indicated that the proposed model allows approximation of the field data with good accuracy. Production data analyses with the model revealed that the pressure-dependent fracture conductivity in the TMS in the Mississippi section declines following a logarithmic mode, with dimensionless coefficient χ varying between 0.116 and 0.130. The pressure-dependent decline of fracture conductivity in the transient flow period is more significant than that in the boundary-dominated flow period.


Introduction
The Tuscaloosa Marine Shale (TMS) across Louisiana and Mississippi has been an attractive unconventional shale oil reservoir since 2012 [1]. The TMS is a sedimentary formation that consists of organic-rich fine-grained sediments deposited during the Upper Cretaceous [2]. The TMS is one part of the Tuscaloosa group consisting of Upper Tuscaloosa, Middle Marine Shale, and Lower Tuscaloosa [3]. The thickness of TMS ranges from 500 ft in the southwestern Mississippi to more than 800 ft in southeastern Louisiana, within a depth range of 11,000 to more than 15,000 ft [4]. Middle Tuscaloosa is composed of a dark grey, fissile, and sandy marine shale. Experimental results indicated that the permeability ranges from less than 0.01 to 0.06 md, and porosity ranges from 2.3% to 8.0% [4].
More than 80 multi-fractured horizontal wells were drilled in the TMS between 2012 and 2014 [2]. Several TMS wells were recorded to have appealing initial oil production rates exceeding 1000 stb/d. Although such significant high production rates were observed, drilling activities stopped in 2014 due to the high cost of drilling and low price of oil. Previous studies estimated an unproven recoverable oil of seven billion barrels [4], but the real production potential from TMS is still poorly understood. It is of significance to assess the TMS well performance through production decline analysis and identify key factors controlling the productivity of TMS wells.
The oil formation is isotropic.

2.
Boundary-dominated flow has been reached within the fracture drainage area.

3.
Linear flow prevails from the shale matrix to the fractures. 4.
Fracture and formation damages are negligible.

5.
No change in fluid composition during production. 6.
Hydraulic fractures have the same geometry. 7.
Reservoir pressure is above the bubble point pressure.
Derivation of the analytical model considering the pressure-dependent decline of fracture conductivity is shown in Appendix A. The analytical model for predicting the productivity of multi-fractured shale oil wells is, where q o is the oil production rate, n f is the number of hydraulic fractures, k m is the matrix permeability, h is the reservoir thickness, p i is the initial formation pressure, p w is the wellbore pressure, B o is the oil formation volume factor, µ o is the oil viscosity, S f is the hydraulic fracture spacing, x f is the fracture half-length, α b is the Biot coefficient, ν is the Poisson's ratio, N p is the cumulative oil production, N i is the original oil in place within the well drainage area, c t is the total compressibility, C f0 is the initial fracture conductivity, and c p is the compressibility of the proppant pack. The analytical model for capturing time-dependent fracture conductivity of shale wells during the boundary-dominate flow period is, where C f is the fracture conductivity. From Equation (4) we can see that fracture conductivity can be predicted if the cumulative oil production data are known.

Flow Regime Diagnosis Method
As stated earlier, the proposed model is only applicable to the wells that have reached the boundary-dominated flow over the production time. Therefore, preliminary analysis of production data should be performed to make sure that the candidate wells are in the boundary-dominated flow regime.
Four flow regimes may exist in multi-fractured reservoirs: (1) early time reservoir linear flow (transient flow), (2) mid-time boundary-dominated flow, (3) late time reservoir radial flow, and (4) very late time bounded flow. For wells in shale plays the third and fourth flow regimes usually are not seen due to the ultralow permeability of shale matrix. To fully understand the performance of multi-fractured horizontal wells, we identify the first two flow regimes between fractures. Figure 1a presents two fractures with an assumption of a virtual boundary between these two fractures. In the transient flow period, pressure propagates outward from the fracture face without encountering the virtual boundary. In the boundary-dominated flow period, the pressure transient has reached the virtual boundary, and the static pressure is declining at the boundary, as shown in Figure 1b. In this section, the log-log diagnostic plots of production rate versus real production time as well as material balance time (MBT) are constructed to identify the flow regime of TMS wells. The procedure for identifying the flow regime is as follows [30].
Construct the production rate versus real production time curve (i.e., PT curve) and production In this section, the log-log diagnostic plots of production rate versus real production time as well as material balance time (MBT) are constructed to identify the flow regime of TMS wells. The procedure for identifying the flow regime is as follows [30].
Construct the production rate versus real production time curve (i.e., PT curve) and production rate versus material balance curve (i.e., MBT curve) on the same log-log plot. Then plot a tangent line L MBT using the end of MBT curve. If the slope k MBT of the tangent line L MBT is close to −1, then we have confidence that the boundary-dominated flow regime has been reached. If the slope k MBT is greater than -1, one can move the line L MBT onto the PT curve and get an approximately tangent point P. We draw another tangent line L PT using data between point P and the end of production data. If the slope k PT of this line is less than or equal to −1, it is believed that the well has reached the boundary-dominated regime, and the real production time corresponding to point P is the estimated switching time from transient flow to the boundary-dominated flow. However, if it is difficult to draw a tangent line L PT with slope k PT less than or equal to −1 owing to few points after point P, then one failed to obtain the conclusion that the well has reached boundary-dominated flow.

TMS Well Description
The study area in Mississippi is presented in Figure 2. The shapes for Louisiana and Mississippi counties were downloaded from the website (https://www2.census.gov). The production data as of May 2018 were gathered from 55 TMS wells in Mississippi from the Alfred C. Moore, Pearl River, and Henry Fields. These data were downloaded from the Mississippi Automated Resource Information System (MARIS) website (www.maris.state.ms.us).
Quality control of the production data was performed to validate the compliance of the field data with the assumptions used in the analysis. Those wells that had abrupt changes in monthly oil production rate during the production history were removed. This helps minimize the influence of changes in the production operations that would have a great effect on the production history. Following this, seven TMS wells that had reached the boundary-dominated flow regime were analyzed. Other wells could not be analyzed because of absence of a definite trend in the reported production data, most likely dictated by the used production strategy for which no details are available.
It should be stated that other wells may also have reached the boundary-dominated flow regime. However, we did not analyze these wells due to lack of completion data and also abrupt changes in production data of these wells over time. Seven TMS wells used in this study are shown as blue dots in Figure 2. The wells' labels and their locations are listed in Table 1. Proppant has a major impact on fracture conductivity. Unfortunately, the type of proppant used in these seven TMS wells could not be identified due to the lack of completion data. The average fracture spacing can be estimated by using the stages.   Figure 3 shows the log-log diagnostic plot of the seven TMS wells. Well 1 has a horizontal length of 9102 ft and an effective lateral length of 8442 ft with 29 stages, and each stage has four clusters. We drew a tangent line of MBT curve and found that the production data at the end of material balance time lies on the tangent line LMBT with a slope of −0.779. As kMBT is less than −1, the tangent line, LMBT, was moved onto the PT curve, and it showed tangency to this curve. If we draw a tangent line LPT using production data vs. real production time. The slope kPT is found to be −1.085, indicating a typical behavior of boundary-dominated flow over 16 months of production.

Flow Regimes of TMS Wells
Well 2 was drilled and completed with an effective lateral length of 6757 ft and a 22 stage hydraulic fracturing operation. We can see from Figure 3(b) that the slopes of MBT curve (kMBT) and PT curve (kPT) are −0.688 and −1.002, respectively, indicating that well 2 should be in the boundarydominated flow regime with an estimated switching time of 9.5 months.
Well 3 was drilled and completed with a 4,508 ft lateral and a ten stage hydraulic fracturing operation. From this plot, we can see that kMBT and kPT are −0.845 and −1.295, respectively. Therefore, we have confidence that well 3 switched from the transient flow to the boundary-dominated flow over 27 months of production.   Figure 3 shows the log-log diagnostic plot of the seven TMS wells. Well 1 has a horizontal length of 9102 ft and an effective lateral length of 8442 ft with 29 stages, and each stage has four clusters. We drew a tangent line of MBT curve and found that the production data at the end of material balance time lies on the tangent line L MBT with a slope of −0.779. As k MBT is less than −1, the tangent line, L MBT, was moved onto the PT curve, and it showed tangency to this curve. If we draw a tangent line L PT using production data vs. real production time. The slope k PT is found to be −1.085, indicating a typical behavior of boundary-dominated flow over 16 months of production.

Flow Regimes of TMS Wells
Well 2 was drilled and completed with an effective lateral length of 6757 ft and a 22 stage hydraulic fracturing operation. We can see from Figure 3b that the slopes of MBT curve (k MBT ) and PT curve (k PT ) are −0.688 and −1.002, respectively, indicating that well 2 should be in the boundary-dominated flow regime with an estimated switching time of 9.5 months.
Well 3 was drilled and completed with a 4,508 ft lateral and a ten stage hydraulic fracturing operation. From this plot, we can see that k MBT and k PT are −0.845 and −1.295, respectively. Therefore, we have confidence that well 3 switched from the transient flow to the boundary-dominated flow over 27 months of production. Figure  The true vertical depth and effective lateral length of well 7 are 11,841 and 5681 ft, respectively. Slope k MBT and k PT are found to be −0.714 and −0.988, respectively. As k PT is approximately close to −1, well 7 should be in the boundary-dominated flow regime in over 13 months from its first production month. The true vertical depth and effective lateral length of well 7 are 11,841 and 5681 ft, respectively. Slope kMBT and kPT are found to be −0.714 and −0.988, respectively. As kPT is approximately close to −1, well 7 should be in the boundary-dominated flow regime in over 13 months from its first production month.

Verification of the Model on Simulated Production Dataset
To verify the proposed model, we compared the model-predicted production rates with simulation results calculated by COMSOL Multiphysics. Figure 4 shows the grid of the model. Fluid

Verification of the Model on Simulated Production Dataset
To verify the proposed model, we compared the model-predicted production rates with simulation results calculated by COMSOL Multiphysics. Figure 4 shows the grid of the model. Fluid flow in the half fracture (x f ) and half spacing (0.5S f ) area is modeled for convenience. In this case, Darcy's law in the subsurface flow module was used to model fluid flow within the porous medium. The hydraulic fracture was represented by a boundary layer with a thickness of 0.5 w (i.e., line AB). The decline of fracture conductivity is modeled using equation (A3). The following assumptions were made in this section: the fracture half-length was 300 ft, the fracture spacing was 60 ft, the initial reservoir pressure was 6,200 psia, the bottom hole pressure was 4,000 psia (at point A), the fracture height was 100 ft, the matrix permeability was 0.0002 md, the number of clusters was 120, the porosity was 8%, the oil formation factor was 1.5 rb/stb, the oil viscosity was 0.5 cp. We compared six cases, case 1: C f0 = 2 md-ft, and c p = 0.0005 psia −1 ; case 2: C f0 = 2 md-ft, and c p = 0.001 psia −1 ; case 3: C f0 = 5 md-ft, and c p = 0.0005 psia −1 ; case 4: C f0 = 5 md-ft, and c p = 0.001 psia −1 ; case 5: C f0 = 10 md-ft, and c p = 0.0005 psia −1 ; case 6: C f0 = 10 md-ft, and c p = 0.001 psia −1 .

Verification of the Model on Simulated Production Dataset
To verify the proposed model, we compared the model-predicted production rates with simulation results calculated by COMSOL Multiphysics. Figure 4 shows the grid of the model. Fluid flow in the half fracture (xf) and half spacing (0.5Sf) area is modeled for convenience. In this case, Darcy's law in the subsurface flow module was used to model fluid flow within the porous medium. The hydraulic fracture was represented by a boundary layer with a thickness of 0.5 w (i.e., line AB). The decline of fracture conductivity is modeled using equation (A3). The following assumptions were made in this section: the fracture half-length was 300 ft, the fracture spacing was 60 ft, the initial reservoir pressure was 6,200 psia, the bottom hole pressure was 4,000 psia (at point A), the fracture height was 100 ft, the matrix permeability was 0.0002 md, the number of clusters was 120, the porosity was 8%, the oil formation factor was 1.5 rb/stb, the oil viscosity was 0.5 cp. We compared six cases,  If the predicted production rates match the simulation results, the variation of fracture conductivity can be captured using the proposed analytical model based on the production data under boundary-dominated flow conditions. R-square (R 2 ) and average absolute relative error percentage (AAREP) were utilized to estimate the proposed model. AAREP is defined as [31,32], where N is the number of data points; Oi is the observed data; Ei is the predicted data. If the predicted production rates match the simulation results, the variation of fracture conductivity can be captured using the proposed analytical model based on the production data under boundary-dominated flow conditions. R-square (R 2 ) and average absolute relative error percentage (AAREP) were utilized to estimate the proposed model. AAREP is defined as [31,32], where N is the number of data points; O i is the observed data; E i is the predicted data. The fracture compressibility can be obtained by fitting the proposed model to the production rate calculated by COMSOL. Figure 5 compares the oil production rate calculated by the proposed model and that of COMSOL Multiphysics software. After 40 months of production, the fracture conductivity was reduced by 53% and 69% for case 1 and case 2, respectively. The comparison between the prescribed fracture conductivity and the predicted value is presented in Table 2. Besides, there is good agreement between the prescribed fracture compressibility and the model-predicted value, with relative differences of fracture compressibility of less than 10%. The fracture compressibility can be obtained by fitting the proposed model to the production rate calculated by COMSOL. Figure 5 compares the oil production rate calculated by the proposed model and that of COMSOL Multiphysics software. After 40 months of production, the fracture conductivity was reduced by 53% and 69% for case 1 and case 2, respectively. The comparison between the prescribed fracture conductivity and the predicted value is presented in Table 2. Besides, there is good agreement between the prescribed fracture compressibility and the model-predicted value, with relative differences of fracture compressibility of less than 10%.  Figure 5. Comparison between the model-predicted production rate and the simulation result. Figure 6 shows the comparison of field data and predicted values for the seven TMS wells studied in this paper. The calculated R 2 and AAREP values are listed in Table 3. For example, well 1 produced 1072 barrels of oil equivalent per day during the first 30 days of production. Even though this well produced considerable amounts of hydrocarbons, its production rate decreased to 120 stb/d over 16 months of production. After 16 months, the cumulative oil production was 148,762 bbl. The R 2 and AAREP are 0.925 and 10.38%, respectively. In general, relatively high R 2 and low AAREP values were observed, indicating a good agreement between the field data and predicted results.  Figure 6 shows the comparison of field data and predicted values for the seven TMS wells studied in this paper. The calculated R 2 and AAREP values are listed in Table 3. For example, well 1 produced 1072 barrels of oil equivalent per day during the first 30 days of production. Even though this well produced considerable amounts of hydrocarbons, its production rate decreased to 120 stb/d over 16 months of production. After 16 months, the cumulative oil production was 148,762 bbl. The R 2 and AAREP are 0.925 and 10.38%, respectively. In general, relatively high R 2 and low AAREP values were observed, indicating a good agreement between the field data and predicted results.     Figure 7 presents the variation of the predicted fracture conductivity with time for the seven TMS wells. The fitted fracture compressibility is listed in Table 3. The numerical method used to calculate the fracture conductivity is demonstrated in Appendix A. The fracture conductivity changes over time owing to depletion of pressure in the hydraulic fractures. Previous studies found that time-dependent fracture conductivity changes in logarithmic or exponential types [28]. It is evident from Figure 7 that the fracture conductivity for these TMS wells varies logarithmically with the production time. Take well 1 as an example. The fracture conductivity at the switching time is about 0.63 times its initial value, which gives an average decline rate of 27.3% per year in the transient flow period. The fracture conductivity was reduced by 52% over 65 months of production. We can see that the variation of fracture conductivity in the transient flow period is more significant than that in the boundary-dominated flow period.   Figure 7 presents the variation of the predicted fracture conductivity with time for the seven TMS wells. The fitted fracture compressibility is listed in Table 3. The numerical method used to calculate the fracture conductivity is demonstrated in Appendix. The fracture conductivity changes over time owing to depletion of pressure in the hydraulic fractures. Previous studies found that timedependent fracture conductivity changes in logarithmic or exponential types [28]. It is evident from Figure 7 that the fracture conductivity for these TMS wells varies logarithmically with the production time. Take well 1 as an example. The fracture conductivity at the switching time is about 0.63 times its initial value, which gives an average decline rate of 27.3% per year in the transient flow period. The fracture conductivity was reduced by 52% over 65 months of production. We can see that the variation of fracture conductivity in the transient flow period is more significant than that in the boundary-dominated flow period.  We also compared the time-dependent fracture conductivity with the empirical model expressed by [28],

Variation of Fracture Conductivity
where t is the production time in months, χ is the dimensionless coefficient describing the time-dependent fracture conductivity. Figure 8 shows the comparison of fracture conductivity of well 4 between our proposed model and the empirical equation. A summary of R 2 and AAREP values for the seven TMS wells is listed in Table 4. From Figure 8 we can see that the empirical equation exhibits a linear trend on the semi-log plot. Our proposed model predicts an S-type decline curve. The dimensionless coefficient χ has been found in the range of 0.116~0.130 with an average of 0.123 for the seven wells. Very high values of R 2 and extremely low AAREP values have been observed, indicating excellent goodness of fit between our proposed model and the empirical equation. and the empirical equation. A summary of R 2 and AAREP values for the seven TMS wells is listed in Table 4. From Figure 8 we can see that the empirical equation exhibits a linear trend on the semi-log plot. Our proposed model predicts an S-type decline curve. The dimensionless coefficient χ has been found in the range of 0.116~0.130 with an average of 0.123 for the seven wells. Very high values of R 2 and extremely low AAREP values have been observed, indicating excellent goodness of fit between our proposed model and the empirical equation.

Discussion
A mathematical model was developed in this study to capture the pressure-dependent decline of fracture conductivity from production data. As stated previously, the newly proposed model is suitable for wells with production in the boundary-dominated flow period. It has been recognized that the production rate declines with material balance time with slopes of −1/2 and −1 during the linear flow and boundary-dominated flow periods, respectively, in the log-log plot when there is no skin effect [33]. Many multi-fractured horizontal wells producing from shale reservoirs have been reported to experience a long-term linear flow (transient flow) period. The presence of decline of fracture conductivity and formation damage is usually interpreted as skin factor, which may change the shape of the diagnostic plot [22]. Nobakht and Mattar [34] stated that the linear flow is no longer a straight line with a slope of −1/2 other than the late-time data. They suggested using the squareroot-of-time plot to identify the linear flow. Sun et al. [28] investigated the influence of time-

Discussion
A mathematical model was developed in this study to capture the pressure-dependent decline of fracture conductivity from production data. As stated previously, the newly proposed model is suitable for wells with production in the boundary-dominated flow period. It has been recognized that the production rate declines with material balance time with slopes of −1/2 and −1 during the linear flow and boundary-dominated flow periods, respectively, in the log-log plot when there is no skin effect [33]. Many multi-fractured horizontal wells producing from shale reservoirs have been reported to experience a long-term linear flow (transient flow) period. The presence of decline of fracture conductivity and formation damage is usually interpreted as skin factor, which may change the shape of the diagnostic plot [22]. Nobakht and Mattar [34] stated that the linear flow is no longer a straight line with a slope of −1/2 other than the late-time data. They suggested using the square-root-of-time plot to identify the linear flow. Sun et al. [28] investigated the influence of time-dependent fracture conductivity on the production decline curves. Their study indicated that both formation damage and time-dependent fracture conductivity could not delay the time that the reservoir entered the boundary-dominated flow. Their influences on production decline are mainly seen at the early stage. Therefore, the log-log diagnostic method proposed by Zhou et al. [30] was used to identify the boundary-dominated flow regime rather than the linear flow regime in the present study. The log-log diagnostic approach is an effective method to avoid misinterpretation and to estimate the switching point.
The switching time from the transient flow to the boundary-dominated flow can also be estimated using the concept of radius-of-investigation [35], where t elf is the switching time in hours, φ is the porosity, c t is the total compressibility in psia −1 . From Equation (7) we can see that the switching time is highly related to the fracture spacing and the matrix permeability. A decrease in fracture spacing will reduce the switching time, provided that the matrix permeability is held constant. As the matrix permeability decreases, the expected switching time increases. Therefore, the long-term transient flow period may be attributed to either the ultralow permeability of the TMS or large fracture spacing.
Besides, it is assumed that the matrix permeability is constant. Besov et al. [2] hold the opinion that the TMS is different from other shale reservoirs because its pores cannot store hydrocarbons. Their experiments showed that the porosity of TMS formation is contained within its inorganic pores. They also indicated that hydrocarbons are more likely to be stored and produced from microfractures of TMS. It is well known that the permeability of microfractures exhibits stress sensitivity, which may delay the switching time from the transient flow to the boundary-dominated flow [28]. To simulate the fluid transport capacity of fractures during the whole life of hydrocarbon production where formation pressure declines and thus effective stress increases, the variation of microfracture permeability values need to be predicted. Our future work will focus on the conductivity of natural fracture and its influence on the production decline behavior of TMS.
Another limitation in our current study is that the proposed model can only estimate the pressure-dependent decline of fracture conductivity in the boundary-dominated flow period. Equation (A3) was developed for fractures filled with different materials such as clay (relatively mechanically unstable) and silica (mechanically stable). In the present study, we modified this equation to predict the conductivity of propped fractures. It is understood that hydraulic fracture conductivity is highly related to the closure stress, formation hardness, and proppant properties, etc. Once pumping is stopped, the fracture conductivity diminishes owing to proppant pack compaction and embedment of proppant particles during the production stage [22]. This means the compressibility of the proppant pack defined in equation (4) may change with production. In this case, the following equation can be used to describe the variation of the compressibility of the proppant pack [36], where α is the decline rate of the propped fracture compressibility with increased effective stress, c p0 is the initial compressibility of the proppant pack, σ e is the effective stress, σ e0 is the initial effective stress. As we stated in Appendix A, the propped fracture compressibility is assumed to be constant to simplify the solution. Results observed from the current study will be checked and further improved.

Conclusions
An analytical model was developed for capturing the pressure-dependent decline of fracture conductivity in the TMS fields from production data. The following conclusions are drawn: (1) The log-log diagnostic approach was used to identify the boundary-dominated flow regime of seven TMS wells. The results indicated that these seven wells exhibited a long-term transient flow period. (2) The R-square and AAREP were utilized in this study to estimate the performance of the proposed analytical well productivity model considering time-dependent fracture conductivity. The results indicated that the proposed model allows approximation of the real production data of the project wells in TMS formation with good accuracy. (3) The proposed fracture conductivity model was verified against simulated production data for a wide range of fracture conductivities/compressibilities. The results indicated that the predicted fracture compressibilities agree well with the prescribed values, with relative differences of fracture compressibility that were less than 10%. (4) The revealed average fracture conductivity decreases over time in the range of 100%-48% of the initial fracture conductivity. The pressure-dependent decline of fracture conductivity can be approximated using a logarithmic function given by equation (6)  initial fracture permeability, md k m matrix permeability, md N number of production data points n f number of hydraulic fractures N i initial oil in place in the drainage area of the well, stb N p cumulative oil production from the well, stb N p1 cumulative oil production at point 1, stb N p2 cumulative oil production at point 2, stb O i observed production rate, stb/d p i initial reservoir pressure, psia p average formation pressure, psia p w wellbore pressure, psia q o oil production rate, stb/d q o1 oil production rate at point 1, stb/d q o2 oil production rate at point 2, stb/d S f fracture spacing, ft t elf time at the end of transient flow, hour w fracture width, inch x f fracture half-length, ft α b Biot coefficient µ o oil viscosity, cp σ e effective stress, psi σ e0 initial effective stress, psi φ matrix porosity α decline rate of the propped fracture compressibility, psia −1 ν Poison's ratio where p is the average formation pressure, k f is hydraulic fracture permeability, and w is the average fracture width. It is generally believed that fracture conductivity plays an essential role in well performance. In the present study, we take the pressure-dependent decline of fracture into account. Chen et al. [36] developed a model to predict the conductivity of fractures filled with some porous materials such as carbonate and silica. It is assumed that the fracture conductivity of the propped fracture (k f ·w) can also be described by the following equation [36], To simplify the model, it is assumed that the compressibility of the proppant pack is constant [37]. Substituting Equation (A3) into Equation (A2) gives, Under the condition that reservoir pressure is above the bubble point pressure, based on the material balance equation, N p N i = e c t (p i −p) − 1 (A5) The average reservoir pressure is expressed as, To determine the values of c p in Equation (A8) and A in Equation (A9) using production data, select two points in the trend where, Solving Equations (A11) and (A13), one can get the values of c p and A. Then the variation of fracture conductivity during production can be determined by substituting Equation (A6) into Equation (A3), If desired, Equation (A9) can be used for predicting future production. The procedure is outlined as follows,