Carbon dioxide geological storage modeling method, device and electronic equipment

Abstract

The application provides a carbon dioxide geological storage modeling method and device and electronic equipment, comprising: constructing a three-dimensional geological grid model of a target area, and giving physical property parameters to each unit of the three-dimensional geological grid model; wherein the three-dimensional geological grid model comprises a conventional fault unit and a hidden fault unit; obtaining experimental parameters of the hidden fault unit and the conventional fault unit determined through indoor experiments, and establishing a fault activation dynamic criterion based on the experimental parameters and actual shear stress of the fault; performing multi-physical field coupling numerical simulation based on the physical property parameters, and generating a fault activation risk evolution space-time cloud map based on the simulation results and the fault activation dynamic criterion. The application improves the accuracy of the activation risk prediction of the hidden fault.

Description

Technical Field

[0001] This invention relates to the field of geological modeling technology, and in particular to a method, apparatus and electronic equipment for modeling carbon dioxide geological sequestration. Background Technology

[0002] In existing geological modeling practices, faults are generally characterized as static geometric interfaces, and their mechanical properties are usually set as constant parameters. Typical approaches include: treating faults as rigid discontinuities to constrain the mesh topology, or assigning them fixed internal friction angles and cohesion values; some models introduce threshold-based activation criteria, such as setting a single critical pore pressure as the fault slip initiation condition.

[0003] In the scenario of carbon dioxide geological sequestration, the long-term injection process is accompanied by significant multiphysics field evolution phenomena. Field monitoring and laboratory core experiments show that as injection continues, carbon dioxide dissolves in formation water to form carbonic acid, leading to changes in the microstructure of fault fillers and a decrease in macroscopic shear strength. Simultaneously, the diffusion of the injection pressure field causes an increase in pore pressure in the surrounding rock, resulting in dynamic evolution of normal and shear stresses on the fault plane. However, current mainstream modeling frameworks do not incorporate such time-varying chemical corrosion effects and real-time stress field feedback mechanisms into their constitutive description systems, leading to a significant underestimation of the fault activation risk under long-term injection. Summary of the Invention

[0004] In view of this, the purpose of the present invention is to provide a method, apparatus and electronic device for modeling carbon dioxide geological storage, so as to improve the accuracy of predicting the activation risk of hidden faults.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows: In a first aspect, the present invention provides a method for modeling carbon dioxide geological storage, comprising: constructing a three-dimensional geological grid model of a target area and assigning physical property parameters to each unit of the three-dimensional geological grid model; wherein the three-dimensional geological grid model includes conventional fault units and hidden fault units; obtaining experimental parameters of hidden fault units and conventional fault units determined by indoor experiments, and establishing a dynamic criterion for fault activation based on the experimental parameters and the actual shear stress of the fault; performing multi-physics coupled numerical simulation based on the physical property parameters, and generating a spatiotemporal cloud map of fault activation risk evolution based on the simulation results and the dynamic criterion for fault activation.

[0006] Optionally, a three-dimensional geological grid model of the target area is constructed, including: acquiring basic data of the target area, and extracting micro-structural anomaly features from imaging well logging, micro-fracture surface features from core samples, and weak reflection phase axis misalignment features from seismic data based on the basic data; identifying hidden fault units based on the micro-structural anomaly features from imaging well logging, micro-fracture surface features from core samples, and weak reflection phase axis misalignment features from seismic data, combined with stratigraphic sedimentary facies analysis; and establishing a three-dimensional geological grid model containing stratigraphic units, reservoir units, conventional fault units, and hidden fault units using the finite element mesh generation method.

[0007] Optionally, indoor experiments include: under simulated formation temperature and pressure conditions, exposing fault rock samples to a carbon dioxide-containing brine environment for different time periods, and measuring their cohesion and friction coefficient; different time periods covering the time scales of the initial, middle and late injection stages; and using a double-shear friction test apparatus to measure the activation critical shear stress of conventional fault samples and latent fault samples under different normal stresses.

[0008] Optionally, the experimental parameters of the latent fault unit and the conventional fault unit determined by indoor experiments are obtained, including: fitting the measured cohesion and friction coefficient with an exponential decay model to obtain the experimental relationship that the cohesion and friction coefficient of the latent fault unit and the conventional fault unit decrease with the carbon dioxide exposure time; fitting the normal stress and the activation critical shear stress to establish the fitting relationship between the critical shear stress and the normal stress.

[0009] Optionally, multi-physics coupled numerical simulations can be performed based on physical property parameters, including: using an implicit iterative solution method to sequentially solve the seepage field, mechanical field, chemical field and temperature field of the target region to obtain the simulation results of each field.

[0010] Optionally, based on the simulation results and the fault activation dynamic criteria, a spatiotemporal cloud map of fault activation risk evolution is generated, including: updating the fault activation dynamic criteria based on the simulation results, and calculating the fault activation probability at different time points based on the updated fault activation dynamic criteria; determining the fault activation risk level based on the fault activation probability, and generating a spatiotemporal cloud map of fault activation risk evolution based on the fault activation risk level.

[0011] Optionally, after generating the spatiotemporal cloud map of fault activation risk evolution, the method also includes: based on the spatial distribution characteristics of high-risk areas in the spatiotemporal cloud map of fault activation risk evolution, outputting a monitoring deployment scheme and an injection pressure control strategy.

[0012] Secondly, the present invention provides a carbon dioxide geological storage modeling device, comprising: a model construction module for constructing a three-dimensional geological grid model of a target area and assigning physical property parameters to each unit of the three-dimensional geological grid model; wherein the three-dimensional geological grid model includes conventional fault units and hidden fault units; a dynamic criterion determination module for acquiring experimental parameters of hidden fault units and conventional fault units determined through indoor experiments, and establishing a fault activation dynamic criterion based on the experimental parameters and the actual shear stress of the fault; and a numerical simulation module for performing multi-physics field coupled numerical simulation based on the physical property parameters, and generating a spatiotemporal cloud map of fault activation risk evolution based on the simulation results and the fault activation dynamic criterion.

[0013] Thirdly, the present invention provides an electronic device including a processor and a memory, the memory storing computer-executable instructions executable by the processor, the processor executing the computer-executable instructions to implement the steps of the method provided in any of the first aspects above.

[0014] Fourthly, the present invention provides a computer-readable storage medium storing a computer program, wherein the computer program, when executed by a processor, performs the steps of the method provided in any of the first aspects above.

[0015] This invention brings the following beneficial effects: The carbon dioxide geological storage modeling method, apparatus, and electronic equipment provided by this invention first construct a three-dimensional geological grid model of the target area and assign physical property parameters to each unit of the three-dimensional geological grid model. The three-dimensional geological grid model includes conventional fault units and hidden fault units. Then, experimental parameters of the hidden fault units and conventional fault units, determined through indoor experiments, are obtained, and a dynamic criterion for fault activation is established based on the experimental parameters and the actual shear stress of the fault. Finally, multiphysics coupled numerical simulation is performed based on the physical property parameters, and a spatiotemporal cloud map of fault activation risk evolution is generated based on the simulation results and the fault activation dynamic criterion. In this method, by constructing a three-dimensional geological grid model containing conventional and hidden fault units, determining the experimental parameters of the hidden and conventional fault units through indoor experiments, establishing a dynamic criterion for fault activation based on the actual shear stress of the fault, and finally generating a spatiotemporal cloud map of fault activation risk evolution through the simulation results of multiphysics coupled numerical simulation and the fault activation dynamic criterion, the accuracy of predicting the activation risk of hidden faults is improved.

[0016] Other features and advantages of the invention will be set forth in the description which follows, and will be apparent in part from the description, or may be learned by practicing the invention. The objects and other advantages of the invention are realized and obtained in accordance with the structures particularly pointed out in the description, claims and drawings.

[0017] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, preferred embodiments are described below in detail with reference to the accompanying drawings. Attached Figure Description

[0018] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0019] Figure 1 A flowchart illustrating a carbon dioxide geological storage modeling method provided in this embodiment of the invention; Figure 2 This is a schematic diagram of a carbon dioxide geological storage modeling device provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] Currently, mainstream modeling frameworks do not incorporate the time-varying chemical corrosion effect and real-time stress field feedback mechanism during long-term injection into the constitutive description system, resulting in models severely underestimating the risk of fault activation under long-term injection.

[0022] Based on this, the present invention provides a carbon dioxide geological storage modeling method, apparatus and electronic equipment to improve the accuracy of predicting the activation risk of hidden faults.

[0023] To facilitate understanding of this embodiment, a detailed description of a carbon dioxide geological storage modeling method disclosed in this invention will be provided first. This method can be executed by electronic devices, such as smartphones, computers, and tablets. See also... Figure 1 The flowchart shown illustrates a method for modeling carbon dioxide geological storage, which mainly includes the following steps S101 to S103: Step S101: Construct a three-dimensional geological grid model of the target area and assign physical property parameters to each element of the three-dimensional geological grid model.

[0024] In one implementation, firstly, basic data of the target area is acquired, and based on the basic data, micro-structural anomaly features of imaging logging, micro-fracture surface features of rock cores, and weak reflection phase axis misalignment features of seismic data are extracted.

[0025] In practice, basic data such as seismic interpretation data (3D seismic profile, pre-stack depth migration data), well logging data (sonic, density, resistivity, imaging logging), drilling core data, formation pressure and temperature data of the target area are collected. Micro-structural anomaly features of imaging logging, micro-fracture surface features of core data, and weak reflection phase axis misalignment features of seismic data are extracted.

[0026] Specifically, for the whole-bore imaging logging data of the target area, the formation dip vector is extracted at 0.125-meter intervals. Weak anomaly vector clusters with dip values ​​between 2° and 8° are screened. The screened weak anomaly vectors are projected onto the Schmidt grid according to depth, and azimuth clustering statistics are performed. The azimuth clustering results are then overlaid with the relative content changes of uranium, thorium, and potassium elements in natural gamma ray spectroscopy logging to obtain the micro-structural anomaly characteristics of imaging logging.

[0027] Identification of microfracture features in core samples includes scanning each field of view under a polarizing microscope, focusing on identifying the following microscopic markers of latent faults: deformation bands, microfracture zones, and grain concentration zones.

[0028] Identifying the weak reflection phase axis fault characteristics of seismic data includes: slightly shifting the seismic trace along the target direction (1-3 sampling points) and calculating the rate of change of the similarity coefficient of the gather before and after the shift. When a hidden fault passes through, the rate of change of the similarity coefficient will show local extreme values, and these extreme value bands often extend linearly along a specific orientation; the coherence anomaly, differential similarity anomaly, and linear anomaly in the curvature attribute are displayed by RGB overlay, and the linear strips with high overlap among the three are the hidden fault trend zones at the seismic scale.

[0029] Secondly, based on the micro-structural anomaly characteristics of imaging logging, the micro-fracture surface characteristics of core samples, and the weak reflection phase axis misalignment characteristics of seismic data, combined with stratigraphic sedimentary facies analysis, hidden fault units are identified.

[0030] In practice, based on the micro-structural anomaly characteristics of imaging logging, the micro-fracture surface characteristics of core samples, and the weak reflection phase axis misalignment characteristics of seismic data, combined with stratigraphic sedimentary facies analysis, the geometric parameters such as the planar distribution, vertical extension, and dip angle of hidden faults are identified, thus distinguishing hidden faults from conventional faults.

[0031] Specifically, the fault trend zone is traced layer by layer on horizontal slices of the seismic coherence body (generally one slice every 20 ms), and the planar position of the anomalous zone on each slice is recorded. The positions of each slice are projected onto the same plane, and the envelope of the anomalous zone is delineated. The planar distribution of hidden faults usually exhibits distinct segmentation, with each segment arranged in a right- or left-order oblique sequence, rather than a single linear pattern. Seismic profiles of interconnected wells are cut along the direction perpendicular to the fault strike. The micro-anomaly well segments determined by imaging logging are calibrated with the weak in-phase axis disturbance zones in the profiles to determine the vertical development interval of the hidden fault. The true dip angle of the micro-fracture surface interpreted by imaging logging at each well point and the apparent dip angle converted from the weak in-phase axis disturbance on the seismic profile are comprehensively statistically analyzed, and the weighted average is taken as the final dip angle estimate.

[0032] Furthermore, the decoded geometric parameters of the hidden fault (segmented planar distribution, stratified vertical extension, and discrete dip distribution) were compared with the corresponding parameters of known conventional faults in the work area to create a comparative chart. Combined with the regional stress field background and the distribution of basement fractures, it was determined whether the hidden fault was related to the activity of the weak zone in the basement. If the strike of the hidden fault was consistent with the conjugate shear direction derived from the regional principal stress direction, and its distribution was controlled by basement fractures, then its hidden nature was further corroborated.

[0033] Finally, a three-dimensional geological grid model containing stratigraphic units, reservoir units, conventional fault units, and hidden fault units was established using the finite element mesh generation method.

[0034] In practical implementation, the finite element method is adopted to establish a three-dimensional geological grid model containing stratigraphic units, reservoir units, conventional fault units, and hidden fault units. The grid is refined in the fault and reservoir contact zone (grid accuracy ≤ 5m), and the grid in non-critical areas is appropriately coarsened (grid accuracy ≤ 50m). Basic physical properties are assigned to each unit. Among them, the fault units need to be distinguished between conventional and hidden types, and the corresponding infill material type (conventional faults: mainly breccia; hidden faults: mainly mudstone smears and microcrystalline mylonite), initial porosity, initial permeability, initial cohesion C0, initial friction coefficient μ0, and other parameters are assigned.

[0035] Step S102: Obtain the experimental parameters of the hidden fault unit and the conventional fault unit determined by indoor experiments, and establish a dynamic criterion for fault activation based on the experimental parameters and the actual shear stress of the fault.

[0036] In one implementation, firstly, experimental parameters such as cohesion and friction coefficient of the two types of fault units are measured through indoor rock mechanics tests to clarify the differences in their mechanical properties; secondly, based on the actual occurrence conditions of the fault, a stress monitoring system is set up to obtain real-time distribution data of normal stress and shear stress on the fault plane; then, the experimental parameters are coupled with the field stress to establish a fault activation discrimination criterion.

[0037] In this embodiment of the invention, the criterion is improved to a dynamic activation criterion that considers the time effect, taking into account experimental parameters, while also adapting to the parameter differences between conventional faults and hidden faults. The criterion formula is as follows:

[0038] in, t The actual shear stress of the fault (MPa) is calculated from the real-time stress field. C ( t The cohesive strength (MPa) decreases with carbon dioxide immersion time, calculated using the experimentally determined chemical weakening function. The values ​​for hidden faults and conventional faults are substituted into the corresponding values. α value; m ( t The coefficient of friction decreases with carbon dioxide immersion time, calculated according to the experimentally determined chemical weakening function. Substitute the values ​​for hidden faults and conventional faults into the corresponding... β value; s n The effective normal stress (MPa) on the fault plane; P The fault pore pressure (MPa) is calculated from the real-time seepage field.

[0039] When the actual shear stress of the fault t When the dynamic critical shear stress on the right is greater than or equal to the stress on the right, the fault is deemed to have a risk of activation. t The higher the ratio of the critical shear stress to the critical shear stress, the higher the activation probability.

[0040] Step S103: Perform multi-physics field coupled numerical simulation based on physical property parameters, and generate a spatiotemporal cloud map of fault activation risk evolution based on simulation results and fault activation dynamic criteria.

[0041] In one implementation, the control equations of the temperature field, seepage field, mechanical field, and chemical field are solved simultaneously in a finite element simulator (preferably COMSOL Multiphysics or FLAC3D, which supports multi-field coupled solution), realizing the full-coupled dynamic simulation of the four fields. The simulation results and fault activation dynamic criteria are combined to generate a spatiotemporal cloud map of fault activation risk evolution.

[0042] The carbon dioxide geological storage modeling method provided in this invention constructs a three-dimensional geological grid model containing conventional fault units and hidden fault units, and determines the experimental parameters of hidden fault units and conventional fault units through indoor experiments. It establishes a dynamic criterion for fault activation by combining the actual shear stress of the fault. Finally, it generates a spatiotemporal cloud map of fault activation risk evolution by using the simulation results of multi-physics coupled numerical simulation and the dynamic criterion for fault activation, thereby improving the accuracy of predicting the activation risk of hidden faults.

[0043] In this embodiment of the invention, after constructing a three-dimensional geological grid model of the target area, the method further includes: during the multiphysics coupled numerical simulation, calculating the equivalent plastic strain of each element based on the changes in pore pressure and stress of each element; when the equivalent plastic strain of any element exceeds a preset strain threshold, marking the element as a failed element, and adaptively re-dividing the grid within a preset range around the failed element to increase the grid density, and mapping the physical property parameters on the re-divided grid nodes from the original grid to the new grid through interpolation.

[0044] In one embodiment, fault core samples from the target area were selected, and standard rock samples (50 mm in diameter and 100 mm in height) of conventional faults (breccia) and concealed faults (mudstone smears, mylonite) were prepared. These samples were dried, vacuumed, and then saturated with simulated formation brine. The aforementioned laboratory experiments included: (1) Carbon dioxide-salt water-rock interaction experiment: Under simulated formation temperature and pressure conditions, fault rock samples were continuously exposed to a carbon dioxide-containing saline environment for different time periods, and their cohesion and friction coefficient were measured respectively; different time periods covered the time scales of the initial, middle and late injection stages.

[0045] In practice, the experimental setup includes: a high-temperature and high-pressure reactor (maximum pressure 50 MPa, maximum temperature 200 °C), a carbon dioxide partial pressure control system, a temperature control system, a rock shear strength testing module, and a data acquisition system.

[0046] Rock samples soaked in saturated brine were placed in a high-temperature and high-pressure reactor. Experimental conditions were set (simulating the formation temperature of the target area and different partial pressures of carbon dioxide: 5MPa / 10MPa / 15MPa). Rock samples were taken out after soaking for 0d / 7d / 30d / 90d / 180d, respectively, and direct shear tests were carried out to determine the cohesion C and friction coefficient μ of the rock samples.

[0047] (2) Using a double shear friction test apparatus, the activation critical shear stress of conventional fault samples and hidden fault samples was determined under different normal stresses.

[0048] In practice, a double-shear friction test apparatus was used to determine the activation critical shear stress of conventional and latent fault samples under different normal stresses (5MPa / 10MPa / 15MPa / 20MPa). t cr .

[0049] Based on this, experimental parameters for hidden fault units and conventional fault units, determined through indoor experiments, were obtained, including: (1) The measured cohesion and friction coefficient were fitted using an exponential decay model to obtain the experimental relationship between the cohesion and friction coefficient of the hidden fault unit and the conventional fault unit, which showed a decay trend with carbon dioxide exposure time.

[0050] In practice, based on different soaking times... C , m The measured values ​​were fitted to the chemical weakening function using an exponential decay model, yielding the following results:

[0051] in, C ( t ), m ( t ) respectively soaking t Cohesion and coefficient of friction after time C 0、 m 0 represents the initial cohesion and the initial coefficient of friction. α The cohesive chemical corrosion rate coefficient, β The coefficient of friction is the chemical corrosion rate coefficient. t Carbon dioxide immersion time (d); Hidden fault α , β The value is 1.5-2.0 times that of conventional faults (because the infill material of hidden faults is more easily corroded by carbon dioxide). For conventional faults: For hidden faults: .

[0052] The higher the partial pressure of CO2, the larger the α and β values, the faster the fault strength decays, and the higher the risk of activation. When the partial pressure is >20 MPa, the α and β values ​​tend to stabilize. As the formation temperature increases, the chemical activity of carbonic acid increases, the α and β values ​​increase, and the fault corrosion rate accelerates. Fine-grained argillaceous and microcrystalline infills (hidden faults) have larger α and β values ​​than coarse-grained breccia infills (conventional faults) and are more easily corroded. The higher the formation water salinity, the lower the CO2 solubility, the smaller the α and β values, and the slower the fault corrosion rate.

[0053] (2) Fit the normal stress and the activation critical shear stress to establish the fitting relationship between the critical shear stress and the normal stress. In specific implementation, the fitting relationship between the critical shear stress and the normal stress is established to provide a threshold reference for dynamic criteria.

[0054] In one implementation, for the aforementioned step S103, that is, when performing multiphysics coupling numerical simulation based on physical property parameters, the following methods can be used, including but not limited to: using an implicit iterative solution method to sequentially solve the seepage field, mechanical field, chemical field and temperature field of the target region to obtain the simulation results of each field.

[0055] In practical implementation, the governing equations for the temperature field, seepage field, mechanical field, and chemical field are solved simultaneously in the finite element simulator to achieve a fully coupled dynamic simulation of the four fields. The initial and boundary conditions for the simulation are set as follows: Initial conditions: original formation temperature distribution, original pore pressure distribution, original geostress field (maximum / minimum horizontal principal stress, vertical stress) in the target area, and initial carbon dioxide concentration of 0; Boundary conditions: The outer boundary of the model is a stress-fixed boundary (maintaining the original geostress), the top and bottom boundaries of the reservoir are water-proof boundaries, the carbon dioxide injection well is a constant flow / constant pressure injection boundary, and the temperature field is an adiabatic boundary (ignoring formation heat exchange).

[0056] Neglecting the heat of convection and chemical reaction during carbon dioxide transport, and considering only the initial temperature distribution of the formation, the temperature field is a steady-state field that satisfies the heat conduction equation:

[0057] The seepage field considers two-phase seepage of carbon dioxide and brine, and uses a modified Darcy's law formula to describe the change in pore pressure P distribution caused by carbon dioxide migration:

[0058] The mechanical field, based on the effective stress principle and considering the porosity elasticity effect, describes the stress field caused by pressure changes. s Updated to satisfy the elasticity equilibrium equations:

[0059] The chemical field considers the diffusion and dissolution of carbon dioxide, describing the control of carbon dioxide concentration distribution on the intensity attenuation of fault units, and satisfies the carbon dioxide concentration diffusion equation:

[0060] The four-field coupling solution employs an implicit coupling iterative method, with the solution order being: seepage field → mechanical field → chemical field → temperature field. Time steps are set to 1 day (initial injection), 30 days (mid-injection), and 90 days (late injection). The convergence criterion is that the iteration error of each field's calculation parameters ≤ 10. -6 .

[0061] In one implementation, for the aforementioned step S103, i.e., when generating a spatiotemporal cloud map of fault activation risk evolution based on simulation results and fault activation dynamic criteria, the following methods may be adopted, including but not limited to: First, the dynamic criterion for fault activation is updated based on the simulation results, and the probability of fault activation at different time points is calculated based on the updated dynamic criterion for fault activation.

[0062] In practice, the simulation results of each field are fed back to the dynamic activation criterion in real time, realizing dual dynamic feedback of pressure perturbation and chemical weakening. The specific logic is as follows: (1) Mechanics-hydrology coupling: Carbon dioxide injection → increase in pore pressure P in the seepage field → effective stress in the mechanical field s Reduction and redistribution of stress field → actual shear stress of fault t Changes → Substitute into the dynamic activation criterion to evaluate stability; (2) Chemical-mechanical coupling: carbon dioxide transport → chemical field carbon dioxide concentration distribution change → fault filling material is corroded by carbonic acid → strength parameters C(t) and μ(t) decay over time → critical shear stress of dynamic activation criterion decreases → fault stability is reassessed. (3) Closed-loop iteration: The evaluation results of each criterion are fed back to the seepage field. If the fault has a high risk of activation, the injection pressure boundary conditions are automatically adjusted to realize the closed loop of simulation-evaluation-feedback.

[0063] Furthermore, based on the results of the fully coupled simulation, the fault activation probability at different time points is calculated, and a calculation model for the fault activation probability is defined according to the dynamic activation criterion:

[0064] in, P a This represents the probability of tomographic activation, ranging from 0 to 1. A higher value indicates a higher risk of activation.

[0065] In this embodiment of the invention, an implicit iterative solution method is adopted to solve the seepage field, mechanical field, chemical field, and temperature field sequentially. This includes: within the current time step, first solving the seepage field to obtain the pore pressure distribution; then solving the mechanical field based on the pore pressure distribution to obtain the stress-strain state; then solving the chemical field based on the updated porosity and permeability to obtain the mineral dissolution / precipitation amount; and finally solving the temperature field to obtain the effect of temperature change on fluid viscosity and rock thermal stress. After each complete four-field fully coupled iteration, the maximum norm of the relative change of each field variable is calculated. When the maximum norm is less than the preset convergence tolerance, the next time step is entered.

[0066] In practice, upon entering the current simulation time step, the pore pressure, stress, mineral concentration, and temperature field converged in the previous time step are used as initial guesses. First, the seepage field solver is started, and the Darcy flow equation is solved using a fully implicit finite element method based on the currently stored permeability tensor and fluid property parameters to obtain the pore pressure distribution across all grid nodes.

[0067] The obtained pore pressure distribution is used as the basis for calculating body forces and effective stress, and the mechanics solver module is invoked. This module solves the momentum balance equations based on the nonlinear elastoplastic constitutive relationship of the rock, outputting the normal stress, shear stress, and strain tensor of each element. Simultaneously, the equivalent plastic strain of the elements is updated according to the stress state, and the corrections for porosity and permeability to the current stress level are recalculated.

[0068] The updated porosity, permeability, and stress-strain state are transferred to the chemical field solver. This solver simulates the carbon dioxide-water-rock interaction, calculates the dissolution and precipitation rates of each mineral phase, and updates the concentrations of each ion and the volume fraction of minerals in the aqueous solution. The chemical field solution results are further fed back to correct the porosity and permeability values, providing a basis for updating physical properties for the next round of seepage calculations, thus achieving two-way coupling between chemistry and seepage.

[0069] Based on the updated fluid density, specific heat capacity, and thermal conductivity, the temperature field solver is activated. This solver simultaneously considers thermal convection and thermal conduction effects, and introduces the heat of carbon dioxide dissolution and the heat of mineral reaction as source terms to calculate the temperature distribution at the current time step. The temperature field results are then used to update the fluid viscosity, rock thermal expansion strain, and reaction rate constant, achieving a complete closed loop across the four fields of heat, fluid, solidification, and chemical processes.

[0070] After solving the four fields sequentially, calculate the maximum norm of the change in the main variables of each field relative to the previous iteration in the current iteration. Compare this maximum norm with the preset convergence tolerance field by field. If the norm of change of all field variables is less than the corresponding tolerance, the iteration at the current time step is considered to have converged, the result is stored, and the process is advanced to the next time step; otherwise, based on the updated results of the four fields, repeat the above steps to start the next round of coupled iteration.

[0071] To ensure iteration efficiency, a field variable relaxation factor is used to control the update amplitude and suppress oscillations and divergence before each iteration. Simultaneously, the iteration count at the current time step is recorded. If the number of iterations is too small and convergence is rapid, the next time step size is appropriately increased; if the number of iterations is too large or a divergence trend appears, the time step size is reduced and the calculation is recalculated.

[0072] Once the iteration converges, a global check of energy and mass conservation is automatically performed to check whether the overall deviation of each field variable is within the engineering tolerance range. If the check passes, the pore pressure isosurface, stress concentration zone, mineral alteration zone, and temperature gradient distribution of the current time step are output as intermediate results. If the check fails, the step size of the time step is forcibly halved and the entire coupled solution process is re-executed until the physical consistency requirements are met.

[0073] Then, based on the fault activation probability, the fault activation risk level is determined, and a spatiotemporal cloud map of fault activation risk evolution is generated based on the fault activation risk level.

[0074] In practice, fault risk is divided into four levels based on the probability of fault activation, corresponding to different colors in the cloud map: No risk: ,blue; Low risk: ,green; Medium risk: ,yellow; High risk: ,red.

[0075] Furthermore, it outputs risk evolution cloud maps for key time nodes such as 1 year, 5 years, 10 years, and 20 years, marking the spatial distribution, area size, and evolution trend of high / medium risk zones at each time node, highlighting the risk changes of hidden faults.

[0076] In one implementation, after generating the spatiotemporal cloud map of fault activation risk evolution, the above also includes: based on the spatial distribution characteristics of high-risk areas in the spatiotemporal cloud map of fault activation risk evolution, outputting a monitoring deployment scheme and an injection pressure control strategy.

[0077] In practical implementation, the output monitoring deployment scheme and injection pressure control strategy include: taking the spatial distribution centroid of the high-risk area as the benchmark, deploying at least two rows of microseismic monitoring sensors and pore pressure monitoring sensors along the fault strike and perpendicular to the fault strike to form a cross monitoring grid; determining the coverage radius of the monitoring sensors according to the boundary range of the high-risk area, and using the particle swarm optimization algorithm to optimize the sensor deployment location, with minimizing the monitoring blind area as the objective function, and outputting the optimized monitoring deployment scheme.

[0078] The output injection pressure control strategy includes: acquiring the spatiotemporal cloud map of fault activation risk evolution at the current time step in real time, and extracting the average activation probability and maximum activation probability of high-risk areas; when the maximum activation probability exceeds a first threshold, reducing the bottom-hole injection pressure of the injection well to 80% of the current pressure; when the average activation probability exceeds a second threshold and the maximum activation probability exceeds the first threshold, alternately executing injection pressure reduction operation and injection rate reduction operation until the maximum activation probability falls back below the first threshold; wherein, the first threshold is higher than the second threshold.

[0079] In another implementation, based on the risk evolution cloud map, specific engineering guidance schemes are provided: Monitoring deployment: Monitoring points (such as microseismic monitoring, pressure monitoring, and carbon dioxide concentration monitoring) are deployed in the fault-reservoir contact zone in high / medium risk areas. The monitoring point density is ≥1 point / km² in high-risk areas and ≥1 point / 5km² in medium-risk areas. Injection pressure control: For high-risk areas, set a dynamic threshold for injection pressure (not exceeding 60%-70% of the formation fracture pressure); for areas where the risk continues to rise, gradually reduce the injection flow rate (reduction ≥20%).

[0080] To facilitate understanding, the embodiments of this invention also provide specific examples, selecting a terrestrial sedimentary CO2 geological sequestration area as the research object. This area covers approximately 100 km², and the sequestration layer is sandstone (burial depth 1500 m, formation temperature 60℃, original formation pressure 15 MPa). Two conventional faults and three hidden faults (mudstone smearing and filling) were identified. The designed CO2 injection rate is 5 × 10⁻⁶. 4 m³ / d, injection cycle 20 years. The method of this invention and a traditional static model were used to assess the risk of fault activation in this region, respectively, to verify the effectiveness of this invention.

[0081] Based on indoor experiments, the core parameter values ​​for this study area were obtained:

[0082]

[0083] As can be seen from the above, the method provided by the embodiments of the present invention can accurately identify the activation risk of hidden faults, predict the delayed activation of faults in advance, and the simulation results are more consistent with the actual underground process. It can also provide specific monitoring deployment and injection pressure control schemes for engineering practice. Compared with traditional static models, it has significant technical advantages, and its feasibility and effectiveness have been fully verified.

[0084] The method provided in this invention introduces a two-way mechanical-chemical coupling in the stability assessment of carbon dioxide storage faults, and introduces chemical corrosion weakening as an independent time variable into the activation criterion, rather than the traditional method that only considers the pressure threshold. At the same time, in response to the difficulty in identifying hidden faults and their easy corrosion characteristics, a coupling parameter system (α, β values) specific to hidden faults is established, filling the gap in the existing technology for assessing the activation risk of hidden faults and realizing coupling simulation of all fault types.

[0085] The core parameters of fault strength (cohesion C and friction coefficient μ) were improved from traditional static constants to dynamic variables controlled by two variables: carbon dioxide concentration and immersion time. An exponentially decaying chemical weakening function was obtained through indoor experiments, clarifying the parameter decay difference between hidden faults and conventional faults, making the strength parameters more consistent with the actual underground chemical corrosion process.

[0086] An implicit iterative solution method with full coupling of four fields in THMC was established, realizing a closed-loop dynamic simulation of carbon dioxide injection → pressure disturbance → stress field update → strength degradation → risk assessment → injection parameter feedback adjustment. This method solves the shortcomings of existing models that can only calculate stress changes in one direction, lack chemical weakening feedback, and lack engineering decision guidance, making the simulation process highly consistent with the actual physical and chemical processes underground.

[0087] In addition to the carbon dioxide geological storage modeling method provided in the foregoing embodiments, this invention also provides a carbon dioxide geological storage modeling apparatus, see [link to relevant documentation]. Figure 2 The schematic diagram shown illustrates the structure of a carbon dioxide geological storage modeling device, indicating that the device mainly comprises the following parts: The model building module 201 is used to build a three-dimensional geological grid model of the target area and assign physical attribute parameters to each unit of the three-dimensional geological grid model; wherein, the three-dimensional geological grid model includes conventional fault units and implicit fault units.

[0088] The dynamic criterion determination module 202 is used to acquire the experimental parameters of the hidden fault units and the conventional fault units determined by indoor experiments, and to establish a dynamic criterion for fault activation based on the experimental parameters and the actual shear stress of the fault.

[0089] The numerical simulation module 203 is used to perform multi-physics field coupled numerical simulation based on physical property parameters, and to generate a spatiotemporal cloud map of fault activation risk evolution based on the simulation results and fault activation dynamic criteria.

[0090] The carbon dioxide geological storage modeling device provided in this embodiment of the invention constructs a three-dimensional geological grid model containing conventional fault units and hidden fault units, and determines the experimental parameters of hidden fault units and conventional fault units through indoor experiments. It establishes a dynamic criterion for fault activation by combining the actual shear stress of the fault. Finally, it generates a spatiotemporal cloud map of fault activation risk evolution by using the simulation results of multi-physics coupled numerical simulation and the dynamic criterion for fault activation, thereby improving the accuracy of predicting the activation risk of hidden faults.

[0091] In one implementation, the model construction module 201 is specifically used to: acquire basic data of the target area, and extract micro-structural anomaly features from imaging logging, micro-fracture surface features from core samples, and weak reflection phase axis misalignment features from seismic data based on the basic data; identify hidden fault units based on the micro-structural anomaly features from imaging logging, micro-fracture surface features from core samples, and weak reflection phase axis misalignment features from seismic data, combined with stratigraphic sedimentary facies analysis; and establish a three-dimensional geological grid model containing stratigraphic units, reservoir units, conventional fault units, and hidden fault units using the finite element mesh generation method.

[0092] In one embodiment, the above-mentioned indoor experiment includes: under simulated formation temperature and pressure conditions, after continuously exposing fault rock samples to a brine environment containing carbon dioxide for different time periods, measuring their cohesion and friction coefficient respectively; different time periods covering the time scales of the initial, middle and late injection stages; and using a double-shear friction test device, measuring the activation critical shear stress of conventional fault samples and latent fault samples under different normal stresses.

[0093] In one embodiment, the dynamic criterion determination module 202 is specifically used to: fit the measured cohesion and friction coefficient using an exponential decay model to obtain the experimental relationship that the cohesion and friction coefficient of the hidden fault unit and the conventional fault unit respectively show a decay trend with carbon dioxide exposure time; and fit the normal stress and the activation critical shear stress to establish the fitting relationship between the critical shear stress and the normal stress.

[0094] In one embodiment, the numerical simulation module 203 is specifically used to: solve the seepage field, mechanical field, chemical field and temperature field of the target region in sequence using an implicit iterative solution method, and obtain the simulation results of each field.

[0095] In one embodiment, the numerical simulation module 203 is specifically used to: update the fault activation dynamic criterion based on the simulation results, and calculate the fault activation probability at different time points based on the updated fault activation dynamic criterion; determine the fault activation risk level based on the fault activation probability, and generate a spatiotemporal cloud map of fault activation risk evolution based on the fault activation risk level.

[0096] In one embodiment, the above-mentioned device further includes: an adjustment module, used to: output a monitoring deployment scheme and an injection pressure control strategy based on the spatial distribution characteristics of high-risk areas in the spatiotemporal cloud map of fault activation risk evolution.

[0097] It should be noted that the device provided in this embodiment of the invention has the same implementation principle and technical effects as the aforementioned method embodiment. For the sake of brevity, any parts not mentioned in the device embodiment can be referred to the corresponding content in the aforementioned method embodiment. The specific numerical values ​​provided in this embodiment are merely exemplary and are not intended to limit the scope of the invention.

[0098] This invention also provides an electronic device, specifically, the electronic device includes a processor and a storage device; the storage device stores a computer program, and the computer program, when run by the processor, executes the method described in any of the above embodiments.

[0099] Figure 3This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention. The electronic device 100 includes: a processor 30, a memory 31, a bus 32 and a communication interface 33. The processor 30, the communication interface 33 and the memory 31 are connected through the bus 32. The processor 30 is used to execute executable modules, such as computer programs, stored in the memory 31.

[0100] The memory 31 may include high-speed random access memory (RAM) or non-volatile memory, such as at least one disk storage device. Communication between this system network element and at least one other network element is achieved through at least one communication interface 33 (which can be wired or wireless), such as the Internet, wide area network, local area network, metropolitan area network, etc.

[0101] Bus 32 can be an ISA bus, PCI bus, or EISA bus, etc. The bus can be divided into address bus, data bus, control bus, etc. For ease of representation, Figure 3 The symbol is represented by a single double-headed arrow, but this does not mean that there is only one bus or one type of bus.

[0102] The memory 31 is used to store programs. After receiving an execution instruction, the processor 30 executes the program. The method executed by the device for defining the flow process disclosed in any of the foregoing embodiments of the present invention can be applied to the processor 30 or implemented by the processor 30.

[0103] Processor 30 may be an integrated circuit chip with signal processing capabilities. In implementation, each step of the above method can be completed by the integrated logic circuitry in the hardware of processor 30 or by instructions in software form. Processor 30 can be a general-purpose processor, including a Central Processing Unit (CPU), a Network Processor (NP), etc.; it can also be a Digital Signal Processor (DSP), an Application Specific Integrated Circuit (ASIC), a Field-Programmable Gate Array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components. It can implement or execute the methods, steps, and logic block diagrams disclosed in the embodiments of this invention. The general-purpose processor can be a microprocessor or any conventional processor. The steps of the methods disclosed in the embodiments of this invention can be directly embodied in the execution of a hardware decoding processor, or executed by a combination of hardware and software modules in the decoding processor. The software modules can reside in random access memory, flash memory, read-only memory, programmable read-only memory, electrically erasable programmable memory, registers, or other mature storage media in the art. The storage medium is located in memory 31. The processor 30 reads the information in memory 31 and, in conjunction with its hardware, completes the steps of the above method.

[0104] The computer program product of the readable storage medium provided in the embodiments of the present invention includes a computer-readable storage medium storing program code. The instructions included in the program code can be used to execute the methods described in the foregoing method embodiments. For specific implementation, please refer to the foregoing method embodiments, which will not be repeated here.

[0105] If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this invention, essentially, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0106] Finally, it should be noted that the above-described embodiments are merely specific implementations of the present invention, used to illustrate the technical solutions of the present invention, and not to limit it. The scope of protection of the present invention is not limited thereto. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that any person skilled in the art can still modify or easily conceive of changes to the technical solutions described in the foregoing embodiments within the technical scope disclosed in the present invention, or make equivalent substitutions for some of the technical features; and these modifications, changes, or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should all be covered within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for modeling carbon dioxide geological storage, characterized in that, include: A three-dimensional geological grid model of the target area is constructed, and physical property parameters are assigned to each unit of the three-dimensional geological grid model; wherein, the three-dimensional geological grid model includes conventional fault units and implicit fault units; The experimental parameters of the hidden fault unit and the conventional fault unit, determined by indoor experiments, are obtained, and a dynamic criterion for fault activation is established based on the experimental parameters and the actual shear stress of the fault. Multiphysics field coupling numerical simulation is performed based on the physical property parameters, and a spatiotemporal cloud map of fault activation risk evolution is generated based on the simulation results and the fault activation dynamic criterion.

2. The method according to claim 1, characterized in that, Construct a three-dimensional geological mesh model of the target area, including: Acquire basic data of the target area, and extract micro-structural anomaly features, micro-fracture surface features of core samples, and weak reflection phase axis misalignment features of seismic data based on the basic data; Based on the micro-structural anomaly characteristics of the imaging logging, the micro-fracture surface characteristics of the core, and the weak reflection phase axis misalignment characteristics of the seismic data, combined with stratigraphic sedimentary facies analysis, hidden fault units are identified. A three-dimensional geological grid model containing stratigraphic units, reservoir units, conventional fault units, and hidden fault units was established using the finite element mesh generation method.

3. The method according to claim 1, characterized in that, The indoor experiments included: under simulated formation temperature and pressure conditions, fault rock samples were continuously exposed to a carbon dioxide-containing brine environment for different time periods, and their cohesion and friction coefficient were measured respectively; the different time periods covered the time scales of the initial, middle and late stages of injection. Using a double-shear friction test apparatus, the activation critical shear stress of conventional fault samples and latent fault samples was determined under different normal stresses.

4. The method according to claim 3, characterized in that, Obtaining experimental parameters of the hidden fault unit and the conventional fault unit determined through indoor experiments includes: The measured cohesion and friction coefficient were fitted using an exponential decay model to obtain the experimental relationship between the cohesion and friction coefficient of the hidden fault unit and the conventional fault unit, showing a decay trend with carbon dioxide exposure time. The normal stress and the activation critical shear stress are fitted to establish a fitting relationship between the critical shear stress and the normal stress.

5. The method according to claim 1, characterized in that, Multiphysics coupling numerical simulation based on the aforementioned physical property parameters includes: An implicit iterative solution method is used to solve the seepage field, mechanical field, chemical field and temperature field of the target region in sequence, and the simulation results of each field are obtained.

6. The method according to claim 1, characterized in that, Based on the simulation results and the aforementioned fault activation dynamic criteria, a spatiotemporal cloud map of fault activation risk evolution is generated, including: The fault activation dynamic criterion is updated based on the simulation results, and the fault activation probability at different time points is calculated based on the updated fault activation dynamic criterion. Based on the fault activation probability, the fault activation risk level is determined, and a spatiotemporal cloud map of fault activation risk evolution is generated based on the fault activation risk level.

7. The method according to claim 1, characterized in that, After generating the spatiotemporal cloud map of fault activation risk evolution, the following steps are also included: Based on the spatial distribution characteristics of high-risk areas in the spatiotemporal cloud map of fault activation risk evolution, a monitoring deployment scheme and injection pressure control strategy are output.

8. A carbon dioxide geological storage modeling device, characterized in that, include: The model building module is used to construct a three-dimensional geological grid model of the target area and assign physical property parameters to each unit of the three-dimensional geological grid model; wherein, the three-dimensional geological grid model includes conventional fault units and implicit fault units; The dynamic criterion determination module is used to acquire the experimental parameters of the hidden fault unit and the conventional fault unit determined by indoor experiments, and to establish a dynamic criterion for fault activation based on the experimental parameters and the actual shear stress of the fault. The numerical simulation module is used to perform multi-physics field coupled numerical simulation based on the physical property parameters, and to generate a spatiotemporal cloud map of fault activation risk evolution based on the simulation results and the fault activation dynamic criterion.

9. An electronic device, characterized in that, The method includes a processor and a memory, the memory storing computer-executable instructions executable by the processor, the processor executing the computer-executable instructions to implement the steps of the method according to any one of claims 1 to 7.

10. A computer-readable storage medium storing a computer program thereon, characterized in that, The computer program, when run by a processor, performs the steps of the method described in any one of claims 1 to 7.

Images