Abstract
A novel optimised back analysis method is proposed in this paper. The in situ stress field of an underground gas storage (UGS) reservoir in a Turkey salt cavern is analysed by the basic theory of elastic mechanics. A finite element method is implemented to optimise and approximate the objective function by systematically adjusting boundary loads. Optimising calculation is performed based on a novel method to reduce the error between measurement and calculation as much as possible. Compared with common back analysis methods such as regression method, the method proposed can further improve the calculation precision. By constructing a large circular geometric model, the effect of stress concentration is eliminated and a minimum difference between computed and measured stress can be guaranteed in the rectangular objective region. The efficiency of the proposed method is investigated and confirmed by its capability on restoring in situ stress field, which agrees well with experimental results. The characteristics of stress distribution of chosen UGS wells are obtained based on the back analysis results and by applying the corresponding fracture criterion, the shaft walls are proven safe.
1. Introduction
In situ stress is an important factor that should be considered during the construction of underground gas storage (UGS). A direct way to get the in situ stress is to measure the principal stress by test methods, such as micro-fracturing test [1].
As the engineering technology develops, kinds of techniques have been selectively accepted for obtaining accurate in situ stress orientation and magnitude. Hydraulic fracturing, when poroelastic effects are considered, seems to be the best option to derive the horizontal in situ stress magnitude and direction [2]. If a vertical fracture is created in the process, then the minimum in situ stress must be in the horizontal direction. Although fracturing has been the most commonly used method in in situ stress measuring, it may not be available or even fail sometimes, for example, in inclined test holes where the fracture direction is not normal to the minimum horizontal stress. On the other hand, the fracturing method will also be ineffective when the measurement requires a higher injection pressure and the terrain, with UGS located, is bearing a very large stress. Some researchers implement methods such as leak-off tests and breakout orientation to the standard well testing and logging; however, these methods are not always reliable. Moreover, the stresses obtained from breakout analysis are confusing as breakout can be induced by different kinds of in situ stresses.
When in situ stress data from test wells is obtained, some techniques can be used to identify the magnitude and direction of the whole in situ stress field [3]. These techniques are organised as an overdetermined system of equations, and the in situ stress state can be solved by using several methods mentioned in some articles [4–6]. Not only can the magnitude of the two principal stresses be solved, but also the direction of the field can be determined. However, the methods assume that the shear stress is negligible as it is zero when the normal stress is a principal stress. Djurhuus and Aadnøy extended the study on in situ stress and presented a new theory [4]. They obtained experiment data from multiple fracturing and induced fractures from imaging logs to calculate the in situ stress state. The magnitude information is mainly calculated from fracturing data, while the orientation information can be obtained from image logs.
However, measurement is difficult to carry out because of high cost on an injection-production well. Therefore, several in situ stress values can only be obtained through experiments. The stress distribution of the test zone can be calculated quickly with the help of finite element optimisation and the results can be more accurate with some optimisation methods implemented. To carry out the inversion calculation, we need to construct the finite element model with proper boundary conditions and loads. Some minor modifications should also be incorporated into the model to eliminate the effect of stress concentration, which will be discussed below.
Boundary conditions and loading methods cannot be input directly because the tectonic movement is unknown and the geologic structure is complicated. These circumstances make in situ stress field analysis difficult to carry out. So boundary load inversion comes to be the crucial way to conduct research on the in situ stress field. From the viewpoint of geomechanics, the tectonic stress field is the in situ stress field forming the tectonic system and type, including its distribution area and inside stress distribution. During the analysis, rock density could be determined relatively accurately in a small variation range with different precision levels of gravity stress and tectonic stress. Instead of being inversed, based on terrain conditions, the finite element method could be used to calculate gravity stress to ensure the precision can satisfy the engineering requirement. However, tectonic stress field calculation is different [7]. During the construction of UGS, the distribution of the in situ stress field along the plane at the top of the salt cavern is needful while the in situ stress information of only a fraction of wells in a certain area is available. Then, inversion and extrapolation need to be performed to obtain the stress field of the whole area. Among the present numerical simulations and analysis methods of the in situ stress field, the most widely used approach is cut-and-try and adjustment of boundary loads [8–13]. Given the complexity of the geological structure and the limited reference of the in situ stress, boundary conditions and loading methods become the difficult problem in calculation and analysis part of a numerical simulation. The major issue is how to improve the precision of the inversion and reliability of the data. Simulation is an effective way to take terrain and geological features into consideration. And only by this can relevant theoretical or numerical analyses such as back analysis, back calculation, and simulation be performed to predict the tectonic stress and its mode based on the measured data [12]. The boundary loads identification and optimisation method presented by this paper adopts damped least squares algorithm to approximate the loads by setting parameters from stress testing references. The initial value for iteration is got by the unit load method to avoid the limits from the structure geometric features and boundary conditions. By controlling the error of observing point, the identification on uncertain and unknown boundary loads of the in situ stress field is achieved. A few ways are attempted and compared to improve the data reliability and the precision of in situ stress field inversion calculation in this paper.
2. Building the Finite Element Model
It is difficult to achieve an analytical solution of boundary loads for a formation due to the complexity of the geological structure and boundary conditions. The finite element method is employed here for mechanical analysis to make sure that the research method can have a general application.
Though salt is viscoelastic, a UGS reservoir stress field can be treated as an elastic problem to simplify the analysis. This study focuses on the original rock salt, without considering the effect of creep and relax. The linear elastic assumption may bring some deviation with the actual situation, but it is high needed for figuring out a specific engineering problem.
The applied loads generate displacement field, stress field, and strain field inside for a certain area. The basic equations of small deformation elastic mechanics are as follows: equilibrium differential equation: geometric equation: physical equation: where is the differential operator and is the elastic matrix depending on the elastic modulus of the formation and Poisson’s ratio .
Generally, the finite element formula iswhere is the global stiffness matrix, is the node displacement matrix, is the boundary load array, and is the parameter vector used to describe the boundary loads of UGS reservoir.
Equation (4) is transformed to
Equation (2) is substituted in (3) to get the following equation:
Then, (5) is substituted in (6); that is,
From (7), the nodal stress and the boundary loads satisfy the following equation:where . is the transfer matrix between the nodal stress matrix and the boundary loads array. For the linear elastic body, could be determined when the finite element model of UGS reservoir is established. is uncertain because the loads are unknown, and only a part of is certain, since the stress values of only a few wells in the reservoir to be studied are known. Therefore, in general, the boundary loads parameter vector cannot be calculated from (8) and is waiting to be discovered.
3. Optimisation Method
In situ stress field inversion is used to find the model or parameters that can fit actual data under certain principles. It is not only a method of data fitting, but also a way of optimising. The inversion result is closely related to the optimal criterions, which is used to evaluate the degree of data fitting. Sometimes, we have prior information about some areas that have already been researched. Based on this prior information, we could restrain the optimal criterion and parameters. Generally, the data going to be used in stress inversion calculation is obtained from the reservoir stress field and contains the stress of logging, fracturing, borehole breakout, or core testing of certain wells. The principle of stress inversion (or objective function) could be described as the following models by using the least squares method as the optimal criterion.
(1) Constraint Model. If the stress value of the test well is , that is, the simulation value of the stress field from the finite element calculation is known, the model could be described aswhere is derived from the finite element calculation, is the mean squared error between the measurement and calculation, is the unknown boundary load vector, and and are the boundary load parameter functions of the upper bound and lower bound, respectively, determined by prior information. and are the boundary load parameters of upper and lower bound, respectively.
(2) Unconstraint Model. An unconstraint model exists when the condition in (9) is not met; that is,where symbols are the same as those in (9).
4. Inversion Algorithm
According to the principle of inversion calculation, back analysis is used to calculate the input data of a system by analysing its output data [13]. Therefore, if the system is well constructed, it is possible to obtain the input data from its output data, which is usually achieved through experiments. The back analysis ensures that a unique solution can be generated, as long as the mechanical model is real and the data is sufficient.
The general way of inversion is to present the parameter vector of boundary load and then calculate the reservoir’s stress response value under . During this process, the finite element method is used to approach the actual data of some well sites based on the least squares method. Usually, the Gauss-Newton method is adopted to solve the least squares problem of (9) or (10). However, ill-conditioned coefficient matrix of the equation, low iterative speed, and spurious convergence may occur sometimes. In this paper, a damped least squares method [Levenberg-Marquardt (LM) method] [14] is implemented to approach the parameter vector of boundary load . The LM method, which has a high optimisation speed and precision, possesses both the advantages of the gradient method and the Newton method. When the initial value of damping factor is small, the step size is set to be equal to that of the Newton method, while it is equal to that of the gradient descent when is large.
When the iteration reaches step , the column vector of stress response and boundary load and are set. Then, series expansion around is performed and the trace above the second-order polynomial is ignored to get the following equation [15]:where is the Jacobi matrix of in , and it is defined from the finite difference method aswhere could be calculated from (7).
The search direction used in LM method is from the solution of a set of linear equalities. The iterative formula iswhere is the search direction vector of and is the unit matrix. Scalar is the damping factor controlling the direction and value of . When is 0, the direction of is the same as the result of the Gauss-Newton method. When tends to infinity, we could derive the following equation from (13):Thus,where is in the direction of the zero vector and has a steep descending slope for , which is large enough; that is,From (13), we derivedThen we got
The boundary load vector could be derived from iterative solving of (13)–(18). The key step is how to choose iterative initial value. If the range of the boundary load is known, we can set the iterative initial value of as the average of the upper and lower bounds. If the range is an unbounded domain, the iterative initial value of could be determined after experiential trials.
From the above discussion, if the measured data of reservoir stress is available, the boundary load could be determined. The load parameter vector and reservoir’s horizontal stress could satisfywhere denotes the transfer matrix, is the point number for the measurement, and is the dimension of the boundary load vector. If , a set of boundary load vectors could be clear based on the principle of controlling the comprehensive error to be minimum for the nonuniqueness of the solution.
The specific steps of finite element inversion of the in situ stress field go as follows: Firstly, choose a proper objective function. Secondly, adjust and search parameters by the optimal method and forward-model the computed value by the finite element numerical algorithm. Thirdly, substitute the measured and calculated values into the objective function and judge the result from the objective function. If the result did not reach the minimum, the adjusting and searching should be continued. Repeat this process until the result from objective function reaches the minimum value and finally, the boundary load parameter is generated. Afterwards, the forward modelling is carried out to acquire the reservoir’s stress field by the boundary load derived from the inversion. (Figure 1 shows the inversion process.) The applied boundary condition is shown in Figure 2. The calculation is performed in the entire circle area to avoid stress concentration. Assuming the strata are boundless and stress is continuously changed, if the stress applied on the circle can make the inversion results close enough to the test results, the stress applied on the rectangular must be right. The stress distribution in the object area is the result we want from the inversion.
5. Example and Application
Based on the aforementioned theory, in order to get the boundary load of the in situ stress field, we could program the finite element iteration by finite element software. In this paper, a Turkey salt cavern is chosen to illustrate the application of reversion technology on resolving in situ stress boundary load. Furtherly, this method is also applied on the in situ stress analysis in Qinshui basin of China [16–18]. The least moving squares method can also be used to deduce the in situ stress field by carefully selecting the primary function and the weight function [19]. The UGS distribution is shown in Figures 3–5. Figure 3 shows the well depths of multiple UGS wells. The minimum depth is 1110 m, while the maximum depth is 1500 m. Figure 4 shows the three-dimensional distribution of these UGS wells. Figure 5 shows the planar structure of the salt cavern. Twelve reservoirs, as many as UGS wells, construct the salt cavern. All sites are in the southeast section, which is indicated by the blue spots. The red line in Figure 5 indicates the fault zone of the research area.
5.1. Building UGS in a Salt Cavern
Salt owns several properties that make it ideal for gas storage [19]. It has a moderately high strength and flows plastically to close fracture that could otherwise become a leak. The porosity and permeability of salt to liquid and gaseous hydrocarbons are close to zero, thereby preventing stored gas escape [20, 21]. Salt caverns can provide high deliverability, and gas can be withdrawn quickly because no pressure from flow is lost through the pores. A salt cavern in south Turkey is used to construct underground natural gas storage to compensate season-related fluctuations of gas consumption. Moreover, by a salt cavern, gas can be provided to customers under optimum economic conditions between buying and selling. To construct UGS by a salt cavern, a well should be drilled first and then, water should be injected to dissolve the salt. Brine is pumped out of the salt cavern, and a cavity is formed which can be treated to store natural gas.
5.2. Calculation Model and Boundary Conditions
Based on the prerequisite of the numerical simulation, a geologic model, the target layer isolation body of a rock block is taken as the simulation object. Two rules should be followed to determine the computational area: (1) extend the geometric range to eliminate or weaken the effects from the boundary conditions on the research area and (2) the geometry constraint conditions on the boundary are easy to get. In general, the layered distribution of underground sedimentary rock could be treated as a plane problem of elastic mechanics in the comprehensive mathematical model. Figure 5 is the mechanical model of the research area, which is a rectangular area with a length of 8 km and width of 3.2 km. Three major fault zones are taken into consideration. The lengths and widths of fault zones are 8 km 1 km, 12 km 0.8 km, and 3.5 km 0.8 km, respectively.
The following eight factors are taken as the basic factors of tectonic force: compressive or tensile forces and shear forces in the east, west, south, and north. The boundary load is adjusted depending on the size of research area. As the stress field is continuously changed, the external stress in the circle can be divided into 12 parts, which is large enough to guarantee the accuracy of calculation. By changing the value of the 12 boundary forces, we could simulate the compressive or tensile forces in the direction of east–west and south–north, as well as shear deformation in the horizontal plane to get the in situ stress field distribution of the partial research area. At the same time, boundary force constraint is applied to make the value of the resultant force and the resultant moment zero. (Figure 6 shows the method to apply the boundary force.) The edge points of the lower boundary of the model are fixed to ensure no rigid translation and rotation in the research area.
Generally, a geologic body is taken as a rock with a uniform block to confirm the rock physical parameters from different tectonic positions. The basic assumption is that each fault zone is homogeneous. The rock strength, Young’s modulus, and other rock mechanical parameters should be decreased according to the rock properties in specific depth. (Table 1 shows information about mechanical parameters chosen by the model.) The cored sample is taken from the spot field and by indoor experiment, elasticity modulus and Poisson’s ratio are obtained. The grids of the model adopted mainly are four-node isoperimetric rectangular elements and for part of the area, triangle elements are used. In addition, the element length of the fault zone and its surrounding area is shorter and the grids are relatively finer than those of other parts.
5.3. Inversion of the In Situ Stress Field of the Salt Cavern
(1)Initial boundary conditions: the initial values of the distributed normal stress in the east, south, west, and north could be set as an arbitrary value with 10 MPa here. (Figure 6 shows how the boundary force acts.)(2)Constraint object: based on the reference about in situ stress measurement, we chose 10 values of in situ stress and five directions from eight survey points as known constraint parameters to inverse 12 unknown boundary forces. (Figure 5 shows the distribution of constraint points, and measured data are shown in Tables 2 and 3.)(3)Searching range: the range of the boundary normal stress constraint value is 1 MPa~30 MPa (lower bound~upper bound) and in the same way the range of the boundary shear stress constraint value −10 MPa~10 MPa. The precision of the iteration termination is 10−5 MPa2. The boundary force constraint should be taken into consideration to make sure that the resultant moment of the boundary force is zero. The relevant constrained inversion model could be derived from (9).(4)Inversion result: based on the conditions above, the calculation results on the survey points t from inversion are shown in Tables 2 and 3, where and are the maximum and minimum principal stress magnitude, respectively. The measurement results from well logging are also presented and compared with the inversion results. Table 4 shows the results of the boundary load inversion.
Tables 2 and 3 show that the maximum relative errors of stress components and are controlled within 8%. The errors of the direction of maximum principal stress are controlled within 4°. All the stress values obtained from the inversion reached or are close to the actual data, and the precision of the inversion results could meet engineering requirements.
We applied the optimised boundary force from the inversion to perform forward calculation. This approach could lead to the distribution of the stress field value and direction of the Turkey salt cavern (shown in Figures 7 and 8).
6. Stress Distribution Character of the Well near the Hole of the UGS
On the basis of the analysis above, the two principal stresses and of the 12 wells are determined. Given that the stress concentration occurs by the hole, which will influence the stability of the shaft wall, the von Mises stress distribution character near the shaft wall (0–2.55 m) should be checked up, as shown in Figure 9. Here well UGS1 is taken as an example, with well diameter of UGS1 27.9 cm, bottom hole pressure 12 MPa, overburden pressure 23.1 MPa, and well depth 1100 m.
(a) 0.00 m around the shaft wall
(b) 0.60 m around the shaft wall
(c) 1.05 m around the shaft wall
(d) 1.65 m around the shaft wall
(e) 2.55 m around the shaft wall
As shown in Figure 9, given that stress nonhomogeneous with the in situ stress exists, the von Mises stress distribution of the shaft wall is in an oval shape. The von Mises stress has a higher value in the shaft wall, which is within 28.1–56.8 MPa. The von Mises stress decreases rapidly away from the shaft wall and the shape changes to a butterfly-like one. The maximum von Mises stress decreases to 17.1 MPa when 0.6 m away from the shaft wall. It is worth mentioning that the result here contains the influence from stress concentration around the shaft wall. As shown in Figures 9(d) and 9(e), the von Mises stress owns the same distribution pattern and tends to be steady far from the shaft wall.
As shown in Figure 10, just like the von Mises stress, the shearing stress distribution is in an oval shape near the shaft wall and changes to a butterfly-like shape from the shaft wall. The shearing stress on the shaft wall is maximum (29.1 MPa) due to the stress concentration and it decreases rapidly away from the shaft wall with 8.9 MPa when 0.6 m away. Shearing stress tends to be stable when it is 1.6 m away from the shaft wall.
(a) 0.0 m around the shaft wall
(b) 0.60 m around the shaft wall
(c) 1.05 m around the shaft wall
(d) 1.65 m around the shaft wall
(e) 2.55 m around the shaft wall
From the stress analysis above, we can draw the conclusion as the following.
The stress concentration phenomenon does exist around the hole. The von Mises stress and shearing stress are maximum on the shaft wall. The stress distribution is in an oval shape around the shaft wall while it is in a butterfly-like shape near the hole. Both stresses decrease rapidly when they are 1 m away from the shaft wall and tend to be stable when 1.3 m away from the shaft wall.
7. The Hole Stability Analysis of UGS
7.1. The Rock Fracture Analysis around the Hole of UGS Based on the Griffith Principle
Griffith defined the inherent defects or micro-fracture in the material as the Griffith fracture. Based on the Griffith criterion [22], the rock fracture condition in the planar stress state is shown below.
When , the fracture criterion is
Equation (20) can also be like , .
When , the fracture criterion iswhere is the tensile stress from the rock unidirectional tensile test. is the angle between the long-axis of the oval fracture and the principal compressive stress axis. The fracture angle is .
If a new fracture occurs, it will be in the direction of oval fracture. The angle between the new fracture and the original one is (Figure 11). The following equation can be obtained by the fracture criterion:
Equation (22) shows that the angle between the new fracture and the original one is in the condition of . The negative sign indicates the direction of the new fracture is clockwise rotation by from the principal compression stress axis (Figure 11).
For the three-dimensional problem, Griffith fracture criterion can be depicted as follows.
When , the fracture criterion is
When , the fracture criterion is
Based on the fracture criterion above, it can be judged only whether the tensile fracture or the shear fracture occurs while the developing degrees of the fracture cannot be determined. In fact, with the asymmetry of the tectonic stress distribution, the fracture developing degrees of tectonic part differ. The concept of fracture factor is introduced to judge the fracture developing degrees of the rock quantitatively. The fracture factor is given bywhere and are tensile fracture factor and shear fracture factor. and are tensile stress and shearing stress in the fracture plane. and are tensile strength and shearing resistance. The shaft rock will not be cracked when the fracture factor is below 1.
7.2. Rock Fracture Analysis Results in the UGS
According to the Griffith fracture criterion presented above, we conduct the tensile fracture analysis on the rock of the sidewall in UGS of a Turkey salt cavern. With the obtained maximum and minimum horizontal in situ stresses of the test wells, the tensile stress and shearing stress can be calculated by the classical elastic theory. The calculation model is shown in Figure 12, where are the maximum and minimum horizontal principal in situ stresses, is the radial normal stress, is the tangential normal stress, and is the shearing stress.
The results are shown in Table 5, showing that the maximum tensile stress is smaller than the tensile strength. Therefore, on the rock of the shaft wall, no fracture will be developed.
Moreover, according to the Coulomb-Navier failure criterion, the shear fracture analysis is conducted on the rock of the shaft wall in UGS of a Turkey salt cavern. Table 6 shows the analysis results. The shear fracture factors, , in the hole are all under 1. Given that the maximum shearing stress is less than the shearing strength, the rock of the shaft wall will not develop a shear fracture.
8. Conclusion
In this paper, we improved the nonuniqueness of the in situ stress field inversion solution by a priori constraints. From the perspective of basic equations of elastic mechanics, an in situ stress numerical inversion model is built and basic factors of in situ stress field are considered in applying the finite element inversion method on the boundary loads of the in situ stress field. An appropriate geometrical model and suitable boundary conditions ensured an accurate in situ stress field is derived. The method is successfully applied in UGS of a Turkey salt cavern. The optimisation results on principal stress are close to those from measurement. By the optimisation, the overall in situ stress field of the research area is obtained and the stress distribution of the shaft wall is presented according to the calculation results. The safety of the shaft wall is confirmed by applying the corresponding fracture criterion.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors wish to acknowledge the financial support of the National Natural Science Foundation of China (Grants nos. 51374228, 51274231, and U1262208) and the CNPC Key Laboratory Fund (Contract no. 2014A-4214).