A numerical simulation method for multi-field coupled heat transfer mechanism in deep carbonate geothermal reservoirs

Abstract

The present invention provides a numerical simulation method for the multi-field coupled heat exchange mechanism of deep carbonate geothermal reservoirs, which relates to the technical field of geothermal development. The present invention improves the real-time monitoring capability of dynamic changes of fractures by using geological three-dimensional models and microseismic monitoring networks; through seismic event positioning algorithms and finite element methods, it can accurately simulate the expansion direction and activity intensity of fractures and evaluate changes in heat flow fields; this method effectively integrates the multi-field coupling effects of heat conduction, fluid flow, solid deformation and chemical reactions, greatly improving the control capability of heat flow in deep carbonate geothermal reservoirs; in addition, by calculating influence coefficients and formulating management strategies, it can significantly reduce potential geological risks caused by geothermal development, optimize the level of geothermal resource development and utilization, and has important scientific value and practical significance for the safe, efficient and economic development of deep carbonate geothermal resources.

Description

Technical Field

[0001] The present invention relates to the field of geothermal development technology, and in particular to a numerical simulation method for the multi-field coupled heat exchange mechanism of a deep carbonate geothermal reservoir. Background Art

[0002] Deep carbonate geothermal reservoirs have abundant heat storage capacity and have become an important target for geothermal resource development. However, traditional geothermal development technologies mainly focus on shallow geothermal systems. Deep carbonate geothermal reservoirs are more difficult to develop due to their complex geological conditions and multi-physics field coupling characteristics. Existing technologies often use simplified models to describe the heat flow behavior of geothermal reservoirs. However, due to the failure to fully consider the expansion of fractures in the formation, chemical reactions and multi-field coupling effects, the model is usually insufficiently accurate and cannot effectively predict and control the changes in the heat flow field during geothermal extraction. This limits the efficient utilization and further development of deep carbonate geothermal resources.

[0003] The existing technology, with publication number CN115310329A, is titled "A Multi-field Coupled Numerical Simulation Method for Geothermal Reservoirs," and includes the following steps: Step 1: Basic geological research; Step 2: Study of thermal reservoir characteristics; Step 3: Construction of a three-dimensional geological model; Step 4: Determination of initial key parameters; Step 5: Establishment of a mathematical model; and Step 6: Optimization of thermal reservoir development and utilization plans. This method, based on a geological model established using the project area's geological structure and drilling data, constructs a multi-field coupled mathematical model for hydrothermal thermal reservoir development. This model simulates the dynamic response of various parameters under large-scale production and irrigation conditions, thereby optimizing the production and irrigation operating parameters for geothermal development in the project area and providing guidance for geothermal development in other regions.

[0004] Existing technologies for deep carbonate geothermal reservoirs have the following main shortcomings: They lack comprehensive consideration of multi-field coupling effects, particularly regarding the coupling of fracture expansion and heat flow, with most studies failing to provide accurate simulation methods. Furthermore, traditional methods for dynamic monitoring of fracture expansion and its impact on geological risk assessment often rely on discrete monitoring data and simple simulation methods, resulting in low reliability and practicality of prediction results.

[0005] The above information disclosed in this Background section is only for enhancement of understanding of the background of the present disclosure and therefore it may contain information that does not form the prior art that is already known to a person of ordinary skill in the art. Summary of the Invention

[0006] The purpose of the present invention is to provide a numerical simulation method for the multi-field coupled heat exchange mechanism of deep carbonate geothermal reservoirs to solve the problems raised in the above background technology.

[0007] To achieve the above object, the present invention provides the following technical solutions:

[0008] A numerical simulation method for multi-field coupled heat transfer mechanism in deep carbonate geothermal reservoirs, comprising the following steps:

[0009] Step S1: obtaining a geological three-dimensional model of the target geothermal area and marking a number of target fractures in the geological three-dimensional model;

[0010] Based on a number of target fractures on the three-dimensional geological model, a microseismic monitoring network is deployed in the target geothermal area to monitor in real time a number of microseismic events occurring in the target geothermal area during geothermal development, and to collect relevant preliminary microseismic event data, including but not limited to time estimation, spatial coordinate estimation point, magnitude estimation value, and focal mechanism estimation information of each microseismic event;

[0011] Step S2: Processing the collected preliminary microseismic event data using an earthquake event location algorithm to accurately determine the time, spatial coordinates, focal mechanism, and magnitude of each microseismic event;

[0012] Based on the spatial coordinates of each microseismic event, determining the associated impact area of ​​each microseismic event and determining the target fracture within the associated impact area;

[0013] Step S3: extracting the seismic phase characteristics of each microseismic event, which include the propagation time, amplitude and frequency of the P wave and the S wave;

[0014] Based on the seismic phase characteristics of each microseismic event, a mathematical model of crack propagation is constructed using the finite element method. This mathematical model is used to analyze the propagation direction and activity intensity of several target cracks in the target geothermal area, thereby obtaining the predicted crack size of each target crack.

[0015] Step S4: obtaining a multi-field coupling model of the target geothermal region, wherein the multi-field coupling model is used to simulate heat conduction, fluid flow, solid deformation, and chemical reaction processes in the deep carbonate geothermal reservoir, and ultimately obtain thermal flow field data for each associated influence area in the target geothermal region;

[0016] Step S5: obtaining the thermal flow field data of each associated influence area in the target geothermal area, and importing the predicted fracture of each target fracture in the associated influence area into the multi-field coupling model for separate simulation, obtaining the thermal flow field change data of each associated influence area, and analyzing the thermal flow field change data of each associated influence area to calculate the influence coefficient. The influence coefficient is used to evaluate the degree of change of the thermal flow field under different fracture extension situations, and finally obtaining the thermal flow field change evaluation result of each associated influence area;

[0017] Step S6: Based on the calculated earthquake magnitude in each associated impact area and the predicted fracture of each target fracture, and through the thermal flow field change assessment results, corresponding management and adjustment strategies are formulated for the geological risk and thermal flow field change status of each target fracture to effectively control the heat flow of the geothermal reservoir.

[0018] Compared with the existing technology, the beneficial effects of the present invention are: by using geological three-dimensional models and microseismic monitoring networks, the real-time monitoring capability of dynamic changes in fractures is improved; through seismic event positioning algorithms and finite element methods, the fracture extension direction and activity intensity can be accurately simulated, and the changes in the thermal flow field can be evaluated; this method effectively integrates the multi-field coupling effects of heat conduction, fluid flow, solid deformation and chemical reactions, greatly improving the control capability of heat flow in deep carbonate geothermal reservoirs; in addition, by calculating the influence coefficient and formulating management strategies, the potential geological risks caused by geothermal development can be significantly reduced, and the level of geothermal resource development and utilization can be optimized, which has important scientific value and practical significance for the safe, efficient and economic development of deep carbonate geothermal resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 Schematic diagram of the overall method of the present invention. DETAILED DESCRIPTION

[0020] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to specific embodiments.

[0021] It should be noted that, unless otherwise defined, the technical or scientific terms used in the present invention should have the usual meanings understood by people with ordinary skills in the field to which the present invention belongs. The "first", "second" and similar words used in the present invention do not indicate any order, quantity or importance, but are only used to distinguish different components. "Include" or "comprise" and similar words mean that the elements or objects appearing before the word include the elements or objects listed after the word and their equivalents, without excluding other elements or objects. "Connect" or "connected" and similar words are not limited to physical or mechanical connections, but may include electrical connections, whether direct or indirect. "Up", "down", "left", "right" and the like are only used to indicate relative position relationships. When the absolute position of the object being described changes, the relative position relationship may also change accordingly.

[0022] Example 1:

[0023] See also Figure 1 , the present invention provides a technical solution:

[0024] A numerical simulation method for the multi-field coupled heat transfer mechanism of deep carbonate geothermal reservoirs is applied to the development process of deep carbonate geothermal reservoirs, and the deep carbonate geothermal reservoirs are used as the target geothermal area. The specific steps include:

[0025] Step S1: obtaining a geological three-dimensional model of the target geothermal area and marking a number of target fractures in the geological three-dimensional model;

[0026] Based on a number of target fractures on the three-dimensional geological model, a microseismic monitoring network is deployed in the target geothermal area to monitor in real time a number of microseismic events occurring in the target geothermal area during geothermal development, and to collect relevant preliminary microseismic event data, including but not limited to time estimation, spatial coordinate estimation point, magnitude estimation value, and focal mechanism estimation information of each microseismic event;

[0027] Step S2: Processing the collected preliminary microseismic event data using an earthquake event location algorithm to accurately determine the time, spatial coordinates, focal mechanism, and magnitude of each microseismic event;

[0028] Based on the spatial coordinates of each microseismic event, determining the associated impact area of ​​each microseismic event and determining the target fracture within the associated impact area;

[0029] Step S3: extracting the phase characteristics of each microseismic event through the determined spatial coordinates, focal mechanism and magnitude of each microseismic event, the phase characteristics including the propagation time, amplitude and frequency of the P wave and S wave;

[0030] Based on the seismic phase characteristics of each microseismic event, a mathematical model of crack propagation is constructed using the finite element method. This mathematical model is used to analyze the propagation direction and activity intensity of several target cracks in the target geothermal area, thereby obtaining the predicted crack size of each target crack.

[0031] Step S4: obtaining a multi-field coupling model of the target geothermal region, wherein the multi-field coupling model is used to simulate heat conduction, fluid flow, solid deformation, and chemical reaction processes in the deep carbonate geothermal reservoir, and ultimately obtain thermal flow field data for each associated influence area in the target geothermal region;

[0032] Step S5: obtaining the thermal flow field data of each associated influence area in the target geothermal area, and importing the predicted fracture of each target fracture in the associated influence area into the multi-field coupling model for separate simulation, obtaining the thermal flow field change data of each associated influence area, and analyzing the thermal flow field change data of each associated influence area to calculate the influence coefficient. The influence coefficient is used to evaluate the degree of change of the thermal flow field under different fracture extension situations, and finally obtaining the thermal flow field change evaluation result of each associated influence area;

[0033] Step S6: Based on the calculated earthquake magnitude in each associated impact area and the predicted fracture of each target fracture, and through the thermal flow field change assessment results, corresponding management and adjustment strategies are formulated for the geological risk and thermal flow field change status of each target fracture to effectively control the heat flow of the geothermal reservoir.

[0034] Further explanation: marking several target fractures in the geological 3D model includes:

[0035] 1.1) Determine the target geothermal area based on the current minimum requirements for heat transfer mechanisms:

[0036] 1.11) Obtain geological exploration data, historical seismic activity data, and current geothermal development status in the current area to define the scope of the target geothermal area. The principles for defining the scope of the target geothermal area are as follows:

[0037] Geophysical exploration: using geological radar, gravity anomaly detection and resistivity imaging to analyze the distribution and boundaries of underground thermal reservoirs;

[0038] GIS technical analysis: Input the collected geological, seismic, and ground temperature information into the Geographic Information System (GIS) and determine the target range through overlay analysis models;

[0039] 1.12) Ensure that the target geothermal area has good heat storage conditions, including heat flow intensity, geothermal temperature, and that the thermal reservoir meets the current minimum requirements for heat exchange mechanism development.

[0040] 1.2) Construct a 3D geological model of the target geothermal area, including the following:

[0041] 1.21) Obtain geological data from geological survey reports, drilling data, seismic data, and geothermal development data for the target geothermal area. Geological data types include, but are not limited to, stratigraphic information, lithologic characteristics, orientation and distribution of each fracture, and rock mechanical properties. Ensure that the spatial resolution of the data meets the modeling requirements below.

[0042] 1.22) For geological 3D models: Select Petrel or Leapfrog geological modeling software, import the collected geological data into the selected geological modeling software, perform data format conversion and verification to ensure data integrity;

[0043] 1.23) Based on the geological data, use the stratigraphic construction tools of the geological modeling software to construct a three-dimensional stratigraphic model to represent the spatial distribution of various rocks and rock layers;

[0044] 1.24) Construct the three-dimensional location of each fracture in the target geothermal area;

[0045] 1.25) Use historical seismic data to verify whether the fracture distribution in the geological 3D model is consistent with the actual situation;

[0046] If there are deviations in the data, parameters are optimized based on historical earthquake activity to adjust the fracture parameters and stratum interfaces in the model. The specific implementation method of parameter optimization is a general technical means in the construction of geological three-dimensional models and will not be described in detail.

[0047] Determine the three-dimensional location of each fracture in the target geothermal area, as well as the fracture characteristics of each fracture, including fracture size, fracture density, and fracture distribution pattern;

[0048] Analyze the fracture characteristics of each fracture in the geological three-dimensional model, including fracture size, fracture density, and fracture distribution pattern;

[0049] The impact of each fracture characteristic on geothermal development is assessed. The total score = size score + density score + distribution score. The larger the assessment value, the greater the corresponding impact. The assessment strategy is as follows:

[0050] Each fracture feature is assigned a score value, and the impact of the fracture on development is determined based on the actual situation. In this embodiment, the scoring rules are set as follows:

[0051] Crack Size Score:

[0052] 1 point: The crack size is very small, with a width of less than 1 mm, a depth of less than 2 m, and a length of less than 10 m.

[0053] 2 points: The crack is small, with a width of less than 1mm, a depth of 2-4m, and a length of 10-20m.

[0054] 3 points: The crack is medium in size, 1-2 mm in width, 4-6 m in depth, and 20-30 m in length.

[0055] 4 points: The crack is relatively large, with a width of 2-3 mm, a depth of 6-8 m, and a length of 30-40 m.

[0056] 5 points: The crack size is very large, width ≥3mm, depth ≥8m, length ≥40m.

[0057] Fracture density score:

[0058] 1 point: The density is very low, with less than 5 cracks per cubic meter.

[0059] 2 points: Low density, with 5-10 cracks per cubic meter.

[0060] 3 points: medium density, with 10-15 cracks per cubic meter.

[0061] 4 points: The density is relatively high, with the number of cracks being 15-20 per cubic meter.

[0062] 5 points: The density is very high, with the number of cracks ≥ 20 / cubic meter.

[0063] Crack distribution pattern score:

[0064] 1 point: The crack distribution is very regular, with more than 80% of the crack area having the same direction, and more than 85% of the crack length having uniform spacing.

[0065] 2 points: The distribution of cracks is relatively regular, the directions of the cracks (70%, 80%) are consistent, and the spacing of the crack lengths (70%, 85%) is uniform.

[0066] 3 points: The cracks are distributed in a moderate (50%, 70%) area with consistent orientation and are evenly spaced (40%, 70%) in length.

[0067] 4 points: The distribution of cracks is relatively irregular. The directions of cracks in the (30% and 50%] regions are consistent. The spacing of cracks in the (20% and 40%] regions is uniform.

[0068] 5 points: The distribution of cracks is very irregular. The cracks are in the same direction in the (0%, 30%) region, and the crack lengths are evenly spaced (0%, 20%).

[0069] Based on these evaluation numerical results, these cracks are prioritized in descending order, and then divided into high-priority cracks, medium-priority cracks, and low-priority cracks according to the size of the total score value, and the cracks that meet the high-priority cracks are selected as target cracks;

[0070] High-priority fractures: with a comprehensive score of 13 or above. These fractures have significant impacts on size, density, and distribution patterns, and have a positive effect on geothermal development.

[0071] Medium-priority cracks: with an overall score of 9 to 12. These cracks have good local characteristics and are suitable as candidate targets.

[0072] Low-priority fissures: Comprehensive scores of 8 or less. These fissures have little impact in all aspects and are generally not preferred. Target fissure selection is based on actual needs, and high-priority fissures are selected as target fissures.

[0073] Assign a number to each target fracture in the geological three-dimensional model, forming {1, 2, ..., i, ..., n}, where i represents the index of the target fracture i and n represents the total number of target fractures; and mark the corresponding fracture characteristics in the corresponding index of each target fracture;

[0074] Deploy a microseismic monitoring network in the target geothermal area, specifically:

[0075] The microseismic monitoring network consists of seismometers and seismic sensors. This embodiment uses three-component seismometers and broadband seismic sensors for deep deployment.

[0076] Deployment plan: Seismometers and seismic sensors are placed in boreholes at different depths within the target geothermal area, and combined with surface and underground equipment to form a multi-layered microseismic monitoring network. Each monitoring location in the microseismic monitoring network is marked as a monitoring point to ensure monitoring of all levels of the fracture system.

[0077] Selection rationale: Densely deployed, high-precision seismic sensor arrays can improve the spatial resolution of microseismic activity detection and are applicable to environments with complex regional geology and multi-level distribution of fracture systems, thereby more accurately reflecting fracture changes.

[0078] Select and test seismic sensors:

[0079] Three-component broadband seismic sensor: Select a broadband sensor that meets geothermal monitoring requirements, such as the GeospaceGS-30 or Geosense3C. These sensors are chosen because they cover a wider frequency range, ensuring sensitivity to small seismic activity.

[0080] Seismograph selection: Choose a high-precision digital seismograph that can accurately capture microseismic events and withstand high temperature and high pressure environments;

[0081] Seismic sensor function test:

[0082] Pre-install seismic sensors in a laboratory environment to test their frequency response, noise rejection, and stability;

[0083] Testing the depth adaptability of seismic sensors to confirm their ability to operate stably in the high-temperature and high-pressure environments of geothermal wells;

[0084] For microseismic monitoring network design:

[0085] Depth distribution design:

[0086] Based on the 3D geological model of the target geothermal area, the high-frequency areas of microseismic activity are identified. Based on the distribution of fractures, the layout of seismic sensors at different depths is designed to ensure that microseismic activity at all levels of the fracture system can be monitored.

[0087] For deep fracture systems, seismic sensors are placed at deeper well locations; for shallow fracture systems, seismic sensors are placed at shallower well locations.

[0088] It is recommended to deploy at least one set of seismic sensors above, below, and on both sides of each target fracture to ensure efficient coverage;

[0089] Surface and underground equipment coordination:

[0090] At least three high-sensitivity seismometers are deployed at the surface installation location as a supplement to ensure sufficient data collection in different spatial dimensions;

[0091] Combining underground seismic sensors with surface monitoring equipment to form a dual monitoring network of surface and underground;

[0092] The equipment is installed as follows:

[0093] For seismic sensor deployment: Use standardized installation tools at the wellhead to install three-component seismic sensors at a specific depth in the well, ensuring that each seismic sensor is stably fixed and connected to the data transmission system;

[0094] After the seismic sensor is installed, a functional test is performed to confirm the working status of each seismic sensor;

[0095] Data transmission system construction:

[0096] Configure a data transmission system to ensure a stable communication link between the underground seismic sensor and the surface data receiving device; in this embodiment, the communication link includes optical fiber and wireless transmission;

[0097] Install data receiving and processing equipment on the surface to ensure real-time reception of data from underground seismic sensors;

[0098] After the microseismic monitoring network is established, preliminary microseismic event data transmitted by seismic sensors are collected in real time. It should be noted that the time estimate, spatial coordinate estimate, magnitude estimate, and focal mechanism estimate collected in the preliminary microseismic event data are all rough measurement results of seismic sensors;

[0099] Preliminary data processing:

[0100] Use data analysis software to preliminarily process real-time data, remove noise and perform filtering;

[0101] Use filtering algorithms or Fourier transform to reduce noise in vibration signals to ensure accurate identification of microseismic events;

[0102] Real-time monitoring system establishment:

[0103] A real-time monitoring system is established at the data receiving center, using visualization tools such as charts or 3D geological displays to display the distribution of current microseismic activity.

[0104] Further explanation: The collected preliminary microseismic event data are processed using an earthquake event location algorithm to accurately determine the time, spatial coordinates, focal mechanism and magnitude of each microseismic event, including:

[0105] The preliminary microseismic event data collected in the microseismic monitoring network also includes: timestamp, arrival time difference and seismic wave information, including:

[0106] The timestamp is the monitoring time point when each microseismic event occurs;

[0107] The arrival time difference represents the propagation time difference from the earthquake source to each monitoring point;

[0108] Seismic wave information is composed of the arrival time and amplitude of P and S waves;

[0109] The data is preprocessed to remove noise, correct errors, and ensure the accuracy of timestamps and epicenter locations. Signal decoding is then performed to determine the preliminary characteristics of the microseismic event, including the type of fault corresponding to the microseismic event and the intensity of the earthquake source.

[0110] Seismic event location algorithms include: using 3D seismic inversion technology to derive and correct preliminary microseismic event data for each microseismic event based on the microseismic monitoring network and 3D geological model, thereby accurately determining the time, spatial coordinates, focal mechanism, and magnitude of each microseismic event;

[0111] The specific steps of 3D seismic inversion technology are as follows:

[0112] 1. Construction of seismic wave propagation velocity model:

[0113] Construct a 3D seismic wave propagation velocity model based on the physical properties of different rock layers in the 3D geological model; the physical properties include rock density, elastic modulus, and P-wave and S-wave velocities;

[0114] The three-dimensional seismic wave propagation velocity model is used to simulate the propagation path of seismic waves in underground media to determine the arrival time difference of seismic waves at each monitoring point;

[0115] 2. Inversion of seismic wave propagation path:

[0116] An inversion algorithm is used to derive the spatial coordinates of each microseismic event, including the epicenter location and focal depth; the inversion algorithm includes the least squares method, genetic algorithm, or simulated annealing algorithm;

[0117] Based on the arrival time difference data and the three-dimensional seismic wave propagation velocity model, the propagation path of the seismic wave is inferred through reverse calculation, and the accurate location of each microseismic event is finally obtained;

[0118] 3. Focal mechanism analysis:

[0119] By analyzing the focal mechanism of microseismic events, we can determine which type of fault slip caused the microseismic event;

[0120] Focal mechanisms include normal faults, reverse faults, and strike-slip faults. They also include non-tectonic focal mechanisms, which include volume expansion, vertical fracture movement, and fluid-induced.

[0121] Isotropic Source: Earthquakes caused by underground volume expansion or compression (such as underground reservoir explosion, steam expansion, etc.). They are generally not accompanied by obvious fault slip.

[0122] For vertical crack movement (Tensile Crack Opening / Closing): It is a micro-earthquake caused by the expansion or closing of vertical cracks, such as volcanic earthquakes caused by magma intrusion leading to crack expansion or micro-earthquakes caused by hydraulic fracturing.

[0123] For fluid-induced seismicity: earthquakes caused by the movement of fissures or faults due to geothermal extraction, oil and gas extraction, waste liquid injection or reservoir induction (such as reservoir earthquakes).

[0124] The fractures generated by normal faults are assumed to be tensile fractures: such fractures are conducive to the flow of thermal fluids and are ideal geothermal reservoir structures.

[0125] The cracks generated by reverse faults are assumed to be compression structures: the permeability of such cracks is poor, which is not conducive to the diffusion of thermal fluids, but there is a certain probability of forming high-pressure areas.

[0126] The waveform inversion method or full-wave inversion technology is used to calculate the direction, sliding mode and stress state of the earthquake source, thereby determining the focal mechanism. This technical means of determining the focal mechanism is common knowledge and will not be described in detail.

[0127] 4. Magnitude calculation:

[0128] The magnitude of each microseismic event is calculated using a seismological formula based on the arrival time and amplitude information of the P and S waves corresponding to the seismic wave information. The seismological formula used is the Richter magnitude formula or the moment magnitude formula. The specific calculation formula of the Richter magnitude formula or the moment magnitude formula is the conventional formula used in existing magnitude calculations and is not described in detail here.

[0129] The Richter scale estimates an earthquake's size by analyzing the amplitude of seismic waves recorded on a seismograph. The larger the amplitude of the seismic waves detected by the seismograph, the greater the energy released by the earthquake. The Richter scale also takes into account the distance from the epicenter to the monitoring point, as seismic waves gradually attenuate as they propagate, with amplitude decreasing with distance. When calculating the Richter scale, the amplitude is corrected for distance, ensuring that earthquakes at different locations and distances can be compared.

[0130] The moment magnitude assesses the size of an earthquake based on the physical properties of the event itself. It considers the area of ​​the fault rupture, the distance of fault slip, and the mechanical properties of the underlying rock. Using these factors, the moment magnitude calculates the total energy released by the earthquake, providing a more accurate reflection of the true intensity of earthquakes of varying sizes.

[0131] The calculation of earthquake magnitude requires combining the distance of the monitoring point and the amplitude of the received P and S waves to quantify the energy released by the microearthquake;

[0132] Derivation and correction of preliminary microseismic event data:

[0133] After obtaining preliminary positioning results through inversion, the time estimate, spatial coordinate estimate, magnitude estimate, and focal mechanism estimate are corrected by combining existing microseismic monitoring network data and geological three-dimensional models. This correction process can also be performed by introducing correction factors, optimization algorithms, or data fusion technology to eliminate errors caused by factors such as seismic wave propagation paths and seismograph accuracy, thereby obtaining more accurate microseismic event location data. The above correction process uses a mature application in existing seismological research and microseismic positioning technology and will not be elaborated on here.

[0134] After obtaining the spatial coordinates of each microseismic event determined by the seismic event location algorithm, these spatial coordinates are mapped to the established geological 3D model; and the microseismic events are associated with the spatial location of the target fracture. The association methods include:

[0135] The spatial coordinates of each microseismic event are obtained, an influence radius RD is set based on the spatial coordinates, and a circular area is constructed based on the influence radius RD to obtain the associated influence area. The target fracture in the associated influence area is associated with the microseismic event. Specifically:

[0136] Each microseismic event is marked to form an event sequence set {1, 2, …, j, …, M}, where j represents the index of the j-th microseismic event and M represents the total number of microseismic events; and the target fracture in the associated impact area corresponding to the j-th microseismic event is represented as ; where i1∈{1,2,…,n};

[0137] Define target cracks based on the Richter magnitude formula or moment magnitude formula The magnitudes of the corresponding j-th microseismic events are or ;

[0138] Earthquakes with magnitudes of 2 to 6 are considered small to moderate: the Richter scale is used.

[0139] Earthquakes with a magnitude greater than 6 are considered large earthquakes: the moment magnitude is used because it more accurately reflects the energy and size of the earthquake.

[0140] Using the Richter magnitude formula:

[0141] If the target rift The associated jth microseismic event is a small to medium-sized earthquake, and the amplitude has been measured near the epicenter. The magnitude of the event is calculated using the Richter magnitude formula:

[0142] ;

[0143] in, is the maximum amplitude of the jth microseismic event, is the distance from the epicenter of the microseismic event j to the observation point, is an empirical constant, which is taken as 1.0 in this embodiment;

[0144] Use the moment magnitude formula:

[0145] For large earthquakes or when studying deep rift faults, the seismic moment is considered, so the moment magnitude formula is used to calculate the magnitude of the jth microseismic event:

[0146] ;

[0147] in, is the seismic moment of the jth microseismic event, is the fault area, is the sliding displacement, is the stiffness of the rock, is the reference seismic moment constant.

[0148] The influence radius RD is determined as follows:

[0149] The influence radius RD is determined based on the Inversion-Based Radius Adjustment Algorithm, as follows:

[0150] The algorithm utilizes microseismic monitoring data and the spatial distribution of known target fractures, combined with inversion techniques, to dynamically optimize the influence radius (RD) of microseismic events to more accurately correlate microseismic events with target fractures. The core of the inversion process is to minimize the spatial error between the associated influence area of ​​the microseismic event and the fracture, thereby ensuring that the influence radius more closely matches the actual underground environment.

[0151] The main goal of the inversion algorithm is to iteratively optimize the influence radius based on the spatial coordinates of different microseismic events, focal mechanisms, and the spatial distribution of target fractures, thereby obtaining the optimal influence range that matches the actual geological conditions. This process includes setting the initial influence radius, applying the inversion algorithm, optimizing the objective function, and finally verifying and correcting the final results.

[0152] The algorithm steps are:

[0153] 1. Initial influence radius setting:

[0154] A preliminary impact radius is set based on the magnitude, focal depth, and geological conditions of the microseismic event. The impact radius is preliminarily estimated based on the magnitude or roughly set using geological attributes in the 3D geological model.

[0155] 2. Establish the objective function:

[0156] An objective function is defined to measure the matching degree between the associated impact area of ​​the microseismic event and the target fracture. The basic form of the objective function is to optimize the value of the impact radius by calculating the spatial distance error between the microseismic event and the fracture so that the associated impact area covers the target fracture.

[0157] 3. Select the inversion algorithm:

[0158] The inversion algorithm includes:

[0159] Least squares method: Applicable to simple optimization problems, it can directly obtain the optimal solution of the influence radius by minimizing the objective function;

[0160] Genetic algorithm: Applicable to nonlinear and multi-peak optimization problems, it simulates the process of biological evolution to find the global optimal solution;

[0161] Simulated annealing algorithm: Applicable to complex global optimization problems, it avoids local optimal solutions by simulating the physical annealing process and finds the global optimal influence radius;

[0162] 4. Perform inversion adjustment:

[0163] Taking the initial impact radius as the starting point, the objective function is optimized through the inversion algorithm to obtain the optimal impact radius. The inversion process will continuously adjust the value of the impact radius and evaluate the matching degree between the microseismic event impact area and the fracture.

[0164] 5. Verification and correction:

[0165] After the inversion adjustment is completed, the optimized influence radius is applied to the geological 3D model to verify the spatial correlation between the impact area of ​​the microseismic event and the target fracture. If the match is not ideal, the inversion adjustment can be performed again to further optimize the influence radius.

[0166] 6. Feedback Mechanism and Iterative Optimization

[0167] After each new microseismic event data and fracture information is added, the inversion adjustment process can be re-executed to gradually optimize the influence radius; as more microseismic data and fracture information are accumulated, the adjustment of the influence radius will become more and more accurate.

[0168] Further explanation: A mathematical model of crack propagation is constructed to obtain the predicted crack length of each target crack, including:

[0169] Fracture properties: Determine the initial fracture size from the fracture characteristics of each target fracture, and determine the physical properties of each target fracture; the initial fracture size includes length, width, and inclination; the physical properties include elastic modulus and fracture toughness;

[0170] External conditions: The temperature field, stress field, and fluid pressure of the target geothermal area are used as external conditions and as driving factors of the mathematical model. Based on the stress field and fluid pressure, the stress state and fluid pressure state of each target fracture are determined.

[0171] The stress state is defined as the distribution of internal and external forces acting on a target fracture under specific conditions. It is expressed as a tensor consisting of three principal stresses: maximum, intermediate, and minimum. In three-dimensional space, these stresses help describe the magnitude and direction of the forces acting on the target fracture in different directions.

[0172] Effect: Stress state affects the opening, closing and sliding of cracks.

[0173] The fluid pressure state is defined as: The fluid pressure state refers to the pressure state exerted by the fluid in the target fracture, which affects the stability and water conductivity of the fracture.

[0174] Function: Fluid pressure affects the water filling capacity, fluidity and response behavior of the fracture under stress field.

[0175] Fusion phase characteristics:

[0176] The extracted seismic phase characteristics are combined with the stress state and fluid pressure state of the target fracture as input parameters of the fracture propagation mathematical model to describe the initial and boundary conditions of fracture propagation.

[0177] Data weighting: Through parameter sensitivity analysis and statistical modeling, the impact of seismic phase characteristics, stress state, and fluid pressure state on fracture propagation behavior is determined, and appropriate weights are assigned to ensure the balance and accuracy of the numerical model;

[0178] Model solving and optimization:

[0179] Numerical simulation: The numerical simulation of crack propagation is performed using finite element software ABAQUS or COMSOL Multiphysics. Based on the mathematical model of crack propagation, the propagation direction and activity intensity of several target cracks in the target geothermal area are analyzed to obtain the predicted crack of each target crack. The predicted crack is the crack characteristics including crack size, crack density, and crack distribution law. The predicted crack of the target crack in the associated impact area corresponding to the jth microseismic event is expressed as The mathematical models of crack propagation in this embodiment include but are not limited to: extended finite element method (XFEM) and nonlinear fracture mechanics (EPFM).

[0180] EPFM combines elastic and plastic behaviors and is suitable for handling the plastic zone during fracture propagation. This is particularly important for fracture propagation in deep carbonate rocks under high temperature and high pressure environments, because under these conditions, rocks typically undergo significant plastic deformation.

[0181] XFEM is an efficient numerical simulation method that can simulate fracture propagation without re-meshing, making it particularly suitable for complex fracture propagation paths. For fracture propagation in deep carbonate reservoirs, XFEM can handle the coupled stress-fluid-heat-deformation problem with great flexibility.

[0182] In terms of parameter optimization: Genetic algorithm or particle swarm optimization algorithm is used to optimize the parameters of the mathematical model of crack expansion to minimize the error between the simulation results and the actual observation data;

[0183] Model Validation:

[0184] Historical data comparison: Compare the model prediction results with historical microseismic event data to verify the accuracy of the model;

[0185] Field verification: Through drilling or other geological survey methods, the fracture propagation predicted by the model is verified on site to further calibrate the model;

[0186] The steps for analyzing the crack extension direction and activity intensity are as follows:

[0187] 1. Analysis of expansion direction:

[0188] Principal stress direction: Determine the principal propagation direction of the target fracture by combining the geological 3D model with stress field analysis;

[0189] Seismic phase feature correlation: Using seismic phase propagation time and amplitude changes, we can identify the potential propagation direction of the crack and adjust the model parameters to reflect the actual situation.

[0190] 2. Activity intensity assessment:

[0191] Amplitude and frequency correlation: By analyzing the amplitude and frequency of P and S waves, the intensity of fissure activity and energy release can be assessed;

[0192] Dynamic monitoring: Real-time monitoring of microseismic events and dynamic update of activity intensity assessment to ensure timely reflection of dynamic changes in fractures;

[0193] 3. Identification of active and stable areas:

[0194] Active areas: Based on activity intensity assessment, identify high activity areas of target fractures as key areas for geothermal development;

[0195] Stable Zones: Identify low-activity areas of the fracture to serve as areas for auxiliary development and safety monitoring;

[0196] The division method of active area and stable area is as follows:

[0197] The intensity of activity is expressed by indicators consisting of displacement rate, stress-deformation response, fluid pressure change, and seismic phase activity / microseismic signals;

[0198] Displacement rate: The opening of the crack surface or the size of the shear displacement per unit time, used to measure the dynamics of crack expansion.

[0199] Stress-deformation response: the stress change rate or energy release rate in the area corresponding to the crack.

[0200] Stress concentration and rapid release often indicate highly active cracks.

[0201] Fluid pressure changes: The magnitude of the fluid pressure gradient within the fracture, specifically the diffusion rate of high-pressure fluid, can reflect the conductivity of the fracture and the strength of fluid coupling.

[0202] Seismic phase activity / microseismic signals: Record tiny signal activities through seismic monitoring and analyze the dynamic change intensity and frequency of cracks.

[0203] Data statistics and analysis: By collecting long-term monitoring data, applying statistical analysis and numerical simulation, we determine the specific thresholds for each indicator and use them to divide active and stable areas.

[0204] Model simulation and verification: Numerical simulation tools are used to conduct sensitivity analysis on the numerical model to verify the rationality of the threshold and accurately divide different sections.

[0205] Historical data comparison: Combined with historical activity data, a comprehensive assessment of the stability and changing trends of fracture activity in different regions is conducted to determine regional boundaries.

[0206] It is further explained that the multi-field coupling model of the target geothermal area is determined using Huang Yonghui's "Multi-field Coupling Simulation Technology and Its Application in Geothermal Exploitation Safety Assessment" or "A Multi-field Coupling Numerical Simulation Method for Geothermal Reservoirs" with publication number CN115310329A;

[0207] The steps for implementing the multi-field coupling model are as follows:

[0208] 1. Data preparation and preliminary analysis:

[0209] Collect basic data on geology, geophysics, and fluid properties of the target geothermal area.

[0210] Conduct preliminary analysis of the carbonate rock layers and reservoir characteristics in the region to determine the physical and chemical parameters required for simulation. Reservoir characteristics include but are not limited to porosity, fracture density, and fracture size;

[0211] Determine model parameters related to thermal, fluid, mechanical, and chemical reactions of geothermal reservoirs, such as thermal conductivity, porosity, permeability, reaction rates, etc.

[0212] 2. Define the coupling equations of the model:

[0213] Heat conduction equation: The heat conduction equation is established based on the thermal conductivity of the rock, fluid temperature, etc. The heat conduction equation is determined based on Fourier's law;

[0214] Fluid flow equation: The flow model of Darcy's law is used to describe the flow process of geothermal fluid.

[0215] Solid deformation equation: Through the stress-strain relationship, a mechanical model of solid deformation is established, taking into account the plastic, elastic or fracture behavior of rocks during geothermal mining.

[0216] Chemical reaction equation: Based on the reaction process of geothermal fluid and carbonate rock, a chemical reaction kinetic model is established, taking into account the reaction rate and product generation, etc.

[0217] 3. Determine boundary conditions and initial conditions:

[0218] According to the geological structure and mining history of the region, the boundary conditions of the simulation area are set, including heat flow, temperature, pressure, and fluid injection.

[0219] Determine the initial conditions, including the initial temperature field, pressure field, chemical composition, etc., to simulate the state of the reservoir at the beginning.

[0220] 4. Constructing numerical grid and discretization:

[0221] Grid the simulation area and select the appropriate grid division strategy based on the reservoir geometry, fracture distribution, physical properties and other information.

[0222] Discretize heat conduction, fluid flow, solid deformation, and chemical reaction equations using the finite element method, finite difference method, or finite volume method.

[0223] 5. Multi-field coupling calculation and solution:

[0224] The numerical solution is performed by coupling the four fields of heat, flow, force, and chemical reaction. The time evolution is performed using explicit or implicit time stepping algorithms.

[0225] Use iterative solvers to perform multi-field interactions to account for the interactions between heat flow, fluid flow, and solid deformation.

[0226] 6. Result analysis and verification:

[0227] Extract simulation results and analyze the distribution of various physical fields such as thermal field, flow field, stress field and chemical reaction field.

[0228] Verify the accuracy and reliability of the model. Compare it with field experimental data, data from existing literature, or theoretical models to confirm the rationality of the simulation results.

[0229] The sensitivity of the model is analyzed to study the impact of different parameter changes on the simulation results.

[0230] 7. Application and evaluation of results:

[0231] Based on the results of the multi-field coupling model, the thermal flow field data of each associated influence area in the target geothermal area are evaluated;

[0232] The heat flow field data include heat flux density, thermal conductivity, heat storage, and heat recovery rate; the heat flux density, thermal conductivity, heat storage, and heat recovery rate of the associated impact area corresponding to the jth microseismic event are defined as 、 、 、 .

[0233] Heat flux density is the heat flow per unit area, measured in watts per square meter (W / m²), and represents the intensity of heat transfer. It is calculated from thermal conductivity and temperature gradient:

[0234] ;

[0235] in, is the heat flux density, is the thermal conductivity, is the temperature gradient.

[0236] Thermal conductivity is the change in thermal conductivity with changes in temperature, pressure, and rock type; the calculation formula is as follows:

[0237] ;

[0238] is the thermal conductivity (unit: W / m·K); is the temperature (unit: K); is the pressure (unit: Pa); is the porosity;

[0239] Thermal Reservoir Capacity is expressed as Heat Capacity, which describes the amount of thermal energy that a reservoir can store. The calculation formula is:

[0240] ;

[0241] is the heat reserve of the reservoir (unit: J, joule); is the density of the rock (unit: kg / m³); is the specific heat capacity of the rock (unit: J / kg·K); is the volume of the reservoir (unit: m³); is the temperature change; this formula calculates the temperature change of the rock on the target fracture The stored heat under the given conditions takes into account the mass, specific heat capacity and volume of the rock.

[0242] Heat Recovery Rate: This parameter is used to evaluate the heat recovery rate under different mining conditions by simulating the mining process. This parameter determines the economic and sustainability of geothermal mining. The calculation formula is:

[0243] ;

[0244] is the heat recovery rate (unit: W, watt); It is the mass flow rate of the fluid (unit: kg / s); it represents the mass of the fluid flowing through the geothermal reservoir per unit time and is obtained by measuring the flow rate.

[0245] is the specific heat capacity of the fluid (unit: J / kg·K), which represents the amount of heat required to change the temperature of a unit mass of fluid; is the temperature change of the fluid (unit: K).

[0246] Further explanation: the predicted fractures of each target fracture are imported into the multi-field coupling model for simulation, and the heat flow field change data of each associated impact area are obtained. The change values ​​of heat flux density, thermal conductivity, heat storage and heat recovery rate of the associated impact area corresponding to the jth microseismic event are expressed as 、 、 、 ;

[0247] The heat flux density change value is defined as the change in heat flow transport intensity due to predicted fracture induction and geothermal mining activities.

[0248] Factors of change:

[0249] Changes in fractures are predicted to affect the path and intensity of heat flow;

[0250] The propagation of cracks is predicted to increase or decrease the local heat flux;

[0251] The calculation formula is ,in and They represent the new heat flux density after the crack is introduced and the heat flux density in the initial state.

[0252] The thermal conductivity change value is defined to reflect the change in heat transfer capacity within the target geothermal area due to the impact of the predicted fractures.

[0253] Factors of change:

[0254] Cracks are predicted to increase overall thermal conductivity by providing more heat conduction paths.

[0255] Conversely, cracks are predicted to reduce local thermal conductivity due to increased thermal resistance;

[0256] The calculation formula is: ,in and It is the thermal conductivity after the crack is introduced and in the initial state.

[0257] Definition: Thermal reserve change is the change in storable thermal energy due to extraction and changes in reservoir conditions.

[0258] Factors of change:

[0259] Predicted fractures will change the volume and structure of the reservoir, affecting the total heat storage.

[0260] The predicted increase in fractures improves the thermal conductivity of the reservoir, potentially increasing heat storage.

[0261] The calculation formula is ,in and It is the heat storage after the fracture is introduced and in the initial state.

[0262] Definition: Changes in thermal recovery rate are adjustments to heat recovery efficiency due to fracture introduction and mining conditions. Thermal recovery rate refers to the ratio of heat energy recovered from a geothermal reservoir to the heat energy initially injected or stored during geothermal energy development.

[0263] Factors of change:

[0264] It is predicted that cracks will change the flow path of the fluid and affect the effective recovery of thermal energy.

[0265] Changes in heat exchange efficiency can result in an increase or decrease in recovery rate.

[0266] The calculation formula is ,in and It is to predict the heat recovery rate after fracture introduction and in the initial state.

[0267] To further illustrate, the comprehensive contribution of the associated impact area corresponding to the target crack i1 of the jth microseismic event to the change of the thermal flow field is marked as , and As the influence coefficient, The specific calculation formula is as follows:

[0268] ;

[0269] in, is the target fracture number of the jth microseismic event Instead of predicting the crack number The absolute difference is used to represent the offset of the crack number; The larger the value, the greater the deviation between the predicted crack and the actual situation;

[0270] The meaning of the crack number offset is: It is the difference between the target fracture numbers, indicating the difference between the actual fracture number and the predicted fracture number when the microseismic event occurs. This offset reflects the error or deviation between the actual fracture and the predicted fracture.

[0271] The larger the offset, the greater the difference between the actual target fracture number and the predicted fracture number, which means that the range of fracture extension is more complex, or the deviation between the predicted model and the actual situation is greater.

[0272] The smaller the offset, the smaller the difference between the actual crack and the predicted crack, the simpler the crack expansion, and the more accurate the prediction result.

[0273] is the weight factor, which corresponds to the influence weight of heat flux, thermal conductivity, heat storage and heat recovery rate respectively. The value range is (0,1) and satisfies ;

[0274] The contribution of different influencing factors in the target geothermal area experiment is statistically fitted and the weight ratio is adjusted. Specifically, the weight determination method combining entropy weight method and fuzzy analytic hierarchy process (FAHP) is adopted to target The weight factor is used to determine the proportion;

[0275] In order to Conduct a detailed and limited analysis and provide quantitative parameters and descriptions based on different scenarios in actual applications. Analyze two typical scenarios:

[0276] Scenario 1 is a scenario with a large impact; Scenario 2 is a scenario with a small impact;

[0277] Large impact scenario: When any parameter among heat flux, thermal conductivity, heat storage or heat recovery rate has a great influence on the geothermal reservoir, the corresponding weight coefficient will occupy a large proportion.

[0278] Small impact scenario: When the influence of any parameter among heat flux, thermal conductivity, heat storage or heat recovery rate is small, the weight coefficient will be correspondingly small.

[0279] Description of scenarios with large impact:

[0280] If, during the development of geothermal reservoirs, the heat flux density has a significant impact on the heat flow field, while the heat recovery rate has a small impact; the weight coefficient is set to ;

[0281] It means that the influence of heat flux is significant and the influence ratio is large;

[0282] It means that thermal conductivity has the second greatest influence on the thermal flow field;

[0283] This indicates that heat storage has a smaller impact, but is still an important factor;

[0284] It indicates that the influence of heat recovery rate is the smallest;

[0285] In this scenario, the change of heat flux density has a stronger dominant effect on the change of heat flow field, so it is given a higher weight coefficient. , the weights of other factors are relatively low.

[0286] Description of scenarios with small impact:

[0287] If the heat flux changes slightly and the change in heat recovery rate plays an important role in reservoir stability, the weight coefficient is ;

[0288] Indicates that the heat flux density has little influence;

[0289] This indicates that the influence of thermal conductivity is also small;

[0290] It indicates that heat reserve is more important to reservoir temperature control and is given a higher weight;

[0291] It indicates that the heat recovery rate has a great influence on the effect of geothermal development;

[0292] In this scenario, the influence of heat reserves and heat recovery rate on the heat flow field is more prominent, so the corresponding weight coefficients are higher.

[0293] The weight coefficients for high-impact scenarios are determined based on the critical role of actual heat flux changes on reservoir temperature and energy recovery. In this implementation, through simulation experiments or historical data, it was determined that changes in heat flux have the most significant effect on reservoir thermal management, and therefore a higher weight was assigned.

[0294] The weight coefficients for scenarios with low impact are set based on the needs of technical optimization. For example, in the early stages of reservoir development, heat flux is relatively stable and the impact of heat recovery on the reservoir is more important. Therefore, in this scenario, the weight of heat flux is lower.

[0295] The optimization strategies for large-scale scenarios are as follows:

[0296] When the heat flux density has a greater impact on the thermal flow field, increase The value of can more effectively evaluate the changes in the thermal flow field and ensure its dominant position in geothermal development.

[0297] Accordingly, reduce , which means that when the heat recovery rate is relatively unimportant, there is no need to rely too much on the effect of the heat recovery rate on the change of the thermal flow field, thereby reducing the computational complexity.

[0298] The optimization strategies for scenarios with small impact are as follows:

[0299] When the heat flux changes slightly, the and The weights can fully reflect the impact of heat reserves and heat recovery rate, improve the accuracy of the assessment results, and are particularly important when optimizing geothermal recovery and reserve management.

[0300] Reduce and The weight of , means reducing the dependence on heat flux and thermal conductivity, so that the evaluation results focus more on the changes in heat reserves and recovery rate.

[0301] is a nonlinear adjustment factor, corresponding to the nonlinear weight of each change, with a value range of [0.5, 2]. The numerical determination rules are as follows:

[0302] When the nonlinear adjustment factor is greater than 1: the relative change will be amplified, indicating that the magnitude of the change in the physical quantity is more important to the overall influence coefficient.

[0303] When the nonlinear adjustment factor is in the interval (0, 1): the relative change will be compressed, indicating that the weight of the change in the physical quantity in the overall influence coefficient is small;

[0304] Example of technical significance: By adjusting the value of the nonlinear adjustment factor, reasonable weights are given to the relative changes of different physical quantities, reflecting the sensitivity of different physical quantities to changes in the thermal flow field during the crack expansion process.

[0305] For example, if the heat recovery rate is a key indicator for evaluating the impact of fractures on thermal flow field, then >1; the effect of thermal conductivity is not as significant as that of heat storage, so it is set <1;

[0306] The determination method is as follows:

[0307] According to historical experimental data or actual monitoring data, the change range of heat flux density, thermal conductivity, heat storage, heat recovery rate and the corresponding heat flow field change results are calculated respectively, and the optimal least square method of nonlinear fitting is used to find value.

[0308] Conduct multiple fracture extension simulation experiments and record the changing trends of heat flux density, thermal conductivity, heat storage, heat recovery rate and corresponding heat flow field influence coefficient in the associated impact area corresponding to the jth microseismic event; use the fitting nonlinear formula Compare the actual impact coefficient results and adjust Select the value that makes the formula fit the best.

[0309] is a normalization constant used to standardize the range of the influence coefficient so that The output value range is limited to (0,1);

[0310] is the distance decay rate of the target crack number difference, which is a dimensionless constant and has a value range of [1,10]; The determination method is as follows:

[0311] Experimental data analysis:

[0312] In laboratory or field observations, the variation in heat flow between different fractures is measured. Through statistical analysis, the relationship between fracture spacing and heat flow variation is evaluated to determine the effective distance of the fracture's influence range.

[0313] By fitting the experimental data, a mathematical model that can describe the variation of heat flow with distance is selected, and appropriate attenuation parameters are extracted from it.

[0314] Field experience and expert opinions: Combined with practical experience in geothermal development and mining engineering, communicated with experts in the field, and determined the reasonable value.

[0315] It means that the effect caused by the offset between the target crack and the predicted crack gradually decreases. The larger the distance, The smaller;

[0316] The calculation logic analysis between the parameters in the formula is as follows:

[0317] 1.1) Relationship between the calculation logic in the formula and the evaluation results of thermal flow field changes:

[0318] The following physical quantities are defined, including heat flux, thermal conductivity, heat storage, and heat recovery rate;

[0319] The relative ratio of each physical quantity change By corresponding weight coefficient and nonlinear factors By combining the above, we can obtain the comprehensive contribution of each physical quantity to the change of the thermal flow field. This part reflects the degree of influence of crack expansion on the thermal flow field, and by adjusting the weights and nonlinear factors, we can flexibly simulate the heat flow changes under different crack expansion situations.

[0320] When the change amplitude of a physical quantity is large, such as When it increases, the corresponding relative change value Increase, and thus the influence coefficient will also increase; by weighting the relative changes of various physical quantities and adjusting them through nonlinear factors, the influence of each target crack on the thermal flow field can be accurately evaluated;

[0321] Crack number offset It will decay according to the change of the crack position; the larger the crack number offset, the smaller the value of the exponential decay term, resulting in a smaller influence coefficient. Conversely, the closer the crack positions are, the larger the influence coefficient.

[0322] When the crack position changes, and As the distance between them increases, the effect of the cracks on the thermal flow field is significantly weakened, and this change is effectively reflected by the attenuation function. This design ensures that the changes in the thermal flow field will not be affected by cracks with large position changes, making the model more stable.

[0323] when When a certain factor increases, the influence of the relevant physical quantity is strengthened, and the degree of change in the thermal flow field also increases; this is because the contribution of the physical quantity to the change in the thermal flow field becomes greater, thereby increasing the influence coefficient, which helps to accurately predict the change in thermal flow caused by crack expansion;

[0324] Increase the nonlinear factor:

[0325] Increase nonlinear factors The value of can amplify the influence of the corresponding physical quantity change, making the formula more sensitive; by reasonably adjusting the value of the nonlinear factor, the thermal flow field changes under different crack expansion situations can be more finely adjusted;

[0326] The increase in the crack number offset leads to a decrease in the attenuation factor and ultimately a decrease in the impact coefficient. This design ensures that when the crack offset is too large, the model's prediction of the thermal flow field changes tends to be stable, thus avoiding overly drastic prediction results.

[0327] The formula is designed to accurately evaluate the changes in the thermal flow field under different crack expansion scenarios. Therefore, its applicable scenarios are mainly based on the following conditions:

[0328] Fracture propagation in geothermal reservoirs: The formula provides a quantitative assessment of the changes in the heat flow field during fracture propagation, which is applicable in geothermal production and reservoir management, especially when predicting the impact of fracture propagation on heat flow;

[0329] Microseismic events: The formula is used to analyze the impact of crack expansion on the thermal flow field based on the monitoring data of microseismic events. It is applicable to areas with active seismic activity and helps to determine the thermal flow impact of crack expansion.

[0330] like The increase in the value of any factor in the equation indicates that the corresponding heat flux, thermal conductivity, heat storage, and heat recovery rate are more important, and the influence coefficient The contribution of

[0331] Will The value range is limited to the interval (0,1), and the evaluation results of the thermal flow field change are as follows:

[0332] The closer it is to 0, the closer it is to the target crack. The smaller the impact on the thermal flow field changes;

[0333] The closer it is to 1, the closer it is to the target rift. The greater the impact on the thermal flow field changes.

[0334] When all parameters In the interval (0,0.2), When it approaches 0, it means that the thermal flow field changes little and the reservoir is in a stable state;

[0335] when In the interval (0.9, 1), it indicates that the thermal flow field changes greatly and the reservoir is in an unstable state.

[0336] Further explanation: According to the geological risk and thermal flow field change status of each target fracture, corresponding management and adjustment strategies are formulated, including:

[0337] For each target rift The impact coefficient value is combined with the magnitude database to analyze geological risks:

[0338] like and , then the target crack is determined The corresponding area is of high geological risk;

[0339] like and , it is judged as low geological risk;

[0340] in is the risk classification threshold;

[0341] The risk classification thresholds are determined through historical data analysis, experience accumulation, simulation experiments or model calculations; these thresholds play a key role in risk management, defining the degree of change or event that should be considered a "risk range" and providing a basis for further management measures and adjustment strategies.

[0342] The specific rules for risk classification thresholds are as follows:

[0343] is the threshold value of the thermal flow field change; it indicates the degree of influence of the target crack on the thermal flow field, and and The calculation rules are as follows:

[0344] The thermal flow field change state is expressed by the heat recovery rate, and the following calculation formula is extracted from the thermal flow field change data:

[0345] The calculation formula is ,in and It is the heat recovery rate after fracture introduction and in the initial state;

[0346] Classification of thermal flow field effects:

[0347] According to different The value range divides the thermal flow field change state into different levels:

[0348] Low impact range: and ;

[0349] Moderate impact range: and ;

[0350] High impact range: and ;

[0351] The value is set by simulation results or historical monitoring data of different crack effects. The common value range is:

[0352] If the impact of the target crack on the thermal flow field is less than 10%, it is considered a low-risk crack area;

[0353] If the impact of the target crack on the thermal flow field is between 10% and 30%, it is considered a medium-risk crack area.

[0354] If the impact of the target crack on the thermal flow field exceeds 30%, it is considered a high-risk crack area.

[0355] In this embodiment, Set to 0.7;

[0356] It is the threshold of magnitude change, and the magnitude is related to the expansion of the crack and geological stability;

[0357] The magnitude classifications are as follows:

[0358] Low magnitude range: ;

[0359] Moderate magnitude range: ;

[0360] High magnitude range: ;

[0361] The value is defined based on the magnitude of microseismic activity and the response characteristics of fracture expansion. Through model analysis of the relationship between magnitude changes and geological risks, the following thresholds are determined:

[0362] If the magnitude is less than 2.0, it means that the microseismic activity has little impact on the reservoir and is considered a low-risk area.

[0363] If the magnitude is between 2.0 and 6.0, it means that the expansion of the cracks may have a certain impact on the thermal flow field and is considered a medium-risk area.

[0364] If the magnitude exceeds 6.0, it means that the cracks will expand on a large scale, accompanied by significant changes in the thermal flow field. Emergency adjustment measures need to be taken and it is considered a high-risk area.

[0365] In this embodiment, Set to 0.3;

[0366] It is the safety magnitude threshold, which indicates the maximum acceptable magnitude of microseismic activity in different fracture areas. Exceeding this value will cause greater geological risks. Its determination rules are as follows:

[0367] It means that during a certain crack expansion process, when the microseismic magnitude exceeds this value, it will cause significant damage to the geothermal reservoir and affect the balance of the heat flow field.

[0368] The value of is set based on the seismic resistance of the geological region, historical magnitude data, and model predictions;

[0369] For example, for relatively stable carbonate formations, according to the systematic experience of the expert group, The value is set to a magnitude of 4.0. If the magnitude is greater than this value, it may cause reservoir structure damage or severe imbalance of thermal flow field.

[0370] For areas with complex geological conditions and easy expansion of cracks, The value can be appropriately lowered based on actual geological exploration results, such as setting it to magnitude 3.0.

[0371] Comprehensive risk classification rules:

[0372] Management and adjustment strategies for low-risk areas: Maintain existing injection and production plans, conduct continuous monitoring, and conduct regular assessments;

[0373] Management and adjustment strategies for medium-risk areas: adjust injection and production flow rates and increase the frequency of dynamic monitoring.

[0374] Management and adjustment strategies for high-risk areas: Immediately adjust the injection and production strategy, reduce injection and production pressure, and increase fracture sealing measures.

[0375] Based on the distribution of influence coefficients and the changing state of the thermal flow field, identify local high-risk crack areas and draw a thermal flow field risk assessment map;

[0376] Based on the assessment results of high-risk fracture areas, a thermal flow field optimization adjustment strategy is formulated, specifically:

[0377] Adjust the injection and production ratio of geothermal fluids in the target area;

[0378] Optimize injection pressure and temperature gradients in the target fracture direction;

[0379] Apply physical compensation measures to areas where thermal flow field changes significantly;

[0380] Continuously monitor the geological risks and thermal flow field changes of each target fracture and dynamically adjust strategies.

[0381] The following is a detailed description. The geological risk judgment rules are as follows:

[0382] according to and Jointly determine the geological risk level using the following rules:

[0383] ;

[0384] Interaction rules and parameter impact analysis:

[0385] Impact of changes on geological risks:

[0386] Rule: Increase or decrease The value of directly affects the geological risk level, but does not change the magnitude .

[0387] Specific description:

[0388] ;

[0389] Impact of changes on geological risks:

[0390] Rule: Increase or decrease magnitude , directly affects the risk level, but does not affect .

[0391] Specific description:

[0392] ;

[0393] Develop management and thermal flow field adjustment strategies:

[0394] Based on the geological risk assessment results and the changing state of the thermal flow field, quantitative management methods and adjustment strategies are adopted to optimize the distribution of the thermal flow field, improve the thermal energy utilization efficiency of the reservoir and reduce the potential geological disaster risks in high-risk fracture areas.

[0395] Regional priority sorting: Based on the geological risk level results, the priority of the reservoir area is determined and divided into high priority, medium priority and low priority:

[0396] High priority areas: fractures have an impact on the thermal flow field of >70% and a magnitude of ≥6.0;

[0397] Medium priority area: The impact of cracks on thermal flow field is 30%-70%, and the magnitude ;

[0398] Low priority areas: The impact of the fracture on the thermal flow field is <30% and the magnitude is <2.0.

[0399] Thermal flow field control strategy design, specific suggestions are as follows:

[0400] For the management and adjustment of high priority areas:

[0401] Strategy 1: Non-uniform injection and production regulation:

[0402] Increase the regional thermal fluid injection ratio by 50%, but reduce the recovery ratio to 60% of the current level, and reduce the temperature gradient of the injected fluid by 25%.

[0403] Technical effect: Enhance crack closure tendency and reduce crack expansion; reduce thermal flow field fluctuation by 30%-40%.

[0404] Strategy 2: Optimizing injection and production pressure by fracture orientation:

[0405] According to the target crack number Adjust the initial injection and production pressure to reduce it by 20%-30%, and control the pressure fluctuation range to ±5%.

[0406] Technical effect: Significantly reduces the damage of high-priority fractures to local reservoir stability, and reduces heat flow deviation by 15%-20%.

[0407] Management and adjustments for medium priority areas;

[0408] Strategy 1: Uniform injection and production flow regulation:

[0409] Adjust the injection and production liquid flow ratio in the target area to 1:1, and control the temperature gradient variation range to <15%.

[0410] Technical effect: Maintaining crack expansion within a controllable range and reducing the impact on the reservoir temperature field by 20%-30%.

[0411] Strategy 2: Dynamic crack monitoring and adjustment:

[0412] Regarding earthquake magnitude In the target fracture area, the dynamic monitoring frequency is increased by 30%-40% to reduce excessive injection and production.

[0413] Technical effect: Significantly reduces fracture shear deformation rate and improves overall reservoir stability by 15%-20%.

[0414] For the management and adjustment of low priority areas:

[0415] The strategy is to maintain the current injection-production strategy:

[0416] The injection and production pressure, flow rate, and temperature gradient in the current fracture area remain stable and do not require dynamic adjustment.

[0417] Technical effect: Ensures maximum reservoir efficiency while reducing management costs by 10%-15%.

[0418] The quantitative effects of the interaction rules and adjustment strategies are as follows:

[0419] Defining interaction rules is to adjust a parameter change in the strategy, and the technical effect on the thermal flow field or risk level must be clear; the following are examples:

[0420] ;

[0421] The beneficial effects of this embodiment are as follows:

[0422] Through the technical solution of this application, the development of deep carbonate geothermal reservoirs will be more precise and stable, realizing the continuous and effective development of geothermal resources while ensuring the sustainable use of the environment and resources; this solution not only provides a scientific basis and technical support for geothermal resource development, but also provides an important reference for the development of deep geothermal resources in other similar areas, and promotes the advancement and application of renewable energy utilization technology.

[0423] The above formulas are all dimensionless and numerical calculations. The formulas are obtained by collecting a large amount of data and performing software simulation to obtain the most recent real situation. The preset parameters in the formulas are set by technicians in this field according to actual conditions.

[0424] The above embodiments can be implemented in whole or in part by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented in whole or in part in the form of a computer program product. Those skilled in the art will appreciate that the units and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are performed by hardware or software depends on the specific application and design constraints of the technical solution.

[0425] The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units, and may be located in one place or distributed across multiple network units. Some or all of these units may be selected to achieve the purpose of this embodiment as needed.

[0426] The above is only a specific implementation method of the present application, but the scope of protection of the present application is not limited thereto. Any technician familiar with this technical field can easily think of changes or replacements within the technical scope disclosed in this application, which should be covered by the scope of protection of the present application.

Claims

1. A numerical simulation method for multi-field coupled heat transfer mechanism in deep carbonate geothermal reservoirs, characterized in that: The specific steps include: Step S1: obtaining a geological three-dimensional model of the target geothermal area and marking a number of target fractures in the geological three-dimensional model; Based on a number of target fractures on the three-dimensional geological model, a microseismic monitoring network is deployed in the target geothermal area to monitor in real time a number of microseismic events occurring in the target geothermal area during geothermal development, and to collect relevant preliminary microseismic event data, including but not limited to time estimation, spatial coordinate estimation point, magnitude estimation value, and focal mechanism estimation information of each microseismic event; Step S2: Processing the collected preliminary microseismic event data using an earthquake event location algorithm to accurately determine the time, spatial coordinates, focal mechanism, and magnitude of each microseismic event; Based on the spatial coordinates of each microseismic event, determining the associated impact area of ​​each microseismic event and determining the target fracture within the associated impact area; Step S3: extracting the seismic phase characteristics of each microseismic event, which include the propagation time, amplitude and frequency of the P wave and the S wave; Based on the seismic phase characteristics of each microseismic event, a mathematical model of crack propagation is constructed using the finite element method. This mathematical model is used to analyze the propagation direction and activity intensity of several target cracks in the target geothermal area, thereby obtaining the predicted crack size of each target crack. Step S4: obtaining a multi-field coupling model of the target geothermal region, wherein the multi-field coupling model is used to simulate heat conduction, fluid flow, solid deformation, and chemical reaction processes in the deep carbonate geothermal reservoir, and ultimately obtain thermal flow field data for each associated influence area in the target geothermal region; Step S5: obtaining the thermal flow field data of each associated influence area in the target geothermal area, and importing the predicted fracture of each target fracture in the associated influence area into the multi-field coupling model for separate simulation, obtaining the thermal flow field change data of each associated influence area, and analyzing the thermal flow field change data of each associated influence area to calculate the influence coefficient. The influence coefficient is used to evaluate the degree of change of the thermal flow field under different fracture extension situations, and finally obtaining the thermal flow field change evaluation result of each associated influence area; Step S6: Based on the calculated earthquake magnitude in each associated impact area and the predicted fracture of each target fracture, and through the thermal flow field change assessment results, corresponding management and adjustment strategies are formulated for the geological risk and thermal flow field change status of each target fracture to effectively control the heat flow of the geothermal reservoir.

2. The method for numerical simulation of multi-field coupled heat transfer mechanism in deep carbonate geothermal reservoirs according to claim 1, characterized in that: Mark several target fractures in the geological 3D model, including: Develop minimum requirements based on current heat transfer mechanisms and identify target geothermal areas: Construct a geological 3D model of the target geothermal area; Determine the three-dimensional location of each fracture in the target geothermal area, as well as the fracture characteristics of each fracture, including fracture size, fracture density, and fracture distribution pattern; The impact of each fracture characteristic on geothermal development is assessed, with the total score = size score + density score + distribution score; the larger the assessment value, the greater the corresponding impact; Based on these evaluation numerical results, these cracks are prioritized in descending order, and then divided into high-priority cracks, medium-priority cracks, and low-priority cracks according to the size of the total score value, and the cracks that meet the high-priority cracks are selected as target cracks; In the geological 3D model, each target fracture is assigned a number to form {1, 2, ..., i, ..., n}, where i represents the index of the target fracture and n represents the total number of target fractures; Deploy a microseismic monitoring network in the target geothermal area, specifically: The microseismic monitoring network consists of seismometers and seismic sensors; Seismometers and seismic sensors are placed in boreholes at different depths within the target geothermal area, and combined with surface and underground equipment to form a multi-level microseismic monitoring network; and each monitoring location in the microseismic monitoring network is marked as a monitoring point.

3. The method for numerical simulation of multi-field coupled heat transfer mechanism in deep carbonate geothermal reservoirs according to claim 2, characterized in that: The collected preliminary microseismic event data are processed using earthquake event location algorithms to accurately determine the time, spatial coordinates, focal mechanism, and magnitude of each microseismic event. Specifically, the following methods are used: The preliminary microseismic event data collected in the microseismic monitoring network also includes: timestamp, arrival time difference and seismic wave information, including: The timestamp is the monitoring time point when each microseismic event occurs; The arrival time difference represents the propagation time difference from the earthquake source to each monitoring point; Seismic wave information is composed of the arrival time and amplitude of P and S waves; Seismic event location algorithms include: using 3D seismic inversion technology to derive and correct preliminary microseismic event data for each microseismic event based on the microseismic monitoring network and 3D geological model, thereby accurately determining the time, spatial coordinates, focal mechanism, and magnitude of each microseismic event; Focal mechanisms include normal faults, reverse faults, and strike-slip faults. They also include non-tectonic focal mechanisms, which include volume expansion, vertical fracture movement, and fluid-induced. The magnitude of each microseismic event is calculated using a seismological formula based on the arrival time and amplitude information of the P and S waves corresponding to the seismic wave information; the seismological formula used is the Richter magnitude formula or the moment magnitude formula; After obtaining the spatial coordinates of each microseismic event determined by the seismic event location algorithm, these spatial coordinates are mapped to the established geological 3D model; and the microseismic events are associated with the spatial location of the target fracture. The association methods include: The spatial coordinates of each microseismic event are obtained, an influence radius RD is set based on the spatial coordinates, and a circular area is constructed based on the influence radius RD to obtain the associated influence area. The target fracture in the associated influence area is associated with the microseismic event. Specifically: Each microseismic event is marked to form an event sequence set {1, 2, …, j, …, M}, where j represents the index of the j-th microseismic event and M represents the total number of microseismic events; and the target fracture in the associated impact area corresponding to the j-th microseismic event is represented as ; where i1∈{1,2,…,n}; Based on the earthquake severity, choose the Richter magnitude formula or the moment magnitude formula to calculate the magnitude. When calculating the magnitude based on the Richter magnitude formula, define the target fracture. The magnitude of the corresponding j-th microseismic event is , when calculating the magnitude based on the moment magnitude formula, define the target crack The magnitude of the corresponding j-th microseismic event is .

4. The method for numerical simulation of multi-field coupled heat exchange mechanism in deep carbonate geothermal reservoirs according to claim 3, characterized in that: Construct a mathematical model of crack propagation to obtain the predicted crack length of each target crack, including: determining an initial fracture size from the fracture characteristics of each target fracture, and determining physical properties of each target fracture; The temperature field, stress field and fluid pressure of the target geothermal area are taken as external conditions and used as driving factors of the mathematical model; and based on the stress field and fluid pressure, the stress state and fluid pressure state of each target fracture are determined.

5. The method for numerical simulation of multi-field coupled heat transfer mechanism in deep carbonate geothermal reservoirs according to claim 4, characterized in that: The extracted seismic phase characteristics are combined with the stress state and fluid pressure state of the target fracture as input parameters of the fracture propagation mathematical model to describe the initial and boundary conditions of fracture propagation. Finite element software ABAQUS or COMSOL Multiphysics is used to perform numerical simulation of crack propagation. Based on the mathematical model of crack propagation, the propagation direction and activity intensity of several target cracks in the target geothermal area are analyzed to obtain the predicted crack of each target crack. The predicted crack is the crack characteristics including crack size, crack density, and crack distribution law. The predicted crack of the target crack in the associated impact area corresponding to the jth microseismic event is expressed as .

6. The method for numerical simulation of multi-field coupled heat transfer mechanism in deep carbonate geothermal reservoirs according to claim 5, characterized in that: Based on the results of the multi-field coupling model, the thermal flow field data of each associated influence area in the target geothermal area are evaluated; The heat flow field data include heat flux density, thermal conductivity, heat storage, and heat recovery rate; the heat flux density, thermal conductivity, heat storage, and heat recovery rate of the associated impact area corresponding to the jth microseismic event are defined as 、 、 、 .

7. The method for numerical simulation of multi-field coupled heat transfer mechanism in deep carbonate geothermal reservoirs according to claim 6, characterized in that: The predicted fractures of each target fracture are imported into the multi-field coupling model for simulation to obtain the heat flow field change data of each associated impact area. The change values ​​of heat flux density, thermal conductivity, heat storage and heat recovery rate of the associated impact area corresponding to the jth microseismic event are expressed as 、 、 、 ; The comprehensive contribution of the associated impact area corresponding to the target crack i1 of the jth microseismic event to the change of the thermal flow field is marked as , and As the influence coefficient, The specific calculation formula is as follows: in, is the target fracture number of the jth microseismic event Instead of predicting the crack number The absolute difference is used to represent the offset of the crack number; The larger the value, the greater the deviation between the predicted crack and the actual situation; is the weight factor, which corresponds to the influence weight of heat flux, thermal conductivity, heat storage and heat recovery rate respectively. The value range is (0,1) and satisfies ; is a nonlinear adjustment factor, corresponding to the nonlinear weight of each change, and its value range is [0.5, 2]; is a normalization constant used to standardize the range of the influence coefficient so that The output value range is limited to (0,1); is the distance decay rate of the target crack number difference, which is a dimensionless constant and has a value range of [1,10]; like The increase in the value of any factor in the equation indicates that the corresponding heat flux, thermal conductivity, heat storage, and heat recovery rate are more important, and the influence coefficient The contribution of Will The value range is limited to the interval (0,1), and the evaluation results of the thermal flow field change are as follows: The closer it is to 0, the closer it is to the target crack. The smaller the impact on the thermal flow field changes; The closer it is to 1, the closer it is to the target rift. The greater the impact on the thermal flow field changes.

8. The method for numerical simulation of multi-field coupled heat transfer mechanism in deep carbonate geothermal reservoirs according to claim 7, characterized in that: Develop corresponding management and adjustment strategies based on the geological risks and thermal flow field changes of each target fracture, including: For each target rift The impact coefficient value is combined with the magnitude database to analyze geological risks: like and , then the target crack is determined The corresponding area is of high geological risk; like and , it is judged as low geological risk; in is the risk classification threshold; Based on the distribution of influence coefficients and the changing state of the thermal flow field, identify local high-risk crack areas and draw a thermal flow field risk assessment map; Based on the assessment results of high-risk fracture areas, a thermal flow field optimization adjustment strategy is formulated, specifically: Adjust the injection and production ratio of geothermal fluids in the target area; Optimize injection pressure and temperature gradients in the target fracture direction; Apply physical compensation measures to areas where thermal flow field changes significantly; Continuously monitor the geological risks and thermal flow field changes of each target fracture and dynamically adjust strategies.

Images