非常规储层复杂裂缝的数据驱动动态反演方法

IF 1.8 4区 地球科学 Q3 GEOCHEMISTRY & GEOPHYSICS Lithosphere Pub Date : 2024-01-12 DOI:10.2113/2024/lithosphere_2023_347
Ruixue Jia, Xiaoming Li, Xiaoyong Ma, Liang Zhu, Yangdong Guo, Xiaoping Song, Pingde Wang, Jiantao Wang
{"title":"非常规储层复杂裂缝的数据驱动动态反演方法","authors":"Ruixue Jia, Xiaoming Li, Xiaoyong Ma, Liang Zhu, Yangdong Guo, Xiaoping Song, Pingde Wang, Jiantao Wang","doi":"10.2113/2024/lithosphere_2023_347","DOIUrl":null,"url":null,"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.","PeriodicalId":18147,"journal":{"name":"Lithosphere","volume":"179 1","pages":""},"PeriodicalIF":1.8000,"publicationDate":"2024-01-12","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":"0","resultStr":"{\"title\":\"Data-Driven Dynamic Inversion Method for Complex Fractures in Unconventional Reservoirs\",\"authors\":\"Ruixue Jia, Xiaoming Li, Xiaoyong Ma, Liang Zhu, Yangdong Guo, Xiaoping Song, Pingde Wang, Jiantao Wang\",\"doi\":\"10.2113/2024/lithosphere_2023_347\",\"DOIUrl\":null,\"url\":null,\"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.\",\"PeriodicalId\":18147,\"journal\":{\"name\":\"Lithosphere\",\"volume\":\"179 1\",\"pages\":\"\"},\"PeriodicalIF\":1.8000,\"publicationDate\":\"2024-01-12\",\"publicationTypes\":\"Journal Article\",\"fieldsOfStudy\":null,\"isOpenAccess\":false,\"openAccessPdf\":\"\",\"citationCount\":\"0\",\"resultStr\":null,\"platform\":\"Semanticscholar\",\"paperid\":null,\"PeriodicalName\":\"Lithosphere\",\"FirstCategoryId\":\"89\",\"ListUrlMain\":\"https://doi.org/10.2113/2024/lithosphere_2023_347\",\"RegionNum\":4,\"RegionCategory\":\"地球科学\",\"ArticlePicture\":[],\"TitleCN\":null,\"AbstractTextCN\":null,\"PMCID\":null,\"EPubDate\":\"\",\"PubModel\":\"\",\"JCR\":\"Q3\",\"JCRName\":\"GEOCHEMISTRY & GEOPHYSICS\",\"Score\":null,\"Total\":0}","platform":"Semanticscholar","paperid":null,"PeriodicalName":"Lithosphere","FirstCategoryId":"89","ListUrlMain":"https://doi.org/10.2113/2024/lithosphere_2023_347","RegionNum":4,"RegionCategory":"地球科学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":null,"EPubDate":"","PubModel":"","JCR":"Q3","JCRName":"GEOCHEMISTRY & GEOPHYSICS","Score":null,"Total":0}
引用次数: 0

摘要

(为了验证本文所提模型的拟合性能,我们以实际区块模型的单级压裂模型为例。为了验证本文所提模型的拟合性能,我们以实际区块模型的单段压裂模型为例,建立了一个同时受微地震数据和水力压裂施工数据约束的裂缝反演模型,并对裂缝形态进行了动态模拟。相关模拟参数见表 1。图 3(a)显示了模型中使用的微震数据。实际水力压裂施工曲线如图 3(b)所示。从微震数据可以看出,裂缝半长约为 98 米,裂缝宽度约为 38 米,裂缝分布较为均匀,没有异常复杂区域。同时,为了获得每个微震数据点的坐标,使用 GetData 软件进行提取。考虑到每个微震数据点的半径比例,考虑了微震强度系数І。在处理构造曲线时,本文直接考虑了裂缝的传播段,忽略了井筒压力上升的过程。裂缝的实时反演结果如图 4 所示。图中描述了不同时间间隔(20 分钟、45 分钟、60 分钟、85 分钟、105 分钟和 135 分钟)的裂缝形态。从图中可以明显看出,在每个阶段,泵压的拟合率都很高,超过 95%。同时,裂缝形态与微地震数据的拟合率也很高,达到 93% 左右。模型反演的断裂长度为 95 米,断裂宽度为 35 米,与微震结果相差 3 米,这在工程应用中是可以接受的。此外,断裂形态充分解释了微震数据的分布范围,并表现出复杂的断裂形态,这是传统数值模拟方法难以实现的。这些结果有力地证明了本文提出的模型可以有效地满足工程要求。本文将该方法应用于实际区块的单井模型,并选择某油田的 W1 井作为研究对象。将该方法应用于实际区块的单井模型,选取某油田的 W1 井作为研究对象,在开发初期对该井进行压裂,并进行微震监测。微震监测结果如图 5 所示。根据 W1 井的微震数据和压裂施工数据,进行了裂缝形态反演。W1 井的模拟参数如表 2 所示,压裂形态结果如图 6 和图 7 所示。图 6 显示了带有微地震点的断裂形态,表明与微地震数据高度匹配,拟合率约为 90%。大部分误差集中在井筒附近,可能是由于井筒边界的重大影响,这在模型中很难考虑。图 7 提供了裂缝结构的清晰视图,平均反演裂缝长度为 189 米,与微震数据的偏差约为 7 米,为微震监测结果提供了合理解释。但必须注意的是,在实际压裂监测过程中,微震数据点是动态的。在本模型中,考虑的是静态微震数据点,这会对压裂方向产生重大影响。在未来的研究中,将探索纳入动态微震数据,并将其与裂缝扩展过程拟合,以实现更高精度的裂缝形态。本文建立了基于微震数据和水力压裂施工数据的多目标约束裂缝反演模型。本文建立了基于微地震数据和水力压裂施工数据的多目标约束裂缝反演模型,并利用概念和实际区块案例验证了模型计算的可靠性和准确性。主要结论如下:构建了基于微地震数据和水力压裂施工数据的多目标约束裂缝反演模型。该模型考虑了裂缝传播的多方向性以及微地震信号强度对裂缝传播方向的影响。微地震数据用于约束断裂传播方向,水力压裂施工数据用于约束断裂节点数量。通过这两种约束方法的耦合,实现了对断裂形态的动态反演。
本文章由计算机程序翻译,如有差异,请以英文原文为准。
查看原文
分享 分享
微信好友 朋友圈 QQ好友 复制链接
本刊更多论文
Data-Driven Dynamic Inversion Method for Complex Fractures in Unconventional Reservoirs
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.
求助全文
通过发布文献求助,成功后即可免费获取论文全文。 去求助
来源期刊
Lithosphere
Lithosphere GEOCHEMISTRY & GEOPHYSICS-GEOLOGY
CiteScore
3.80
自引率
16.70%
发文量
284
审稿时长
>12 weeks
期刊介绍: The open access journal will have an expanded scope covering research in all areas of earth, planetary, and environmental sciences, providing a unique publishing choice for authors in the geoscience community.
期刊最新文献
Integrated Simulation for Microseismic Fracture Networks with Automatic History Matching in Tight Oil Development: A Field Case from Block Y2 in Ordos Basin, China Insight into the Evolution of the Eastern Margin of the Wyoming Craton from Complex, Laterally Variable Shear Wave Splitting Re−Os Isotope and PGE Abundance Systematics of Coast Range Ophiolite Peridotites and Chromitite, California: Insights into Fore-Arc Magmatic Processes Indirect Tensile Strength Test on Heterogeneous Rock Using Square Plate Sample with a Circular Hole Complex Segment Linkage Along the Sevier Normal Fault, Southwestern Utah
×
引用
GB/T 7714-2015
复制
MLA
复制
APA
复制
导出至
BibTeX EndNote RefMan NoteFirst NoteExpress
×
×
提示
您的信息不完整,为了账户安全,请先补充。
现在去补充
×
提示
您因"违规操作"
具体请查看互助需知
我知道了
×
提示
现在去查看 取消
×
提示
确定
0
微信
客服QQ
Book学术公众号 扫码关注我们
反馈
×
意见反馈
请填写您的意见或建议
请填写您的手机或邮箱
已复制链接
已复制链接
快去分享给好友吧!
我知道了
×
扫码分享
扫码分享
Book学术官方微信
Book学术文献互助
Book学术文献互助群
群 号:481959085
Book学术
文献互助 智能选刊 最新文献 互助须知 联系我们:info@booksci.cn
Book学术提供免费学术资源搜索服务,方便国内外学者检索中英文文献。致力于提供最便捷和优质的服务体验。
Copyright © 2023 Book学术 All rights reserved.
ghs 京公网安备 11010802042870号 京ICP备2023020795号-1