Research on calibrating rock mechanical parameters with a statistical method

Research on the modeling of rock mechanics parameters is of great significance to the exploration of oil and gas. The use of logging data with the Kriging interpolation to study rock mechanics parameters has been proven to be effective in reservoir prediction and other oilfield applications and can provide additional data. However, there will sometimes be a great deviation due to the limited samples and the strong heterogeneity of a layer. To solve this problem, a new approach was proposed to calibrate rock mechanical models through the statistical analysis of logging data. A module was developed to calibrate rock mechanics parameters automatically, which was then applied to the Wangyao area of the Ansai oilfield. This method significantly improved the accuracy of rock mechanics modeling.


Introduction
In the exploration and development of lithological oil-gas reservoirs, rock mechanics parameters, such as Young's modulus, Poisson's ratio, the effective stress coefficient, and the angle of internal friction, are crucial to geo-stress simulation, wellbore stability analysis, sweet spot prediction, and simulation techniques.In recent years, as the use of unconventional resources have developed rapidly, high-accuracy rock mechanics modeling has drawn a great degree of attention and been studied by experts worldwide, and many research results have been put into field application.
Considering the unbiased optimal advantages, the Kriging method is one of the most widely used methods in rock mechanics modeling.However, due to some limitations of Kriging in practical applications, many other methods have been developed to improve the applicability of Kriging, such as simple Kriging, ordinary Kriging [1], co-Kriging [2], universal Kriging [3], and indicator Kriging [4].In recent years, Roustant optimized the method to solve a covariance function by proposing a nonnegative solution of linear equations to eliminate part of the subjective impacts [5].Erum [6] proves that Bayesian universal Kriging fits better than the universal Kriging in predicting sodium concentrations.Hu improved the results of Kriging, which can be impacted by scaling, with the Bayesian-based collocated co-Kriging method [7].The above methods can improve the accuracy mathematically but fail to consider the spatial distributions of rock mechanics parameters and are unable to optimize the search radius and range.Therefore, these methods can cause a large error in the simulation results.To solve this problem, this study analyzed the behavior around the wellbore, constrained the Kriging results by setting up the law of geological parameters' distribution characteristics, and improved the mechanics modeling of rock with a calibration method.The module was also developed to constrain the rock mechanics parameters.This proposed method was applied successfully to the Wangyao area of the Ansai oilfield, and the results showed that the accuracy of rock mechanics modeling improved significantly.

Materials and characterization
The research reservoir is located in the mid-east region of the Ordos basin, the stratigraphic occurrence is gentle [8], the local structure is stable, and there are no major fault activities.The dipping magnitudes of the strata are approximately 0.5 degree; the average gradient is 8-10 m/ km.Differential compaction effects form a low angle nose-like uplifted structure.The delta sand-mud interaction affected the accumulation of hydrocarbons in this area.The reservoir has low porosity and permeability but is highly fractured [9][10][11], which results in high heterogeneity and complexity in the distribution of rock mechanics properties [12].Therefore, it is critical to create a high-resolution rock mechanics model based on the abundant well log data.
In this study, rock mechanics parameters at the wellbore were calculated using conventional logging data from more than 500 wells in the Wangyao area, together with some cross-dipole acoustic logging data (X-MAC).The file format of logging data is ".Las"; these data include the version information, using "~VERSION" as the identifier.Well information is separated by the "~WELL" logo.The curve information has a "~CURVE" flag; ASCII data are identified by "~A", which includes the actual data on logging.

Structure of rock mechanics model
Traditional rock mechanics modeling fails to consider the whole spatial distribution trend of rock mechanics parameters by simply using the neighborhood properties and local locations.Therefore, this study restricted rock mechanics modeling by analyzing Young's modulus and Poisson's ratio at wellbores from the macroscopic view to achieve practical results.
This study consists of four parts (Fig 1).First, data were collected and rock mechanics parameters at wellbores were calculated, such as Young's modulus and Poisson's ratio, using correlation analysis and the regression approach to address cross-dipole acoustic logging and density logging.Then, the initial rock mechanics model was built using ordinary Kriging, including searching the neighborhood, solving the covariance function, and meshing.The original data for Kriging come from the rock mechanics parameters at wellbores.Statistical analyses of Young's modulus and Poisson's ratio were carried out.Then, a rock mechanics parameter restriction model was built to better calibrate the model.Finally, a module was developed with the programming language C++ to calibrate the Young's modulus and Poisson's ratio and improve the accuracy of the rock mechanical model.

Calculation of rock mechanics parameters at the wellbore
The current methods to collect and calculate rock mechanics parameters are mature and mainly include the measurement method [13] and logging operation method [14].During the exploration and development of oil fields, conventional logging plays an important role in extracting rock mechanics parameters due to its rich data, low cost, and higher accuracy than seismic data.It also contains important information, such as the reservoir structure, lithology, and physical properties.
In this study, rock mechanics parameters at the wellbore were calculated using X-MAC logging [15] together with density logging, and Young's modulus and Poisson's ratio were calculated with the following equations: where V refers to Poisson's ratio, dimensionless; E refers to Young's modulus, MPa; ρ is the density, g/cm3; and Δt s and Δt p denote the shear wave slowness and compressional wave slowness, respectively.Using X-MAC logging data to calculate Young's modulus and Poisson's ratio brings feasibility to revising the final results, which were calculated with acoustic logging, gamma ray logging, and spontaneous logging, and eventually the Young's modulus and Poisson's ratio were established through a regression (Fig 2).
According to the correlation analysis, there is a strong correlation between Young's modulus/Poisson's ratio and acoustic logging/gamma ray logging, which can bring good regression results.As shown in Fig 3, the value increases as the color scale darkens.Based on acoustic logging and gamma ray logging data, Young's modulus and Poisson's ratio were calculated with the following regression equations: where AC refers to the wave slowness, μs/m, and GR refers to the gamma ray data, API.The correlation coefficients from the above two equations are 0.996 and 0.998, respectively, which confirm the strong correlation between Young's modulus/Poisson's ratio and acoustic logging/gamma ray logging.After calculating the dynamic Young's modulus and Poisson's ratio of wells in the Wangyao area using the above equations, they were converted to static parameters using equations provided by Du [16].Then, they can be used to change the stress field in a low-permeability reservoir and secondary exploitation.

Constructing the initial rock mechanics model
Kriging interpolation is the core of geostatistics, and it is widely used in the field of geoscience.Kriging was first introduced by the South African mining engineer D.G. Krige in 1951 to resolve the problems of deposit reserve computation and error estimation [17].Considerable research on Kriging has been published globally since then.The basic idea of Kriging is to predict the value by assigning weights to known sample points according to their spatial locations and correlations relative to the target point and calculating the moving weighted average of known values.The Kriging method takes advantage of property correlations and spatial locations between neighboring sample points, and it provides high interpolation estimation accuracy; therefore, this method has been used in the field of rock mechanics parameters modeling [18][19][20][21].
The theoretical model of the Kriging interpolation algorithm is expressed as follows: where λ i (x 0 ) are the weight coefficients, x is the location of the unknown point, x k is the location of a known point, n is the number of measurement points that have a relationship with x, G Ã (x 0 ) is the estimate of the unknown point, and G(x k ) is the measurement value of the known point k.The equations for calculating weights λ i (x 0 ) are as follows: where γ(x i −x j ) is the variogram between known points i and j and μ is a Lagrange multiplier.In this paper, the steps for using the Kriging method to construct rock mechanics parameters model were proposed and include inputting wellbore data, generating a 3D geo-model mesh, fitting the variogram model, searching the neighborhood data, and performing interpolation.The method for creating an initial rock mechanics model using Kriging is as follows: 1. Acquire the wellbore data.Wellbore rock mechanics parameters, including Young's modulus and Poisson's ratio, are calculated from conventional logging data and prepared in standard data formats for Kriging interpolation.In addition, due to the existence of abnormalities in log data, outliers will be processed or eliminated.
2. Construct the geological mesh.The stratigraphic model [22] is represented on a triangular surface mesh [23][24][25][26][27][28], but the Kriging interpolation of rock mechanics parameters is performed on volume mesh elements, so it is necessary to convert the triangular mesh geomodel into a volume mesh geomodel [29][30][31][32].In this paper, the hexahedral element is selected as the basic unit, and each stratigraphic layer of the geomodel is constructed by mesh voxelization.
3. Fit the variogram model.Given the same source data and geomodel mesh, variograms influence the accuracy of Kriging interpolation.For layers with high heterogeneity, it is necessary to compute variograms in multiple directions to achieve higher interpolation accuracy.In this paper, the empirical variograms were computed with sample data.We select the spherical model as the variogram model to simplify the procedure and use least square regression to find the parameters in the variogram model.4. Search the neighborhood.To calculate weights in Eq 6, it is necessary to search for grid cells in the neighborhood of the unknown point to obtain the Young's modulus and Poisson's ratio of the sample cell.However, when too many grid cells are found, the search efficiency is low; therefore, an algorithm was developed to find the closest grid cells to the point of interest to improve the modeling efficiency.

Building the rock mechanics parameter constraints model
There are more than 500 conventional wells in the whole Wangyao area.With this high well density, Young's modulus and Poisson's ratio at the wellbore can reflect the spatial distribution of rock mechanics properties, which provides a solid source to build a three-dimensional rock mechanical model.The Young's modulus and Poisson's ratio were analyzed mathematically and statistically, and they were used to calibrate the model.In this study, the statistical factors have been considered, including the maximum value, minimum value, average, variance, and probability distribution.The variance, average and histogram of Young's modulus and Poisson's ratio at the wellbore and from the initial model were compared to better construct the theoretical model.Assume that the maximum, minimum, average, variance, frequency of the wellbore Poisson's ratio and initial Poisson's ratio are (max 1 , min 1 , m 1 , σ 1 , t 1 ) and (max 2 , min 2 , m 2 , σ 2 , t 2 ), respectively.The differences Δm = m 1 − m 2 and Δσ = σ 1 − σ 2 were calculated.Most importantly, the following linear approximation equation was utilized to optimize the distribution of initial rock mechanics parameters: where Pois 1 (x) is the optimal Poisson's ratio at location x; Pois 2 (x) is the initial Poisson's ratio at location x; and k 1 (k 1 >0) and k 2 (k 2 >0) are coefficients.

Development of calibration module
To calibrate the rock mechanics model efficiently, a module was developed to realize the automatic calculation of properties.This module consists of reading data from the model, defining the constraint model, calibrating the model and managing the data.Due to the varying properties in a geological grid, such as lithological and geophysical properties, an efficient data reading module needs to be established that allows us to select different geophysical properties.
After loading the data from the last step, the calibration model was applied to modify all the properties.Different areas and different properties correspond to different calibration models.Therefore, the module was designed to customize the calibrated model so that users can input the coefficients of the model.This module also offers the feature of property extraction through the interface after calculation.Some simple calculations and value inputs are provided by the "property calculator" as shown in Fig 6, which can easily calibrate the model with a single property or multiple properties.
One goal of model calibration is to convert the model and properties into a format that can be read by computers.When there are multiple geological mesh models, properties cannot be read at once.Therefore, the model was calibrated block by block.As shown in  After the initial rock mechanical model has been calibrated, it is necessary to save a new model and export property data.In this study, the mesh coordinates and rock mechanics properties were saved in text format, and the Young's modulus and Poisson's ratio were saved at different physical locations.

Analysis
Based on the frequencies t 1 and t 2 , the corresponding p 1 and p 2 can be calculated, which are the probability of the wellbore Poisson's ratio and initial Poisson's ratio at the same Poisson's ratio value.Therefore, under each Poisson's ratio, there will be two sets of probability.Part of the values are shown in Table 1.
Two sets of probability are plotted in Fig 9 .The left picture is the probability distribution of the wellbore Poisson's ratio, and the right picture is the probability distribution of the initial Poisson's ratio.The vertical coordinates represent the probability, whereas the horizontal coordinates represent the Poisson's ratio value.
The linear relationship between the two sets of probabilities was established with a regression.As shown in Fig 10, the correlation coefficient is 0.8872.Assuming the slope in the linear regression equation is K, k 1 can be solved by k 1 ¼ K= ffiffiffiffiffiffi Ds p , and then, k 2 can be solved by The accuracy of the initial Poisson's ratio model largely increases after the calibration, and the initial Young's modulus model can also be calibrated in the same way.The workflow of the model calibration can be summarized as follows:  1. Collect all the data related to rock mechanics properties from the wellbore and the initial model.
3. If the differences (m 1 -m 2 ) and (σ 1 -σ 2 ) are within the acceptable error, stop the workflow, otherwise, continue.The error set in this study is 10%.
4. Apply correlation analysis between the probability of the wellbore Poisson's ratio and the initial Poisson's ratio given the same Poisson's ratio value 5. Solve coefficients k 1 and k 2 using the regression method, and then, build a constraint model of rock mechanics parameters.The Poisson's ratio of the wellbore, the initial and calibrated models are summarized in Table 2.The maximum value, average, and variance at maximum probability of Poisson's ratio after calibration are closer to those of the wellbore model.Therefore, the proposed method can calibrate the initial Poisson's ratio effectively and improves the model accuracy.Fig 12 shows the models before and after the calibration, demonstrating the uneven distribution of the Poisson's ratio obviously in this area, because more accurate geological information can be captured after calibration.
A data interface was developed to import the calibrated parameters and geological model to finite element software ANSYS, and the geostress was calculated successfully on a 3D finite model.Fig 13 shows the minimum and maximum geostress distribution.

Conclusions
In this study, using conventional logging data and X-MAC logging data, rock mechanics parameters were calculated using a regression, a geological mesh model was built, and the  initial three-dimensional rock mechanics model was constructed using Kriging.Based on this work, a new approach for calibrating three-dimensional models of rock mechanics parameters was proposed.Through analysis of the rock mechanics parameters of the wellbore and initial model, the optimal coefficients were calculated, which were important in building the constraint model.Then, the rock mechanics parameters could be calibrated using the constraints model; in this way, the accuracy of the initial mechanics model using the Kriging interpolation method was significantly improved.

Fig 3 .
Fig 3. Correlation between conventional log curves and rock mechanics parameters in the Ansai oilfield.https://doi.org/10.1371/journal.pone.0176215.g003 Fig 4 is a graphic interface of the variogram fitting function developed with the language C. Fig 4 shows that the sphericalvariogram is selected because it can be easily converted to a linear variogram to simplify calculation and because it is convenient for finding parameters such as nugget and sill[33- 34].The lag distance, lag tolerance and number of lags are required inputs for variogram fitting.Meanwhile, parameter isotropy or anisotropy is selected based on the actual scenario.The initial values of nugget, sill and range are calculated from the variogram fitting.These values are usually optimized according to the interpolation results.The Kriging calculation can be done after all required settings are provided.

Fig 4 .
Fig 4. Fitting of variogram function.https://doi.org/10.1371/journal.pone.0176215.g004 Fig 7, the whole model was first divided into slices, and then each slice was divided into columns to improve the efficiency of data reading.Fig 8 demonstrates the procedure of the dividing models.

Fig 11 (
Fig 11(b) shows the histogram of the Poisson's ratio probability distributions after calibration.The probability distribution after calibration moves to the right to a certain extent compared with the probability distribution of the initial model in Fig 11(a).The Poisson's ratio corresponding to the maximum probability in Fig 11(b) is 0.288, which coincides with the Poisson's ratio in Fig 11(c), which shows the histogram of the Poisson's ratio probability distribution at the wellbore.In addition, the distributions of Poisson's ratios on the left side of the maximum probability in Fig 11(b) and 11(c) seem to be similar.In total, the Poisson's ratio probability distribution after calibration matches better with the distribution of the wellbore parameters than the initial model.The Poisson's ratio of the wellbore, the initial and calibrated models are summarized in Table2.The maximum value, average, and variance at maximum probability of Poisson's ratio after calibration are closer to those of the wellbore model.Therefore, the proposed method can calibrate the initial Poisson's ratio effectively and improves the model accuracy.Fig 12 shows

Table 2 . Comparison of Poisson's ratio in different models. Poisson's ratio of different situations Poisson's ratio at maximum probability Minimum value Maximum value Average Variance
https://doi.org/10.1371/journal.pone.0176215.t002