Data-Driven Dynamic Inversion Method for Complex Fractures in Unconventional Reservoirs

Abstract

Hydraulic fracturing is a crucial technology for enhancing the recovery of oil and gas from unconventional reservoirs. Accurately describing fracture morphology is essential for accurately predicting production dynamics. This article proposes a new fracture inversion model based on dynamic data-driven methods, which is different from the conventional linear elastic fracture mechanics model. This method eliminates the need to consider complex mechanical mechanisms, resulting in faster simulation speeds. In the model, the fracture morphology is constrained by combining microseismic data and fracturing construction data, and the fracture tip propagation domain is introduced to characterize the multi-directionality of fracture propagation. The simulated fracture exhibits a multi-branch fracture network morphology, aligning more closely with geological understanding. In addition, the influence of microseismic signal intensity on the direction of fracture propagation is considered in this study. The general stochastic approximation (GSA) algorithm is employed to optimize the direction of fracture propagation. The proposed method is applied to both the single-stage fracturing model and the whole well fracturing model. The research findings indicate that in the single-stage fracturing model, the inverted fracture morphology aligns closely with the microseismic data, with a fitting rate of the fracturing construction curve exceeding 95%, and a microseismic data fitting rate exceeding 93%. In the whole well fracturing model, a total of 18 sections were inverted. The fitting rate between the overall fracture morphology and the microseismic data reached 90%. The simulation only took 5 minutes, demonstrating high computational efficiency and meeting the needs of large-scale engineering fracture simulation. This method can effectively support geological modeling and production dynamic prediction.The world has abundant shale gas reservoir resources; however, due to the influence of reservoir rock properties, its development poses significant challenges [1-4]. Hydraulic fracturing technology can effectively enhance the physical properties of reservoirs and form complex fracture networks within the reservoir, thereby promoting oil and gas production [5-7]. In order to assess the development impact of shale gas reservoirs and devise appropriate development plans, it is necessary to establish a numerical model that is specific to the shale gas reservoir in question. Accurately describing the post-fracturing fracture morphology is crucial for model construction and subsequent flow simulations, as it is a key factor in ensuring the accuracy of model calculation results [8]. Moreover, the morphology of fractures post-fracturing is often highly complex, characterized by a network structure of fractures [9, 10]. Many existing fracture propagation models only consider a simplified quasi-three-dimensional or three-dimensional straight fracture structure. However, these models fail to adequately explain the real distribution of fractures observed during hydraulic fracturing operations [11]. To achieve accurate simulation of fracture morphology, it is crucial to find a complex fracture inversion method that can effectively combine the constraints of fracturing monitoring data for synchronous inversion of fracture morphology.The microseismic data points are dynamically monitored fracture data during the hydraulic fracturing process and provide a relatively accurate reflection of the distribution of underground fractures [12]. While microseismic data points alone may not provide a clear representation of fracture morphology continuity, integrating numerical simulation methods for fracture propagation with microseismic data can be a feasible approach to determine the direction of fractures [13]. Existing methods for hydraulic fracturing fracture propagation mainly include finite element methods, extended finite element methods, boundary element methods, and unconventional fracture propagation models [14-16]. However, these methods assume a single direction for fracture propagation, determined by criteria such as the maximum circumferential tensile stress criterion or the Griffith criterion [16]. As a result, these methods may not be suitable for accurately simulating the complex simulating the morphologies of multi-branch fracture networks when constrained by microseismic data points.Some scholars have explored the use of fractal theory in constructing a multi-branch fracture network structure and fitting microseismic data, which has shown significant progress [17-20]. Sheng et al. utilized fractal theory to establish a complex fracture model for fitting microseismic data. They employed fractal dimension to characterize parameters such as fracture porosity, permeability, and compressibility [17]. Cui used a tree structure generated by an L-system, similar to fractal theory, as an approximately to represent the fracture morphology and fit the microseismic data [21]. Fractal methods offer the advantage of creating complex multi-branch structures that can fit widely dispersed microseismic data. However, these methods tend to construct fracture structures with high similarity and may not fully capture the diversity of complex fracture structures.Furthermore, scholars have developed the discrete fracture network model (DFN) as an alternative approach. Unlike fractal methods, DFN models do not display the branch structure. Instead, they utilize multiple intersecting long straight fractures to fit microseismic data [22-25]. The DFN method is widely applied due to its simplicity and ease of integration with flow simulation methods. Warpinski et al. proposed a box model based on the DFN method to approximate the distribution of microseismic data and evaluated the stimulated reservoir volume [13]. Many scholars have improved this method [24, 25]. However, whether using fractal methods or DFN models, they have not taken into account the influence of geological parameters and microseismic intensity changes on fracture structures. In the hydraulic fracturing process, complex fractures are formed in the formation by fracturing fluid, so the conservation of fracture space and the injected volume of fracturing fluid should be considered. In these methods, volume conservation is often overlooked. Therefore, to accurately simulate the distribution of complex fractures, it is essential to not only constrain the model with microseismic data points but also incorporate actual hydraulic fracturing data, such as pump pressure parameters and injection volume data, to jointly constrain the fracture morphology.In this article, based on the method proposed by Zhao et al. [26], the concept of extended domain is introduced to characterize the multi-directionality of fracture propagation. In addition, the model considers the influence of microseismic signal intensity on the direction of fracture propagation and uses microseismic data to constrain the direction of fracture propagation. Furthermore, the model considers the volume conservation of fractures and injection volume, utilizing hydraulic fracturing construction data to constrain the number of fracture nodes. By coupling these two methods through the local fracture morphology, the model achieves dynamic inversion of complex fracture networks. The method proposed in this article takes into account the constraints by microseismic data points, which are assumed to accurately reflect the actual fracture distribution in the reservoir. However, it should be noted that the model does not currently consider the propagation of fractures in the height direction, assuming a constant height value and no fluid loss in the fracture. Therefore, future research is required to investigate the impact of fluid loss and the propagation of fractures in the height direction.In this section, drawing inspiration from the complex fracture description method proposed by Zhao et al. [26], the stochastic characteristics of fracture propagation are considered. Building upon this, the anisotropy of reservoir geological conditions is taken into account, and dynamic geological parameters are obtained by introducing well-logging data, enabling the characterization of reservoir heterogeneity.To obtain the mechanical parameters of the reservoir, it is necessary to use well-logging data to obtain specific parameters such as shear-wave travel time, compressional-wave travel time, and reservoir density at particular locations, thereby calculating the formation stress. Shear-wave travel time and compressional-wave travel time can be obtained through dipole shear-wave logging, while reservoir density is acquired from density logging data [27]. Therefore, the Poisson’s ratio and Young’s modulus in the reservoir can be expressed as:where Ed,I is the elastic modulus of the reservoir, GPa; νd,i i is the Poisson’s ratio ; ξi and ςi are defined as:where ρb,i is reservoir density, kg/m3; Δts is the shear wave travel time, s; Δtc is the compressional wave travel time, s.Given the obtained Poisson’s ratio and elastic modulus data, considering the rock as a porous elastic medium, the stress parameters of the rock can be determined based on linear elastic fracture mechanics theory, represented as:where σh is the minimum horizontal principal stress, MPa; σH is the maximum horizontal principal stress, MPa; σv is vertical stress, MPa; ps is pore pressure, MPa; α is the effective stress coefficient; αT is the coefficient of linear expansion; ΔT is the formation temperature difference, K; ζH and ζh are the strain coefficients of the maximum principal stress and the minimum principal stress, respectively.The induced stress can be solved using the analytical solution proposed by Green and Sneddon [28]. Employing the superposition principle allows the determination of the induced stress components for the global fracture, which can be expressed as:where σx,ij is the x-axis induced stress component, MPa; σy,ij is the y-axis induced stress component, MPa; τxy,ij is the shear stress component, MPa; N is the number of fractures; pnet is the net pressure of fracture, MPa; a is the half length of the fracture, m; lm is the distance between the space point (i, j) and the midpoint of the fracture; l1,m, l2,m are the distance between the space point (i, j) and the two ends of the fracture, m; θm is the distance between the space point (i, j) and the midpoint of the fracture; and θ1,m, θ2,m are the angles between the spatial point (i, j) and the two ends of the fracture, respectively.Considering the multi-directionality of fracture initiation, unlike traditional fracture criteria, within a propagation domain Γ at the fracture tip, fractures will initiate, as shown in Figure 1. The size of the fracture tip propagation domain is related to the fracture toughness of the reservoir and the tip stress intensity factor. In the model presented in this article, propagation occurs when the circumferential stress in the polar coordinate system at the fracture tip exceeds the critical initiation stress. Therefore, the range covered by the propagation domain Γ in the polar coordinate system with the fracture tip as the origin is expressed as:where r is the polar radius of the propagation domain Γ, m; θ is the polar angle of the propagation domain Γ, rad; K1 and K2 represent the fracture stress intensity factors of type I and type II, respectively, MPa m0.5; σcr is the critical fracture stress strength, MPa; and σfr is the residual stress strength, MPa.In this section, the evolution process and monitoring intensity of microseismic data are considered. The fracture orientation is treated as a variable, and a model for fitting the fracture morphology to microseismic data points is established. The basic idea of the model is to establish a function of the distance between fracture points (Nf) and microseismic points (Mf) based on error estimation, with the microseismic point intensity as the weight. The mathematical model can be represented as:where Ei and Ei−1 are the error metrics of microseismic and fracture in different stages; χ and ψ are the weight coefficients of microseismic fitting values in stage i − 1 and stage i, χ = (i − 1)/i, ψ = 1/i; ωi is the microseismic intensity weighting factor; Іi is the micro-seismic intensity coefficient; m is the number of fracture micro-elements in stage i; n is the number of micro-seismic points in stage i; and Nf,i and Mf,i are the location of stage fractures and microseismic data points, respectively.The microseismic data fitting model in equation (7), can be considered as an optimization problem, where the propagation direction of the fracture at the ith stage is treated as the independent variable, aiming to minimize the error value. This can be expressed as:where J is the objective function; u is the independent variable and represents the fracture direction in this article.For this optimization problem, the independent variable cannot be expressed with a specific function, making conventional gradient-based algorithms unsuitable for this model. In this section, the general stochastic approximation (GSA) algorithm is employed for optimizing and solving the model. The GSA algorithm is an improvement over the simultaneous perturbation stochastic approximation (SPSA). In the conventional SPSA algorithm, the control variables are randomly perturbed simultaneously to obtain an approximate gradient. This process only involves the objective function value J and ensures that the objective function consistently moves uphill in complex problems. Its approximate gradient is represented as:Then, the control variable under the iteration step l + 1 is:where εl is the perturbation step size; Δl is a vector of Bernoulli distribution satisfying ±1; and λl is the iteration step.There is a drawback in the SPSA algorithm due to the randomness of the approximate gradient, the search direction exhibits high uncertainty. This makes the algorithm prone to local convergence or difficulty in escaping the current optimization solution. Therefore, the GSA algorithm has been partially improved upon the SPSA algorithm, generating N-perturbed control variables around the control variable ul at the lth iteration step. This is expressed as:Therefore, the optimal control variables at the l + 1 iteration step are determined by selecting the optimal values from the N-perturbed variables. This method can effectively solve the difficulty of the SPSA method falling into local optimal solution and obtain the best optimization solution.Hydraulic fracturing construction is a real-time process, and in the absence of proppants, the most crucial parameters for hydraulic fracturing are the injection rate and pump pressure. At a certain time t, when a certain amount of fracturing fluid is injected into the fracture, the fluid at the fracturing perforation will produce a certain pressure due to the fluid-solid coupling effect, which is referred to as the bottom hole perforation flow pressure pwp in this article. At the same time, the fracturing pump at the wellhead will provide pumping pressure to the fluid. Ignoring the flow pattern of the fluid in the wellbore, assuming that the fluid in the wellbore is filled, the bottom hole pump injection flow pressure pwb can be calculated by pump pressure, fracturing fluid properties, well depth, and other parameters. The bottom hole perforation flow pressure pwp and the bottom hole pump injection flow pressure pwb can be converted, expressed as:where ppf is the perforation friction of the fluid, MPa; and pcf is the frictional resistance of the fluid along the way, MPa.In equation (12), the calculation of the bottom hole perforation flow pressure pwp is directly related to the fracture morphology and injection rate data. On the other hand, the bottom hole pump injection flow pressure pwb can be directly calculated based on on-site pump pressure data. A mathematical model to constrain the fracture morphology through the injection rate and pump pressure will be established.Under the conditions of known pump pressure and fracturing fluid properties, the bottom hole pump injection flow pressure can be expressed as:where pb is the pump pressure at time t + dt, MPa; ρ is the density of fracturing fluid, kg/m3; and g is the local gravity acceleration, m/s2.For the frictional pressure drop along the fluid flow and the perforation hole frictional pressure drop, the methods established by Crump and Conway through physical experiments can be used for calculation, and are expressed as:where np represents the number of perforations; dp represents the perforation diameter, m; kd represents the orifice flow coefficient, which is between 0.5 and 0.95, the flow coefficient of the perfect hole is equal to 0.5 μ is fluid viscosity, mPa s; D is the wellbore diameter, m; xj is the distance between the perforation point and the bottom hole, m; qin is the construction displacement, m3/min; and qk is the flow rate of the kth perforation, m3/min.Assuming that at time t, the fracture morphology is the constrained ideal fracture morphology, the bottom hole perforation flow pressure at the perforation point can be calculated using the injection rate data at time t + dt and the current fracture morphology. It can be expressed as:where n is the rheological index of the fluid, when the fluid is Newtonian fluid, n = 1; k is the viscosity coefficient of the fluid, mPa s; E is Young’s modulus, GPa; H is the height of the fracture, m; and L is the length of the fracture, m.In this model, only the length of the fracture needs to be adjusted to complete the fitting of the bottom hole perforation flow pressure and the bottom hole pump injection flow pressure. When the two are equal, it is considered that the fracture morphology at the current time is the real formation fracture morphology.In the previous sections, methods for constraining complex fractures using microseismic data points and hydraulic fracturing construction parameters are respectively proposed. In the computations, these two methods need to be coupled, with the fracture propagation direction and morphology serving as the bridge between them. In the model, we make several assumptions: (1) In the microseismic constraint part, the microseismic data points can truly reflect the distribution of reservoir fractures. (2) In the hydraulic fracturing construction parameter constraint part, due to the low permeability of shale reservoirs, fluid loss is not considered. (3) The fracture height considered in the model is a fixed value, which can be set as the reservoir height.The simulation calculation process is shown in Figure 2. The specific steps are as follows: (a) Under given geological conditions, establish the initial model, including properties, mechanical parameters, and perforation parameters. (b) Assume that at time t, the fracture morphology represents the constrained ideal fracture morphology. Utilize pump pressure at time t + dt to calculate the bottom hole pump injection flow pressure. (c) Accounting for fluid resistance effects, calculate the bottom hole perforation flow pressure using injection rate data at time t + dt and the fracture morphology at the current time step. (d) Based on injection rate and wellbore data, calculate the fluid perforation friction and fluid along-path friction. (e) Utilize equation (12) for fitting calculation. If the error is less than the limit error, proceed to the next iteration step or output the fracture morphology. (f) If not fitted, increase the number of fracture nodes, and use microseismic data points to constrain the fracture propagation direction. Finally, output the new fracture morphology. (g) Repeat steps (b) to (f) until the final fracture morphology is output.To validate the fitting performance of the model proposed in this article, we take the single-stage fracturing model of the actual block model as an example. A fracture inversion model with simultaneous constraints from microseismic data and hydraulic fracturing construction data is established, and the fracture morphology is dynamically simulated. The relevant simulation parameters are presented in Table 1.The microseismic data utilized in the model are illustrated in Figure 3(a). The actual hydraulic fracturing construction curve is depicted in Figure 3(b). It can be seen from the microseismic data that the half length of the fracture is about 98 m, the fracture width is about 38 m. The fracture distribution is relatively uniform, and there is no abnormal complex area. Simultaneously, to obtain the coordinates of each microseismic data point, GetData software is used for extraction. The microseismic intensity coefficient І is considered, taking into account the radius proportion of each microseismic data point. In processing the construction curve, the paper directly considered the fracture propagation segment, disregarding the process of wellbore pressure build-up.The real-time inversion results of the fracture are illustrated in Figure 4. The figures depict the fracture morphology at different time intervals: 20, 45, 60, 85, 105, and 135 minutes. It is evident from the images that at each stage, the fitting rate of the pump pressure is high, exceeding 95%. Simultaneously, the fitting rate of the fracture morphology with microseismic data is also high, reaching around 93%. The fracture length of the model inversion is 95 m, and the fracture width is 35 m, which is 3 m different from the microseismic results, which is deemed acceptable in engineering applications. Furthermore, the fracture morphology adequately explains the distribution range of microseismic data and exhibits complex fracture morphology, which is challenging to achieve with conventional numerical simulation methods. These results strongly demonstrate that the proposed model in this article can effectively meet engineering requirements.The method is applied to the single well model of the actual block, and the well W1 of an oilfield is selected as the research object. The well is fractured in the early stage of development, and microseismic monitoring is carried out. The microseismic monitoring results are shown in Figure 5. A total of 18 sections are fractured in well W1, with an average length of about 60 m. The microseismic monitoring results show that the average fracture length of well W1 is 196 m and the fracture width is 45 m. According to the microseismic data and fracturing construction data of well W1, the fracture morphology inversion is carried out. The simulation parameters of well W1 are shown in Table 2.The results of the fracture morphology are shown in Figures 6 and 7 . Figure 6 displays the fracture morphology with microseismic points, indicating a high degree of matching with the microseismic data, with a fitting rate of around 90%. Most of the errors are concentrated near the wellbore, potentially due to significant influences from the wellbore boundary, which are challenging to consider in the model. Figure 7 provides a clear view of the fracture structure, with an average inverted fracture length of 189 m. The deviation from the microseismic data is around 7 m, providing a reasonable interpretation of the microseismic monitoring results. However, it is essential to note that in actual fracturing monitoring processes, microseismic data points are dynamic. In this model, static microseismic data points are considered, which can significantly impact the fracture direction. In future research, incorporating dynamic microseismic data and fitting them with the fracture propagation process will be explored to achieve a higher precision in the fracture morphology.In this article, a multi-objective constrained fracture inversion model based on microseismic data and hydraulic fracturing construction data is established. Conceptual and practical block cases are used to verify the reliability and accuracy of the model calculations. The main conclusions are as follows:A multi-objective constrained fracture inversion model based on microseismic data and hydraulic fracturing construction data is constructed. The model considers the multi-directionality of fracture propagation and the influence of microseismic signal strength on the direction of fracture propagation. Microseismic data is used to constrain the direction of fracture propagation, and hydraulic fracturing construction data is used to constrain the number of fracture nodes. The dynamic inversion of fracture morphology is achieved by coupling these two constraint methods.A single-stage single-cluster fracture propagation model is established to simulate the fitted morphology of fractures at different times. The model results show that the fitting rate of the construction curve reaches 95%, and the microseismic fitting rate is around 93%.Fracture fitting is conducted for an 18-stage well, and the average length of the fractures differed by approximately 7 m compared to the microseismic data. The fitting rate is around 90%, meeting engineering requirements and providing technical support for subsequent geological modeling and numerical simulations.A multi-objective constrained fracture inversion model based on microseismic data and hydraulic fracturing construction data is constructed. The model considers the multi-directionality of fracture propagation and the influence of microseismic signal strength on the direction of fracture propagation. Microseismic data is used to constrain the direction of fracture propagation, and hydraulic fracturing construction data is used to constrain the number of fracture nodes. The dynamic inversion of fracture morphology is achieved by coupling these two constraint methods.A single-stage single-cluster fracture propagation model is established to simulate the fitted morphology of fractures at different times. The model results show that the fitting rate of the construction curve reaches 95%, and the microseismic fitting rate is around 93%.Fracture fitting is conducted for an 18-stage well, and the average length of the fractures differed by approximately 7 m compared to the microseismic data. The fitting rate is around 90%, meeting engineering requirements and providing technical support for subsequent geological modeling and numerical simulations.The data involved in this paper can be obtained in the manuscript.This article has not been submitted elsewhere for publication, in whole or in part, and all the listed authors have approved submission of the manuscript. We declare that there are no conflicts of interest.This study was supported by the 13th Five-Year National Science and Technology MajorProject, grant number 2016ZX05061-009.