A Quantitative Prediction Method for Fracture 1 Density Based on the Equivalent Medium Theory 2

Fracture density, a critical parameter of unconventional reservoirs, can be used to evaluate 11 potential of unconventional reservoirs and location of production wells. Many technologies, such 12 as amplitude variation with offset and azimuth (AVOA) technology, vertical seismic profiling (VSP) 13 technology, and multicomponent seismic technology, are generally used to predict fracture of 14 reservoirs. they can qualitatively predict fracture by analyzing seismic attributes, including seismic 15 wave amplitudes, seismic wave velocities, which are sensitive to fracture. However, it is important 16 to quantitatively describe fracture of reservoirs. In this study, based on a double-layer model, the 17 relationships between fracture density and the double-layer model’s physical parameters, such as 18 velocity of fast shear-wave, velocity of slow shear-wave, and density, were established, and then a 19 powerful quantitative prediction method for fracture density was proposed dramatically. 20 Afterwards, the Hudson model for crack was used to test the applicability of the method. The result 21 shown that the quantitative prediction method for fracture density can be applied suitable to the 22 Hudson model for crack. Finally, the result of validation models indicated that the method can 23 predict fracture density effective, in which absolute relative deviation (ARD) were less than 5% and 24 root-mean-square error (RMSE) was 4.88×10-3. 25


Introduction
Unconventional resources have received widespread interest in recent years because of their emergence as a source of clean energy, and they are now being developed and produced in North American shales, Southern North China shales, etc. [1][2][3].As self-sourced reservoirs, the fracture space controls not only the potential of unconventional reservoirs but also location of production wells [4].Thus, the accurate prediction of fracture density plays an important role in the exploration and development of unconventional reservoirs.In particular, knowing the fracture density variation between zones within a reservoir can help in the accurate determination of well locations [5].
Many achievements have been made with the continuous improvements by researchers for predicting the fracture in sedimentary strata [6][7][8][9].Sedimentary strata, such as unconventional reservoirs, are usually described as anisotropic mediums on the scale of seismic waves, and the anisotropy of a medium is produced by skeleton orientated arrangement and fracture development [10][11].Ruger (1997Ruger ( , 2002) ) researched the relationship among P-wave reflection coefficients, offset, and azimuth named AVOA technology, which was found to be suitable for analyzing the anisotropy of a medium [12][13].While, some researchers used the AVOA technology to predict the anisotropy of fractured reservoirs for analyzing the fracture development, and the results in these studies were presented by anisotropic parameters to reflect the fracture development [9,[13][14][15].Thus, the AVOA technology was a non-direct way to predict the fracture development of reservoirs clearly.
Using the vertical seismic profiling (VSP) technology to predict the fracture development of reservoirs is realized by the splitting property of the S-wave, the S-wave split into fast S-wave and slow S-wave when passing through fracture with an angle [16][17], which is obviously different from AVOA technology.The difference between the velocities of these waves can be estimated by measuring the increase of the time delay between them with the depth [18].Meadows et al. [19], Pevzner et al. [20], and Shevchenko et al. [21] used VSP technology to predict the fracture of reservoirs.The prediction results shown that the VSP technology can estimate the S-wave anisotropy of reservoirs by the difference between the fast S-wave velocity and the slow S-wave velocity, which was useful for representing the fracture development of reservoirs.However, this technology requires a borehole that is available for the exploration area.Thus, it only has the ability to reflect the fracture development near the borehole, actually.Ata et al. [22], Gaiser [23] and Vetri et al. [24] used multicomponent seismic technology to predict the fracture development of reservoirs by Swave splitting property.Compare with the VSP technology, multicomponent seismic data did not have the advantage of recording data near the reservoirs, which can bypass complications introduced by the overburden.Fortunately, with the improvement of multicomponent seismic data acquisition and processing technology [25][26][27][28][29][30][31], the application of multicomponent seismic technology has gradually increased, providing favorable technical support for identifying the fracture of reservoirs [32][33].However, Using the splitting property of the S-wave to predict the fracture of reservoirs is also a non-direct way yet.
During the past few decades, the equivalent medium theory, a method for digitizing rocks, has attracted the attention of many researchers [34][35][36].Backus [37] demonstrated that a multilayer transversely isotropic medium was equivalent anisotropic if the wavelength of seismic wave is much larger than the monolayer of the multilayer transversely isotropic medium, and then the equivalent elastic tensor of this medium was given.Eshelby [38] studied the situation of a single elliptical inclusion in isotropic media, and then Cheng [39] gave the equivalent elastic tensor of the situation.
Hudson [40] researched the condition that fracture is represented as a gap or inclusion with a thin coin-shaped ellipsoid, and provided the equivalent elastic tensor.Schoenberg [41] ignored the shape and microstructure of fracture and considered the fracture was an infinitely thin and very soft stratum satisfying the linear sliding boundary condition.While, a linear sliding model was proposed.
Based on the equivalent medium theory, this study established the relationships among fracture density, multicomponent seismic response, density.As a result, a powerful quantitative prediction method for fracture density was proposed.Furthermore, the applicability of the quantitative prediction method for fracture density was verified by the Hudson model for crack.In this study, Section 2 addresses details to construct a quantitative prediction method for fracture density based on an equivalent medium model.Section 3 verifies the quantitative prediction method for fracture density with the Hudson model for crack.The data to test the proposed method and discussion in this study are presented in Section 4. Finally, conclusions are presented in Section 5.

Method Assumptions and Establishment
In this section, firstly, we make some idealized assumptions in order to achieve a powerful method for predicting fracture density.Secondly, the relationships between the fracture density and the physical parameters of target layer, such as fast S-wave velocity, slow S-wave velocity, root mean square (RMS) velocity, and density, are developed under the idealized assumptions.Finally, calculation Eqs. for the fracture density are given in detail, and then a quantitative prediction method for fracture density is established.

Assumptions
As a result of their genesis, sedimentary strata, such as unconventional reservoirs, can be exhibited in a thin interbedded medium.In this work, we assume that the thin interbedded medium, shown in Fig. 1, is composed of two components: rock skeleton layer (matrix) represented by black in Fig. 1 and rock fracture layer (including filling, water, oil, and gas) denoted by gray in Fig. 1.Moreover, following assumptions and simplifications are made to generate a practical method for predicting the fracture density which is defined by the volumetric proportion between the fracture layer and the total layer.(2) There are no sources of intrinsic energy dissipation, such as friction or viscosity, between the rock skeleton layer and the rock fracture layer.

Preprints
(3) The thickness of the rock skeleton layer and the rock fracture layer must be much smaller than a seismic wavelength, which the seismic wavelength must be at least ten times a layer thickness.
According to the aforementioned assumptions, we can divide the thin interbedded medium into two parts: a skeleton layer and a fracture layer, and then the thin interbedded medium can be equivalent to a double-layer medium.Fig. 2 shows the double-layer medium, in which the black part represents the skeleton layer and the gray part indicates the fracture layer.

The relationships between S-wave velocities and fracture density
Notable, Backus (1962) demonstrated that if the wavelength of seismic wave is much larger than the monolayer of the thin interbedded medium, it will be equivalent anisotropic [37].Thus, the stiffness tensor ( C ) of the double-layer medium can be represented facile by Backus Average, which depicted with the P-wave velocity ( P V ), the S-wave velocity ( S V ), the density (  ), and the volume of each layer.Tab. 1 shows the physical parameters of the double-layer medium, and then the stiffness tensor ( C ) of the double-layer medium can be described through Backus average [37] and Levin formula [47], which has the following form: The brackets  indicates averages of the enclosed properties weighted by their volumetric proportions, which is often called the Backus average.As a result, we can represent the fast S-wave velocity ( Fast V ) and the slow S-wave velocity ( Slow V ) of the double-layer medium as follows: is the fracture density with the range [0, 1), dimensionless.

The relationships between root mean square (RMS) velocities and fracture density
On the basis of above-mentioned, we can develop the relationships between root mean square (RMS) velocities and the fracture density.In this study, the P-wave RMS velocity ( P V ) in the double- layer medium can be expressed as: represents the vertical propagation time of P-wave in the skeleton layer, s; denotes the vertical propagation time of P-wave in the fracture layer, s.Moreover, the physical parameters of the skeleton layer and the fracture layer, shown in Table 1, are substituted into Eq.( 4), we can obtain: Through further simplification, the Eq. ( 5) can be rewritten: In the same way, the S-wave RMS velocity ( S V ) can be obtained as follows:

The relationship between P-wave velocity and density
Gardner [48] and Castagna [49] studied the relationship between the P-wave velocity and the density in different lithology, and then determined their relationships.Gardner [48] suggested an empirical relation between the P-wave velocity ( G V ) and the density ( G  ) that represents an average over many rock types: where G V is in ft/s.

The quantitative prediction method of fracture density
Through aforementioned Eqs. on the relationships between the fracture density and the physical parameters of the double-layer medium, such as the fast S-wave velocity, the slow S-wave velocity, the root mean square (RMS) velocity, and the density, the key issue required to be solved is how to establish a practical method for predicting fracture density.
After the analysis on the correlations of the Eqs.above, we consider to build the relationship between the P-wave velocity and density of a target layer firstly, which is according to the studies by Gardner [48] and Castagna [49].Actually speaking, a certain amount of core analysis and logging data will be finished before development of an unconventional reservoir, and as a result the nonlinear relationships between the P-wave velocity and the density of the target layer can be obtained: where coefficient a reflects the relationship between the P-wave velocity and the density of the skeleton layer, dimensionless; coefficient b indicates the relationship between the P-wave velocity and the density of the fracture layer, dimensionless.Meanwhile, the average density ( all  ) of the target layer can also be obtained: Furthermore, taking Eq. ( 10) and Eq. ( 11) into Eq.( 6), and its expression Eq. was as follows: Notable, the RMS velocities about the P-wave and S-wave, and the velocities about the fast Swave and the slow S-wave of the prospecting stratum are acquired from a conventional data processing of multicomponent seismic exploration.Moreover, the average density can be obtained from the core analysis and logging data of survey area.Finally, the fracture density can be achieved through solving the Eq. ( 2), Eq. ( 3), Eq. ( 7), Eq. ( 12), and Eq. ( 13), and the flowchart of the calculation is illustrated in Fig. 3.According to the aforementioned analysis, we established the quantitative prediction method of fracture density for the first time in this way.
The key for applying and verifying the new method to the Hudson model is how to represent the Eq.
(   Skeleton part According to above analysis, we consider that the 1  , and  of the double-layer model are equal to , and H  of the Hudson model, respectively.By adjusting the 1

Unique solution set determination
In order to obtain the unique solution set of the The Eqs. ( 14) can be expanded into the following expression: For further simplification, we have defined some intermediate variables shown in Eqs.(17), and then Eqs. ( 15) can be rewritten: where There are two sets of solutions about the variables including the 1 S , 2 S , 1 P , and 2 P which has a certain relation with the V , and 2 S V , respectively, through solving the Eqs.(16), and its expressions are as follows: where apostrophe denotes the first set of solutions; asterisk represents the other set of solutions.Substituting Eqs.(17) ) According to the Appendix A, Eqs. ( 20) is not the reasonable solution for the  V under the requirement that the S-wave velocity of the skeleton part (  V .Moreover, it is obviously that there will be have 16 sets of solutions for the  V , due to the four parameters at the left of the Eqs.( 21) are quadratic.Finally, the Eqs.( 22) is the only solution about the Eqs.( 21) with the positive number, owing to the other solutions are eliminated under the actual situation that the seismic wave V )must be real numbers larger than zero, which is as follows: Based on the above analysis, there is a set of the P-wave and the S-wave velocities ( 1 ) that make the stiffness tensor of the double-layer model and the stiffness tensor of the Hudson model equal.Meanwhile, the density of the skeleton part, the density of the fracture part, and the fracture density of the two models are equal.Thus, the quantitative prediction method of fracture density can be applied to the Hudson model.

Validation model
In order to examine the effective of the quantitative prediction method of fracture density, ten groups of test models are established to validate the prediction method of fracture density in Tab. 3.
Due to the assumptions that the rock skeleton layer and rock fracture layer of a medium are isotropic, homogeneous, and linear elastic, etc.Ten groups of commercial numerical models are performed to validate the application effect of the proposed method in this paper.Each group of test model contains two different core samples.Moreover, the core sample with high velocity and density is defined as a skeleton layer, and the low is defined as a fracture layer, which is conform to the assumptions in Sec.2.1.As a result, each group of model can be used as a double-layer model.Tab.

283
shows the fast S-wave velocity ( Fast V ), the slow S-wave velocity ( Slow V ), the P-wave RMS velocity ( 284 P V ), the S-wave RMS velocity ( S V ), the average density ( all  ), the coefficient a, and the coefficient b, 285 which are calculated with the Eq. ( 2), Eq. (3), Eq. ( 6), Eq. ( 7), Eq. ( 12), Eq. ( 10), and Eq.(11), respectively.coefficient b can be achieved from the core analysis and logging data of survey area.Finally, the fracture density (  ) is predicted by the method of fracture density through solving the Eq. ( 2), Eq.
Table 4 The physical parameters of models

Results
According to the quantitative prediction method of fracture density in Sec.2.5 and the aforementioned physical parameters of models in Tab. 4, we can predict the fracture density (  ).
Moreover, we used the quasi-Newton method to solve the Eq. ( 2), Eq. ( 3), Eq. ( 7), Eq. ( 12), and Eq. ( 13) for predicting the fracture density in this study.Fig. 5 shows the predicted values of the fracture density ( '  ).It can be seen that the predicted fracture density is 0.24-0.25 for the ten models, which indicates the effective of the predicted values for the fracture density.The minimum value located in the seventh model, with a value of about 0.241, and the maximum value, 0.249 is located in the fourth model.In addition to the predicted values for the fracture density, the root-mean-square error (RMSE) and absolute relative deviation (ARD) were used as indexes to evaluate the prediction effects.The smaller the root-mean-square (RMSE) and ARD were, the more favorable the prediction effect was.
The root-mean-square error (RMSE) was used to weight the deviation between the prediction result and the actual value, and its calculation Eq. was Eq. ( 23).The ARD was usually used to evaluate the accuracy of method through comparing the prediction value and the actual value.Its calculation Eq. is Eq. ( 24), as follows: In Eq. ( 23) and Eq. ( 24),  Fig. 6 shows the results of RMSE and ARD.It was shown that ARD were less than 5% and RMSE was 4.88×10 -3 .It indicated that the fracture density predicted by the method developed in this study were in agreement with the actual data.The results demonstrated that the quantitative prediction method proposed in this study was capable for predicting fracture density data.
In this study, for establishing a practical method to predict the fracture density quantitatively, we assume that a thin interbedded medium is composed of two components: rock skeleton layer (matrix) and rock fracture layer (including filling, water, oil, and gas), and the rock skeleton layer and rock fracture layer of the thin interbedded medium are isotropic, homogeneous, and linear elastic, etc., which is an idealized hypothesis.As a result, the thin interbedded medium can equivalent to the double-layer medium under these assumptions.Actually speaking, these assumptions will involve in deviation to this method, and the effect of these assumptions on the prediction method of fracture density will be the focus in the future research.However, it is significant that the new method can provide a scale in predicting the fracture density.Besides, the new method can be applied to the Hudson model, and then the scope of the prediction method of fracture density is extended, which is verified in Sec. 3.
Based on the above assumptions, the relationships between fracture density and the seismic response is established.In Sec.2.2-2.3,we use the fast S-wave velocity, the slow S-wave velocity, the P-wave RMS velocity, and the S-wave RMS velocity as constraints to predict the fracture density, which can be obtained from a conventional data processing of multicomponent seismic exploration.
Therefore, the processing quality of the multicomponent data will affect the application of the prediction method for fracture density.Besides, only concerns the wave velocity as a function of density in Sec.2.4 is imprecise, which is to simplify the prediction method.In fact, the velocity is related to several factors among which for instance the Young's modulus.Moreover, there will be an error in predicting the average density of a target layer by the core analysis and logging data.The influence of these factors on the prediction results will be reported in the later article.

Conclusions
The accurate prediction of fracture density plays an important role in the exploration and development of unconventional reservoirs which have received widespread interest in recent years.In this study, we considered that a thin interbedded medium is composed of the rock skeleton layer and rock fracture layer, besides the rock skeleton layer and rock fracture layer are isotropic, homogeneous, and linear elastic, etc.Under the above assumptions, the thin interbedded medium can be equivalent to a double-layer medium, and then the relationships between the fracture density and the physical parameters of a target layer, such as the fast S-wave velocity, the slow S-wave velocity, the root mean square (RMS) velocity, and the density, were established, which makes up the gap between the fracture density and the seismic response.While, a quantitative prediction method for fracture density was proposed.Afterwards, the Hudson model was used to test the applicability of this method.The result shown that the new method can be applied to the Hudson model.Moreover, the quantitative prediction method of fracture density was used to the ten groups of validation model, in which RMSE and ARD were as the index of evaluation.The RMSE=4.88×10 - 3 demonstrated that the fracture density predicted by the method developed in this study are in satisfactory agreement with the actual data.Meanwhile, the ARD of the fracture density were less than 5%, which indicated that this method can accurately predict fracture density.All in all, a powerful method for predicting fracture density is proposed in this study, and the method can be used to predict the fracture density of unconventional reservoirs, which provided a new idea for the study of reservoirs fracture.

Figure 1 .
Figure 1.Thin interbedded medium (1) The rock layer and rock fracture layer of the thin interbedded medium are assumed isotropic, homogeneous, and linear elastic.

Figure 2 .
Figure 2. Double-layer medium the skeleton layer and the total layer with the range (0, 1], dimensionless; 2 f   indicates the volumetric proportion between the fracture layer and the total layer with the range [0, 1), dimensionless;

Figure 3
Figure 3 Flowchart of the calculation ) of the double-layer model by the Hudson model parameters.The section will discuss how to solve Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 16 April 2018 doi:10.20944/preprints201804.0208.v1 the velocities about the double-layer model by the Hudson model parameters according to the Eq.(1).

Figure 4 .
Figure 4. Hudson model for crack

and 2 SV 1 S V , 2 PV , and 2 SV
to make the stiffness tensor of the double-layer model equals to the stiffness tensor of the Hudson model, and then the equivalent RMS velocities about the P-wave and S-wave of the Hudson model can be obtained, which is the key for applying and verifying the new method to the Hudson model.As a result, the above issue is converted to another form that is whether an unique solution set of the 1 P V , is existence to make the stiffness tensor of the double- layer model equal to the stiffness tensor of the Hudson model. and

V
) should be larger than the S-wave velocity of the fracture part (

V
) in the actual situation.Hence, the Eqs.(21) is the Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: and

279
in situ indicates in-situ measurement.

280 4 . 2
Validation model forward 281Tab.3 provides ten groups of model, which are regard as the double-layer model.Thus, the 282 physical parameters of the models can be calculated according to the above Eqs. in Sec. 2. Tab. 4


represents the prediction fracture density; i  denotes the actual fracture density; and n indicates the number of the validation models.

Figure 6 .
Figure 6.The RMSE and ARD of the predicted fracture density.
. (A1), and then the both sides of Eq. (A1) are divided by both side, we can obtain: further simplification, the Eq.(A3) can be rewritten:

2 SV
are the S- wave velocities about the skeleton part and the fracture part, respectively, according to the actual conditions.

Table 1
Physical parameters of the double-layer medium

Table 2
Physical parameters of the Hudson model for crack and 2

Table 3
The core sample parameters of models 278