Carbonate rock crack modeling method based on digital outcrop and geomechanical simulation

Abstract

The invention discloses a carbonate rock crack modeling method based on digital outcrop and geomechanical simulation. According to the method, firstly, a high-precision digital outcrop model is constructed through unmanned aerial vehicle oblique photography, outcrop surface crack parameters are obtained, and crack groups are divided; and determining the ancient stress field characteristics of multiple fracture forming periods by combining filler isotope test and acoustic emission experiment. And further, carrying out geomechanical simulation based on the recovered paleostructure, and quantitatively calculating the density and opening parameters of the crack in each stage by utilizing a two-stage fracture criterion and a strain energy density theory. And then, a deterministic modeling technology and a random modeling technology are fused, a multi-stage and multi-scale discrete fracture network model is constructed, and fluid simulation is performed based on the model so as to analyze the control effect of the fracture on the ancient karst. And finally, through comparison verification with outcrop measured data, the corrected crack-matrix coupling prototype geologic model is output, quantitative prediction of the cracks in the outcrop is realized, and the precision and reliability of the crack prototype model are remarkably improved.

Description

Technical Field

[0001] This invention belongs to the field of petroleum development, specifically relating to a method for constructing a prototype model of carbonate rock fractures driven by digital outcrops and geomechanical simulation. Background Technology

[0002] Carbonate reservoirs in the Tarim Basin are currently a hotspot for oil and gas exploration and development, accounting for approximately 38% of the basin's total oil and gas resources. Tectonic fractures and karst fracture-cavity systems, as important structural components of carbonate reservoirs, not only provide migration pathways for oil and gas but also serve as crucial storage spaces. Therefore, clarifying the development patterns and distribution laws of fracture-karst evolution is essential for oil and gas exploration in carbonate reservoirs. With the continuous advancement of carbonate oil and gas reservoir development, the demand for accurate geological models in oilfield exploration and development operations is increasing daily. However, due to the heterogeneous nature of reservoirs, relying solely on seismic data and well logging data interpretation is insufficient to accurately characterize reservoir fractures. To overcome the limitations of subsurface data in providing detailed descriptions of reservoir fractures, geologists have recently attempted to use methods such as outcrop observation, testing and analysis, and the construction of prototype models to provide reasonable parameters for reservoir fracture modeling. A prototype model refers to a detailed reservoir model of outcrops, densely developed well networks in mature oilfields, or modern sedimentary environments that are similar in characteristics to the target reservoir area. Practice has proven that outcrops are the most authentic and direct source of data. Building prototype models and geological knowledge bases based on outcrop data is the most suitable choice for assisting in reservoir fracture modeling. However, traditional outcrop prototype models generally employ measurement descriptions, sampling analysis, and ground-penetrating radar (GPR). These methods are mostly two-dimensional representations, offering limited assistance in reservoir fracture modeling. With the development of surveying technology, three-dimensional digital outcrop technology based on UAV scanning has emerged, which can effectively help identify and statistically analyze fracture parameters on the outcrop surface. However, digital outcrop models cannot reflect the fracture characteristics within the rock mass. Current methods such as GPR have low detection accuracy and limited application conditions, making it difficult to detect and identify fractures throughout the entire outcrop rock mass. Therefore, research is needed from the perspectives of geological analysis, geomechanics, and tectonic stress fields to predict and characterize fractures within the outcrop rock mass, and to combine this with digital outcrop modeling to constrain and guide the construction of carbonate rock prototype models. Summary of the Invention

[0003] To address the aforementioned problems, this invention provides a method for constructing a prototype model of carbonate rock fractures driven by digital outcrops and geomechanical simulation. Specifically, for the problem of outcrop prototype model construction, a three-dimensional fracture modeling method driven by both digital outcrop models and geomechanical theory is developed. Fracture data on the outcrop surface is measured using the digital outcrop model, and fracture data inside the outcrop is predicted using geomechanical theory. Based on the three-dimensional digital model of a typical outcrop, modeling data and fracture parameters are extracted. Through comprehensive geological analysis and testing, the geometric distribution and developmental stages of the fractures are analyzed. Numerical simulation methods are used to simulate the stress field of the study area in multiple phases and further calculate fracture parameters. With the data support of the digital outcrop model and the guidance of geomechanical theory, a multi-phase, multi-scale fracture discrete model construction technique is employed. Based on this, multi-phase seepage simulation is performed to reconstruct paleofluid dissolution and seepage characteristics, restore the development and distribution patterns of fractures and paleocavities, and establish a prototype geological model of outcrop fractures.

[0004] The present invention adopts the following technical solution:

[0005] A modeling method for carbonate rock fractures based on digital outcrops and geomechanical simulation, such as Figure 1 The process, as shown, includes the following steps:

[0006] S1. Constructing a digital outcrop model of a typical outcrop: Based on UAV oblique photogrammetry technology, the selected outcrop area is scanned and three-dimensionally reconstructed to obtain a digital outcrop model, and geometric information of strata, faults and outcrop surface cracks is extracted from it.

[0007] S2. Statistical analysis of outcrop crack parameters: Combining digital outcrop models with field identification, the orientation, length and aperture parameters of outcrop surface cracks are statistically analyzed, and different crack groups are classified according to their orientation and mechanical characteristics.

[0008] S3. Determine the main fracture formation period of the outcrop crack: Determine the filling period and the formation sequence of the crack by isotope testing and structural trace analysis of the crack filling material, and obtain the direction and magnitude of paleostress in each period by combining acoustic emission experiments.

[0009] S4. Paleostress field simulation during major fracture formation periods: A geomechanical model is established based on the restored paleotectonic model and rock mechanical parameters. The paleostress data obtained in S3 is used as boundary conditions to simulate the distribution of paleostress fields during each major fracture formation period.

[0010] S5. Quantitative characterization of multi-stage fracture parameters: Based on the stress field results of S4, the rock fracture criterion is used to identify fracture, and the density and aperture of each stage of fracture are calculated according to the quantitative relationship model between stress-strain and fracture parameters.

[0011] S6. Multi-stage, multi-scale crack discrete modeling: Deterministic modeling based on digital outcrop data is used for large-scale cracks, and stochastic modeling with crack density obtained from S5 as constraint is used for medium and small-scale cracks. Finally, the models are superimposed to form a multi-stage crack discrete network model.

[0012] S7. Fluid simulation based on discrete fracture network: Fluid flow simulation is performed based on the discrete fracture network model in S6 to calculate fracture permeability and analyze the control effect of multi-stage fractures on the development of paleokarst.

[0013] S8. Verification and correction of the outcrop prototype model: Compare and correct the predicted attributes of the fracture model with the measured outcrop data, and output the verified fracture-matrix coupled prototype geological model.

[0014] Preferably, step S1, constructing a digital outcrop model of a typical outcrop, includes the following steps:

[0015] First, through field investigation, typical outcrop areas with similarities to the oilfield reservoir and suitable for constructing a prototype model were selected. To construct a prototype geological model of the outcrops, it is necessary not only to collect complete and large-scale outcrop data, but also to collect outcrop data with sufficient accuracy to identify fractures. Considering the needs of outcrop fracture research, two scanning methods were employed based on UAV oblique photogrammetry: constant-altitude aerial photography and close-range aerial photography. First, constant-altitude scanning was used to cover the distribution area of ​​the main structures in the outcrop area. Then, close-range scanning was conducted on valley-type profiles with fracture development to construct a high-precision digital model of the outcrops for fracture statistical analysis.

[0016] The fixed-altitude scanning method involves the UAV acquiring data at a predetermined fixed flight altitude. By planning a reasonable flight path for the target area, holographic data of the outcrop is scanned. The model constructed using this method can perform three-dimensional feature analysis of the overall topography, strata, and fault structures of the outcrop area, with a resolution of 0.1–0.3 m. The close-up scanning method employs a serpentine flight path close to the profile, maintaining a certain and relatively short distance between the target geological body and the camera. To obtain higher quality image data, the distance between the UAV and the geological body is generally within 15 meters, with variable flight altitude. This method is suitable for typical profiles with well-developed fractures, acquiring higher quality image data, with a model resolution of 3–6 mm. The specific process of three-dimensional data acquisition of outcrops includes: (1) field survey of outcrops and preparation of instruments, manual survey and point determination of outcrops in the study area; (2) instrument assembly and testing, including checking the camera parameters of the UAV and the camera gimbal, etc.; (3) field outcrop image acquisition, using flight control software DJI GO 4 to monitor the status of the UAV and the image acquisition situation, and flying and shooting according to the planned design path; (4) data testing and point cloud data processing, recovering the UAV and transmitting the acquired data, checking the image quality, and performing post-processing of the acquired data; (5) establishing a digital outcrop model through professional software.

[0017] Digital outcrop model data processing mainly includes data extraction, point cloud data processing, and assignment of a 3D spatial coordinate system. The model data collected by the UAV primarily includes 2D point cloud data, 3D point cloud data, and image data. The UAV continuously takes pictures during its cruise, recording coordinate information and storing it on a memory card, ensuring a good match between the collected image data and the point cloud coordinate information. Using Context Capture software, the processed point cloud is transformed into a 3D TIN grid model, and surface texture mapping is applied to assign color information, establishing a prototype digital outcrop model. The digital outcrop model can be viewed from any angle, overcoming the limitations of human visual observation of cross-sections. This facilitates the observation of lithological characteristics and structural topography, as well as the extraction of data such as crack length and density, significantly improving the efficiency and quantification of field investigations.

[0018] Preferably, step S2, the statistical analysis of outcrop crack parameters, includes the following specific steps:

[0019] First, during field investigations, the basic characteristics of cracks, such as their type, morphology, and mechanical properties, are identified. Natural structural cracks generally have the following characteristics: (1) Structural cracks generally exist in groups. Cracks generated under the same tectonic stress are generally parallel or conjugate and have regional distribution characteristics. Cracks caused by human destruction are often centered on the point of destruction and exist in local areas; (2) Natural cracks have been preserved for a long time and are mostly filled to varying degrees; (3) Natural structural cracks usually extend for a long time and extend into the interior of the rock. They are generally divided into tensile cracks and shear cracks. Shear cracks are the most common cracks in outcrops. They are relatively straight and there is generally slight displacement between their fracture surfaces. Some fracture surfaces have traces such as striations and steps. Tensile cracks, on the other hand, have irregular serrated traces. Based on the experience of recognizing outcrops in the field, three-dimensional visualization software such as Acute3D Viewer is used to identify cracks on digital outcrop models. The software's measurement and positioning tools are used to identify and quantify the structural features of faults and fractures, and to identify and statistically analyze the attitude, length, aperture, and filling status of outcrop fractures.

[0020] A fracture system refers to a collective term for fractures or fracture groups generated under the same tectonic activity. Based on the fracture orientation and mechanical characteristics, outcrop fractures are divided into multiple fracture systems, each representing a different period of tectonic movement. There are generally relationships such as cut-off and termination between superimposed fracture systems; later-stage fractures may terminate their development at earlier-stage fractures, and earlier-stage fractures may be cut off and displaced by later-stage fractures. Taking early and late-stage fractures as examples, the superimposed fracture relationships can be summarized into four modes: oblique intersection, displacement, confinement, and slip. The oblique type refers to a fractured rock mass that fractures again during intense tectonic activity in later periods, resulting in new fractures that cut through existing fractures, but the shear stress is insufficient to displace the original fractures. The displacement type refers to new fractures generated by later tectonic activity that show obvious shear slip marks, cutting through and displaced from the original fractures. The confined type refers to a situation where, when the existing fractures are large, the stress generated by later tectonic activity is blocked, and the developing new fractures are confined to one side of the existing fractures, failing to penetrate them and developing intermittently on both sides. The sliding type refers to a situation where, under strong later tectonic activity, although the late-stage fractures are confined by the existing fractures, the residual tectonic stress causes the existing fractures to reactivate, resulting in sliding failure along the early fracture surface. Through the analysis of fracture combination relationships, the superposition sequence of fractures in the outcrops is further analyzed.

[0021] Preferably, step S3, determining the main crack formation period of the exposed crack, includes the following steps:

[0022] Carbonate rock fractures are generally highly infilled, suggesting they may have undergone infilling processes at different times. The development time of fractures can be inferred through analysis of the fracture infill material. Based on the analysis of the superposition relationship of fracture strata, isotopic analysis was performed on fracture infill materials from different strata. Following the research of Wei Yi and Hofs (1976), etc., the δ¹⁸O⁻ of Ordovician marine limestone generally... 13 C PDB Values ​​range from -5 to 5‰, δ 18 O PDB The values ​​range from -6 to -10‰. Combining the method used by Keith and Weber (1964) to calculate paleosalinity using carbon and oxygen isotopes, a carbon and oxygen isotope cross-section diagram was plotted based on the test results to zonate the infill material. A paleosalinity greater than 120 indicates a marine sedimentary environment, while less than 120 indicates a freshwater sedimentary environment. Sr isotopes can also reflect the infill formation environment. Typically, during the infiltration of freshwater from the atmosphere, the fluid carries a certain amount of... 87 Sr, leading to 87 Sr / 86 Sr value increases and δ 13 The C-value is negative. Deep hydrothermal fluids may carry [something] during their ascent. 87 Sr caused contamination of the filling material and 87 Sr / 86 The Sr value increased. By comparing the composition of the protolith and paleofluid environment, the sedimentary environment of the fracture infill in different zones was analyzed. Combined with the tectonic background, the formation sequence of the fracture infill was clarified. Based on the results of the fracture assemblages superposition sequence analyzed in step 2, the formation period of the outcrop fractures was determined.

[0023] Acoustic emission (AE) technology, also known as the acoustic emission method, can record the stresses previously experienced by rocks. It is highly practical and applicable to various geological conditions, including outcrops and downholes. When a rock is subjected to external forces, its internal structure develops corresponding microcracks. When the rock is subjected to external forces again, if the force is less than the magnitude of the initial force, the existing microcracks will not further fracture. However, if the force exceeds the magnitude of the initial force, the existing microcracks will expand again. This phenomenon is called the Kaiser effect, and the stress at this point represents the magnitude of the external force that previously caused the microcracks. Based on the characteristics of the Kaiser effect, a triaxial compression tester can be used to continuously apply loads to the rock. An acoustic emission receiver and stress sensor are then used to receive the crack expansion and the emitted acoustic emissions, thereby inferring paleostress. Fresh, crack-free rock samples were taken for acoustic emission experiments. Four standard plunger samples conforming to national standards were collected. Figure 2Three plunger samples were parallel to the horizontal plane with axial angles of 45 degrees, while one plunger sample was perpendicular to the horizontal. For paleostress measurement, the response characteristics of the cumulative acoustic emission count curve need to be statistically analyzed. Each significant steep increase in slope on the response curve needs to be statistically analyzed. To ensure experimental accuracy, multiple repeated experiments are required at the same sampling point. In addition to comparing multiple experiments, statistical comparisons of different tectonic locations on the surface are also necessary. Under the influence of a certain tectonic movement, the stress field will exhibit consistent characteristics over a large range; therefore, sampling at multiple locations within the study area is required. Based on the slope variation characteristics of the cumulative acoustic emission count curve, Kaiser effect points and corresponding memory stress values ​​are identified and read. The rock sample test results are statistically analyzed, the frequency of different memory stress values ​​is calculated, and a histogram of memory stress frequency is plotted. Different peak intervals of memory stress represent different tectonic movements. Based on the memory stress in different directions, the average value is taken to calculate the magnitude and direction of the principal stresses of different tectonic movements. Underground rock masses are typically subjected to principal stresses in three directions. Based on the memory stresses in each direction obtained from acoustic emission experiments, the average value is taken, and the magnitude and direction of the principal stresses are calculated using the following formula:

[0024]

[0025] in, , and These are the stress values ​​(MPa) obtained from horizontal sampling tests at 0°, 45°, and 90°. The angle between the direction of the maximum principal stress and the due north direction is expressed in degrees. It is the maximum principal stress in the horizontal direction, in MPa; The minimum principal stress in the horizontal direction is denoted as MPa. Based on the formation law of conjugate shear cracks and the orientation characteristics of outcrop cracks, and matching the direction of the maximum principal stress calculated from the acoustic emission test results, the stress field staging results of cracks at various stages and the acoustic emission test results are coupled to analyze the tectonic stress field characteristics of the main crack-forming period, and determine the magnitude of the principal stress and the direction of the maximum principal stress.

[0026] Preferably, step S4, paleostress field simulation during the main fracture formation period, includes the following steps:

[0027] Based on the digital outcrop model scanned by UAV aerial photography in step 1, point cloud coordinate information is extracted by forming a 3D TIN grid model from the point cloud. This point cloud is then imported into Petrel modeling software to construct the current outcrop surface model. Combined with lithological identification through field investigation, the interface locations of various strata in the outcrop area are determined. Relevant coordinate information of the stratigraphic interfaces is identified and extracted through digital outcrop observation. Based on the coordinate information of the stratigraphic layers, models of each stratigraphic interface are constructed in the modeling software. Combined with the outcrop surface model, a surface reflecting the current real-world topographic relief of the outcrop is constructed, restoring the current outcrop stratigraphic characteristics. Based on the stratigraphic faulting and valley morphology on the digital outcrop model, and combined with the results of field investigation, faults and large-scale fractures are identified. Based on the stratigraphic relief and fault faulting relationships, the weathered and eroded stratigraphic model is restored and completed. The original stratigraphic model before weathering and erosion is then restored in the modeling software. Based on the outcrop-restored stratigraphic model and fault data, 3DMove is used to reconstruct the tectonic evolution, determine the fault tectonic evolution characteristics of different fracture formation stages, and combine the stratigraphic model with faults of different stages to construct tectonic geological models of different periods for paleostress field simulation modeling of different stages.

[0028] Paleostress field reconstruction is fundamental to fracture prediction. After determining the main fracture formation period, it is necessary to simulate the paleostress field at the time of fracture formation. Using the reconstructed paleotectonic surfaces and faults, a corresponding paleostress field finite element model for the fracture formation period is constructed on ANSYS finite element simulation software. A combination of dense and sparse meshing is used, with dense meshing in key areas of complex structures such as fractures and slightly coarser meshing in gentler, less critical areas. Combining lithological analysis and rock mechanics experiments, the lithological characteristics and rock mechanics parameters of each stratum are determined. The rock mechanics parameters of each stratum in the model are then divided and assigned values, and combined with the finite element model to construct an outcrop geomechanical model. After the outcrop geomechanical model is constructed, based on the paleostress characteristics analyzed by acoustic emission experiments in step 3, boundary conditions for stress field simulation are set, and paleostress field simulation of the main fracture formation period is performed. Finally, stress field data and stress distribution cloud maps characterized in the form of nodes or elements can be output. Each node and element includes coordinate, stress, and strain data.

[0029] Preferably, step S5, quantitative calculation of crack parameters, includes the following steps:

[0030] The characteristics of crack density distribution are important parameters for constraining discrete crack modeling. Rock fracturing under stress generally follows certain patterns. The formation and evolution of tectonic cracks refer to the process by which the internal structure of a rock gradually fractures under external stress loads due to shearing or tension. In nature, different stress fields with different characteristics, such as direction and magnitude, produce cracks with varying orientations and types. Previous experiments on numerous rocks have shown that the fracturing process in rocks under external compressive loads generally goes through several stages, typically including a compaction stage, an elastic deformation stage, and a plastic deformation stage. Carbonate rocks, having undergone long-term diagenesis, have low porosity and are stronger and more brittle than dense sandstone, resulting in slightly different fracturing evolution patterns. Based on comprehensive analysis of stress-strain curves, acoustic emission signals, and permeability variation curves, the fracturing evolution stages of carbonate rocks are divided into an initial compaction stage, an elastic strain stage, a plastic deformation stage, and a fracture propagation stage. The threshold ranges for each stage are determined through experimental results. Taking the Ordovician carbonate rocks of the Tarim Basin as an example, the initial compaction stage is generally when the loaded stress is less than 10% of the load limit. The existing cracks in the rock mass close, and the fractures have not yet developed. The elastic strain stage is when the loaded stress is between 10% and 75% of the load limit. The rock undergoes elastic strain, and a small number of microcracks develop, which have little effect on the transformation of oil and gas porosity and permeability. The plastic deformation stage is when the loaded stress exceeds 75% of the load limit. Microcracks begin to develop in large numbers, and the rock undergoes irreversible plastic deformation, and the permeability of the rock begins to increase. The fracture propagation stage is when the loaded stress exceeds a certain proportion of the load limit (83% under low confining pressure or 91% under high confining pressure). The number of microcracks reaches its peak, begins to merge and expand into larger-scale fractures, and finally forms macroscopic fractures, leading to rock mass fracturing and instability.

[0031] Generally, due to different stress field states, the mechanical properties of cracks also differ. The most common structural cracks are divided into tensile cracks and shear cracks. Therefore, in actual crack prediction, it is necessary to select an appropriate fracture criterion for prediction and judgment. Based on the basic principles of rock fracture criteria, a fracture criterion for carbonate rocks considering confining pressure was determined through multiple trials. Under triaxial compressive stress, the two-stage Coulomb-Mohr criterion is used as the judgment criterion for shear cracks, while the Griffith criterion is used when tensile cracks occur under extensional stress.

[0032] According to the Coulomb-Mohr criterion, when rock undergoes shear fracture, the shear stress on the fracture surface must reach the shear strength of the rock mass, and the friction between the fracture surfaces is also relevant. The magnitude of this friction arises from the compressive effect of the normal stress. Therefore, it can be said that under triaxial compression, when rock undergoes shear fracture at a certain angle, there is a certain functional relationship between the shear stress and the normal stress:

[0033]

[0034] In the formula, τ is the shear stress on the fracture surface, MPa; C is the cohesion of the rock, MPa; f is the internal friction coefficient of the rock, which can be defined as f=tan φ, where φ is the internal friction angle of the rock; σ is the normal stress on the fracture surface, MPa. Based on the point of tangency between the stress Mohr's circle and the fracture envelope, the formulas for calculating the normal stress and shear stress on the fracture surface when shear fracture occurs can be derived:

[0035]

[0036] In the formula, α is the angle between the fracture surface and the minimum principal stress, in °; σ1 is the maximum principal stress on the rock as a whole, in MPa; σ3 is the minimum principal stress on the rock as a whole, in MPa. There is a certain angular relationship between α and φ:

[0037]

[0038] Under triaxial compression conditions, rock fracture satisfies the two-stage Mohr-Coulomb criterion. Combining the results of triaxial mechanical experiments on carbonate rocks, and simplifying equations 5 and 6 into equation 4, we can obtain:

[0039]

[0040] Because carbonate rocks exhibit different properties under different confining pressures, stress Mohr's circles and their envelopes were plotted based on triaxial mechanical experiments. The envelope morphology was extracted and presented as a two-segment Mohr-Coulomb curve. Two internal friction angles, φ1 and φ2, can be obtained from the curve, with their slopes (i.e., internal friction coefficients) being k1 and k2, respectively. The confining pressure boundary value is σ0. When 0 < σ3 < σ0, the internal friction angle is φ = φ1; when σ3 ≥ σ0, the internal friction angle is φ = φ2. The two-segment Mohr-Coulomb criterion can reasonably determine whether rock fracture has occurred during shear fracture and can indicate the fracture direction.

[0041] Griffith proposed the concept of the Griffith fracture, believing that rocks inherently possess weak points and micro-fractures, which he abstracted as elliptical fractures. When rocks are subjected to tensile stress, stress concentration occurs at the ends of the Griffith fracture. When the local stress exceeds the tensile strength of the rock, it begins to propagate at its ends. Under sustained tensile stress, the fracture gradually expands, leading to rock instability and ultimately the formation of a tensile crack. Therefore, according to Griffith's strength theory, the tensile stress state can be summarized into the following two cases:

[0042] when At that time, the rupture criterion is:

[0043]

[0044] when At that time, the rupture criterion is:

[0045]

[0046] In the formula: σ t Tensile strength of rock, expressed in Pa; This is the tensile fracture angle.

[0047] Based on the paleostress field simulation results, under compressive stress (i.e., when σ3 ≥ 0), the two-stage Coulomb-Mohr fracture criterion is used to determine the stress state and fracture status at the nodes in the finite element simulation. Under tensile stress (i.e., when σ3 < 0), the Griffith fracture criterion is used to determine the fracture status at the nodes. If the node is determined to be unfractured, the crack parameter is set to 0. If fracture is determined, the crack parameter is calculated using the following crack parameter calculation model considering confining pressure.

[0048] According to relevant theories of elasticity and the maximum strain energy density theory, brittle rocks accumulate elastic strain energy under external stress. When the rock fractures, part of this accumulated energy is used to generate the crack surface, while the remainder is released as elastic waves. The stress reaches the load limit σ. c Microcracks begin to develop extensively when the elastic strain energy density ω needs to be overcome before crack formation occurs, which corresponds to 85% of the total elastic strain energy density. e The total strain energy density ω of the rock can be expressed by the following equation:

[0049]

[0050] In the formula, ω is the total strain energy density of the rock, J / m 3 ;ω f The strain energy density required for crack initiation, J / m 3 σ1, σ2, and σ3 represent the maximum, intermediate, and minimum principal stresses experienced by the rock, respectively, in MPa; ε1, ε2, and ε3 represent the strain of the rock in the directions of the maximum, intermediate, and minimum principal stresses, respectively. Compared with dense sandstone and carbonate rocks, the Ordovician carbonate rocks in the Tarim Basin are more brittle and stronger, with more easily developed fractures. Their fracture characteristics and key stress values ​​differ from those of brittle sandstone. Based on the analysis of carbonate rock fracture evolution, after the elastic strain stage, the rock can be further subdivided into a plastic deformation stage and a fracture propagation stage. During the plastic deformation stage, microcracks begin to appear in large numbers, corresponding to the 0.85σ1 stress of brittle sandstone. c The stress critical point is the stage after the formation of microcracks, but this node appears earlier in carbonate rocks, with the stress critical point at approximately 0.75σ. c .

[0051] In summary, the deformation and fracture development characteristics of carbonate rocks differ under low and high confining pressure conditions. Furthermore, considering the two-segment Mohr-Coulomb curves obtained from rock mechanics experiments, the mechanical properties of the rock change significantly before and after the confining pressure reaches σ0. Therefore, the calculation of fracture parameters in carbonate rocks will be discussed in two separate cases:

[0052] (1) Low confining pressure condition (σ3<σ0)

[0053] At this stage, the strain energy generated by the plastic deformation is relatively small, and most of it is related to the final rock mass fracture. Therefore, this portion of the strain energy can be attributed to the surface energy of the fracture released during the final fracture, hence the following relationship:

[0054]

[0055] Under triaxial compressive stress, based on elastic theory and rock fracture evolution analysis, 0.75σ can be calculated. p ω at time e :

[0056]

[0057] In the formula, σ p It is the magnitude of the stress that causes rock fracture, in MPa; σ 0.75p , ε 0.75p The axial stress of the rock is at 0.75σ. p The magnitude of stress and strain.

[0058] Combining equations 11 and 12 above, we can simplify to:

[0059]

[0060] In the formula, E is the elastic modulus of the rock (GPa); μ is the Poisson's ratio of the rock.

[0061] Under tensile stress, rock fractures instantaneously after reaching peak strength, thus requiring an elastic strain energy density ω to overcome before cracks can form. e To achieve the tensile strength σ t The strain energy density corresponding to this time is:

[0062]

[0063] (2) High confining pressure condition (σ3≥σ0)

[0064] At this point, crack formation requires overcoming not only elastic strain energy but also plastic strain energy. According to rock fracture evolution analysis, the stress in the plastic strain stage is at 0.75σ. p ~0.91σp Since this stage is short-lived, the formation of cracks requires overcoming the strain energy of this stage.

[0065]

[0066] In the formula, ω p The strain energy density consumed by the rock during the plastic strain stage, in J / m³. 3 Through stress-strain curve analysis, this curve can be approximated as a straight line, and its slope can be defined as E. b The strain energy density ω during the plastic strain stage p This can be represented by the following formula:

[0067]

[0068] In the formula, σ 0.91p , ε 0.91p The axial stress of the rock is at 0.91σ. p The magnitude of stress and strain.

[0069] Substituting equation 13 into equation 12, we get:

[0070]

[0071] Under triaxial compressive stress (σ3 ≥ 0), based on the paleostress field simulation results, the rock fracture is assessed, and the crack density and aperture are calculated using the following formulas:

[0072]

[0073] Under the condition of tensile stress (σ3<0), the calculation methods for parameters such as crack density and aperture are as follows:

[0074]

[0075] In the formula, D lf The calculated crack linear density parameter is given in terms of cracks / m; D vf For the calculated crack density parameters, m 2 / m 3 C0 represents the cohesion of the rock, in MPa; J represents the energy required to generate a crack per unit area, in J / m³. 2 J0 is the surface energy of the crack under zero confining pressure, in J / m³. 2 b is the crack opening, in meters; L1 and L3 are the lengths of the characterizing unit along the directions σ1 and σ3, respectively, in meters; ε is the maximum tensile strain under the current stress state, dimensionless; ε0 is the maximum elastic tensile strain, dimensionless; θ is the fracture angle, in degrees.

[0076] In Equation 19, there are two cases regarding the calculation of crack linear density parameters:

[0077] when When, then the fracture angle is

[0078]

[0079] when When, then the fracture angle is

[0080]

[0081] That is, the volume density and linear density parameters of the crack are the same.

[0082] The calculated fracture parameters can be imported into the outcrop strata model in the form of node distribution to construct fracture density and fracture aperture models, thereby showing fracture development characteristics and providing constraint properties for the construction of discrete fracture models.

[0083] Preferably, step S6, multi-stage, multi-scale crack discrete modeling, includes the following steps:

[0084] While fracture prediction based on tectonic stress fields can characterize the distribution features of fracture parameters, it lacks intuitiveness and cannot meet the quantitative and intuitive needs of prototype models. In contrast, the Discrete Fracture Network (DFN) model is more intuitive, combining fracture parameters from geomechanical simulations as constraints for fracture modeling, and can better represent fracture distribution and development characteristics. Different fracture modeling methods were adopted to model fracture grids of different phases within the outcrop area, considering the development characteristics of fractures at different scales. Large fractures exceeding 100 meters in size and penetrating the entire stratum are easily identified in outcrops. Therefore, deterministic modeling techniques were used for large fracture modeling. High-precision coordinate data of large-scale fractures at different phases were systematically extracted from the digital outcrop model using image processing and 3D reconstruction techniques, including parameters such as fracture development location, attitude, and length. Within the framework of the existing outcrop strata feature model, the coordinates and geometric feature parameters of large-scale fractures were inserted, directly transforming them into corresponding large-scale fracture sections, thus completing the modeling of large-scale fractures at different phases.

[0085] Because small and medium-sized fractures have a small spatial scale, stochastic modeling methods can be used to characterize their geometry and distribution. This method requires constraints based on fracture attribute distribution patterns to randomly generate the center locations of small and medium-sized fractures, and then assigning fracture attributes to construct corresponding fracture models. Therefore, geostatistical understanding is needed to integrate constraints such as fracture intensity (density), size distribution, dip / strike, and spatial clustering patterns to generate a statistically representative fracture network. Fracture measurements based on field outcrops and digital models can determine the range of key geometric parameters, including attitude (strike, dip, azimuth) and size (length, aperture), for subsequent DFN modeling. However, fractures observed in outcrops only represent the development characteristics of exposed fractures on the outcrop surface and are not suitable for determining the spatial distribution of fractures inside the outcrop. Geomechanical simulation methods can effectively characterize the spatial distribution attributes of small and medium-scale fracture density and can serve as constraint models for fracture modeling. In the stochastic modeling process, multiple mathematical models and attribute data need to be integrated to systematically constrain key geometric parameters. Fracture density, as a constraint model defining fracture intensity distribution, provides a three-dimensional spatial distribution of fracture density through geomechanical simulation. The accuracy of the results is ensured through comparison with field measurements. Based on outcrop fracture statistics, the specific proportion of medium- and small-scale fractures can be determined. This proportion is used to determine the number and density of medium- and small-scale fractures in the model, allowing for separate modeling. A Fisher distribution model is employed, providing probabilistic quantification parameters to characterize the azimuth distribution of fracture attitude (such as strike and dip). Simultaneously, a power-law function is used to characterize the size distribution of fracture length for modeling. By comprehensively representing fracture networks at different scales and directions using multiple factors and models, a discrete model of medium- and small-scale fractures integrating multivariate data is established. Within the grid framework of the outcrop geological model, large, medium, and small fractures are superimposed to construct a multi-scale discrete fracture model. Following this approach, fractures of different phases can be modeled separately using geomechanical simulation results and then superimposed to obtain the current outcrop fracture network model.

[0086] Preferably, step S7, crack network seepage simulation, includes the following steps:

[0087] Fluid simulation based on crack discrete network

[0088] Fluid simulation based on the DFN model yields the permeability distribution characteristics of the fracture network. Combining this with the permeability characteristics of the fracture network at different stages allows for further analysis of the control effect of fractures on karst development at different times. The cubic law is a fundamental model in subsurface fluid mechanics, idealizing fractures as parallel plates and deriving flow equations under the assumption of laminar flow through planar pores. The basis for permeability calculations primarily stems from the geometric characteristics of fractures and fluid flow simulations. Based on the Oda method and the cubic law, a permeability tensor is constructed. The calculations describe the laminar flow of fluid along a smooth crack system:

[0089]

[0090] in, The symbol for Kronecker; , , , represent the projections of the crack normal vector in the i and j directions, respectively. By coarsening the discrete crack network (DFN) model, key crack properties, including crack line density and connectivity, can be systematically extracted from the scaled-up mesh. Subsequently, single-phase flow simulations using the finite element method (FEM) or finite volume method (FVM) can calculate the equivalent permeability tensor in different coordinate directions, providing crucial input for studying the heterogeneity of fluid flow.

[0091] The methods described above can be used to determine the permeability distribution of fracture networks in different directions, thereby describing the heterogeneity of fluid flow at different times. Using fluid flow modeling to study the genetic relationship between fractures and paleokarst aims to provide a basis for predicting complex karst reservoirs in the subsurface environment.

[0092] Preferably, step 8, verification and correction of the prototype model, includes the following steps:

[0093] Based on the existing outcrop geological model mesh, the fracture discretization model is coarsened to obtain fracture-matrix coupling coefficient and fracture density attributes. The fracture discretization model and matrix model can serve as prototype models for outcrop fractures. The fracture mesh and attribute model data are projected onto the existing outcrop surface, and typical outcrop profiles are selected for statistical analysis of fracture density and orientation. This data is then compared and validated with the fracture density and orientation data in the prototype model. The outcrop prototype model is corrected by adjusting the boundary conditions of the geomechanical simulation to maintain basic consistency with the outcrop fracture characteristics. Furthermore, the reliability of fluid seepage simulation and fracture network can be further verified by comparing the outcrop paleokarst characteristics with the permeability characteristics of fracture networks from different periods.

[0094] Compared with the prior art, the beneficial effects of this invention are:

[0095] This invention discloses a modeling method for carbonate rock fractures based on digital outcrops and geomechanical simulation. Based on geomechanical theory, it forms a quantitative prediction method for carbonate rock fractures that considers the influence of confining pressure, improving the accuracy of carbonate rock fracture parameter characterization under complex stress conditions. By integrating digital outcrop models and geomechanical simulation methods, it develops a multi-stage, multi-scale fracture discrete grid model construction technology, which extends the quantitative characterization of outcrop fractures from the surface to the interior of the rock, solving the problem of difficult characterization of fractures inside carbonate rock outcrops. This provides a new approach for digital characterization of outcrops and construction of fracture prototype models, and provides technical support for the construction of a carbonate reservoir geological knowledge base and the exploration and development of oil and gas reservoirs. Attached Figure Description

[0096] Figure 1 Flowchart of the method for constructing a prototype model of fractures in carbonate rocks;

[0097] Figure 2 Schematic diagram of acoustic emission experiment sampling;

[0098] Figure 3 Statistical analysis of parameters for exposed cracks in a single room;

[0099] Figure 4 The superposition relationship of cracks in different groups;

[0100] Figure 5 A cross-sectional diagram of carbon and oxygen isotope content;

[0101] Figure 6 A cross-sectional diagram of strontium and carbon isotope content;

[0102] Figure 7 Statistical histograms of memory stress in rocks from different periods;

[0103] Figure 8 Boundary conditions were set for the fourth phase of paleostress field simulation.

[0104] Figure 9 The distribution characteristics of the crack density simulation results in four phases;

[0105] Figure 10 The diagram shows the current multi-stage discrete model of cracks, which is a superposition of four crack models.

[0106] Figure 11 Figure showing the simulation results of seepage in a fracture network;

[0107] Figure 12 Comparison chart for verifying crack patterns and crack density in the prototype model. Detailed Implementation

[0108] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below. However, it should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of the invention. Furthermore, descriptions of well-known structures and technologies are omitted in the following description to avoid unnecessarily obscuring the concept of the invention.

[0109] Example:

[0110] This embodiment uses an outcrop of Ordovician carbonate rocks in northwestern Xinjiang as an example to illustrate the specific technical solution of the present invention:

[0111] The Ordovician strata in the Tarim Basin of Xinjiang are characterized by multiple phases of tectonic superposition. Intense karstification results in extremely heterogeneous reservoirs, making the construction of reservoir fracture models and prototype models very difficult, posing significant challenges to oil and gas reservoir exploration and development. The widespread outcrops of the Ordovician strata in northwestern Tarim Basin, with their prevalent multi-phase fractures, solution pores, and caves, provide the most direct data samples for the study of carbonate rock fractures and dissolution, offering new insights for the research of multi-phase fractures and the construction of prototype models. The Yijianfang outcrop, a typical example, was selected as the research sample for an outcrop fracture prototype geological model. This outcrop exhibits multi-phase fractures and significant karstification, and establishing a prototype model based on it can provide important reference value for reservoir fracture model construction and the evolution of fracture-karst systems.

[0112] First, construct digital outcrop models of typical outcrops.

[0113] The northwestern Tarim Basin, nestled against the southern edge of the Tianshan Mountains, is a typical outcrop area of ​​Ordovician carbonate strata. Through multiple field investigations in the northwestern Tarim Basin, the Yijianfang outcrop area was selected as the research target area. It shares a similar tectonic location and dynamic background with the North-Central Tarim Oilfield, exhibiting the same Ordovician carbonate strata. The outcrop fractures and karst development patterns also provide valuable references for studying underground reservoir fractures, making it suitable for constructing a prototype geological model. Ordovician hydrocarbons in the Tarim Basin are primarily preserved in reef-shoal deposits of carbonate platform facies. Based on field outcrop observations, the reef-shoal deposits in the Yijianfang area mainly consist of the Yingshan Formation, the Yijianfang Formation, and the Tumushuke Formation, which are the primary target strata for this modeling study.

[0114] To construct a prototype geological model of the outcrop, it is necessary not only to collect complete and extensive outcrop data, but also, considering the needs of outcrop fracture research, to collect outcrop data with sufficient accuracy to identify fractures. Based on UAV oblique photogrammetry technology, two scanning methods were employed: constant-altitude aerial photography and close-range aerial photography. First, constant-altitude scanning was used to cover the distribution area of ​​the main structures in the Yijianfang outcrop area, constructing a centimeter-level precision digital model of the outcrop for geological modeling. Then, a close-range scanning was conducted on a valley-type profile with fracture development in the Yijianfang outcrop to construct a millimeter-level precision digital model of the outcrop, used for statistical analysis of fracture and cave characteristics.

[0115] Second, statistical analysis of outcrop crack parameters.

[0116] Based on experience in identifying outcrops in the field, 3D visualization software such as Acute3D Viewer was used to identify cracks on digital outcrop models. The software's measurement and positioning tools were employed to identify and quantify the structural features of faults and cracks, and to identify and statistically analyze the attitude, length, aperture, and filling status of outcrop cracks. Figure 3 The cross-section of Yijianfang shows a zigzag valley formed by faults trending approximately 210°, 270°, and 250°. Statistical analysis of the cracks within these valleys revealed that cracks trending approximately 210° and 250° were the most numerous, with dip angles around 70°, classifying them as high-angle cracks and the primary cause of valley formation. Cracks trending approximately 270° and 330° were the next most numerous, with dip angles ranging from 80° to 85°, classifying them as vertical cracks. These cracks extend in the same direction as the caves and are closely related to their dissolution and formation. In terms of crack density, the linear density of cracks in the Yijianfang area ranges from 0.5 to 2.6 cracks per meter. Overall, the crack development exhibits segmented characteristics, with higher linear density in the southwest and lower in the northeast. The southwest segment has the most cracks trending 210°, 250°, and 270°, while the northeast segment is dominated by cracks trending 330° and 210°. In terms of crack scale, most of the exposed cracks are less than 10 meters in length and less than 2 mm in diameter. Overall, small and medium-sized cracks dominate, while large-scale cracks are present but not numerous. In terms of crack filling characteristics, non-structural cracks are basically fully filled, more than half of the structural cracks show a high degree of filling, exhibiting full-fill to half-fill, while some late-developing cracks are basically unfilled.

[0117] For cases of multiple overlapping fractures, the stage can be determined by factors such as truncation, displacement, and the degree of filling. A fracture system refers to a collective term for fractures or fracture groups generated under the same tectonic activity. Based on the fracture's attitude and mechanical characteristics, outcrop fractures are divided into multiple fracture systems, each representing a different period of tectonic movement. There are generally truncation or termination relationships between the superposition of different fracture systems; later-stage fractures may terminate their development at earlier-stage fractures, and earlier-stage fractures may be cut off and displaced by later-stage fractures. Through analysis of fracture combination relationships, the superposition sequence of multiple fracture systems in the outcrop can be further analyzed.

[0118] Based on the crack orientation, the outcrop cracks in the northwestern part of the Tarim Basin can be basically divided into groups of 330°, 270°, 250°, 210°, 300°, and 0° cracks. Among them, the number of cracks with orientations of 210° and 250° is roughly equal, their dip angles are similar, and their crack characteristics are similar, forming a conjugate group of cracks. The 300° and 0° cracks have low filling degrees, and through multi-stage crack intersection analysis, they also form a conjugate group of cracks. An outcrop in one room shows that the 330° cracks are mostly interrupted by the 270° and 250° cracks, indicating that this is a relatively early-developing group of cracks. Figure 4 -a); the 270° crack was interrupted by conjugate cracks at 250° and 210° ( Figure 4 -b) indicates that the 270° crack developed earlier, representing the second-stage crack series after the 330° crack series. Based on the combined analysis of multiple crack superpositions in the outcrop, the crack series development is preliminarily ranked from ancient to modern as follows: 330° crack series (first-stage cracks), 270° crack series (second-stage cracks), 250° and 210° crack series (third-stage cracks), and 300° and 0° crack series (fourth-stage cracks).

[0119] Step 3: Determine the main crack formation period of the exposed crack.

[0120] Carbonate rock fractures are generally highly infilled, and these fractures may have undergone infilling processes at different times. The development time of fractures can be inferred by analyzing the fracture infill. The Ordovician strata in the Tarim Basin have experienced multiple phases of tectonic uplift and subsidence. After the Late Permian, the Tarim Plate underwent large-scale marine transgression and regression, and then transitioned to a terrestrial sedimentary stage. The aquatic environment changed, and the characteristics of the infill within the fractures of the strata also differed.

[0121] Based on the analysis of the superposition relationship of fracture assemblages, isotopic tests were conducted on fracture filling materials from different assemblages. According to previous studies by Wei Yi and Hofs (1976), the δ¹⁸O values ​​of Ordovician marine limestone are generally... 13 C PDB Values ​​range from -5 to 5‰, δ 18 O PDBValues ​​range from -6 to -10‰. Carbon and oxygen isotope analysis was performed on Ordovician carbonate rocks and cavity fillings collected from the Yijianfang area. The δ¹⁰ values ​​of the protoliths... 13 C PDB The value ranges from -0.5 to 0.3‰, δ 18 O PDB The values ​​range from -7.9 to -9.2‰, and the elemental characteristics of some cavern and fissure fillings are similar to these values. Combining the paleosalinity calculated using carbon and oxygen isotopes by Keith and Weber (1964), a carbon and oxygen isotope cross-section diagram was drawn based on the test results, dividing the fissure fillings into four parts (…). Figure 5 The paleosalinity of zones I, II, and III is greater than 120, indicating a marine environment. Zone IV has a paleosalinity less than 115, and the infill material is freshwater-originating calcite veins. The 330° and 270° fracture series infill materials are mainly distributed in zones I, II, and III, which are marine sedimentary environments, while the 250° and 210° fracture series infill materials are distributed in zone IV, which is a freshwater environment. The infill material in zone I is similar to the protolith, with only the 330° fracture infill material, representing the earliest developed fracture. Based on the tectonic background, it is further inferred that during the Caledonian Orogeny, the northwestern Tarim Basin was a marine carbonate platform, and the 330° and 270° fracture assemblages were filled after their appearance. During the Hercynian Orogeny, the Bachu Uplift gradually rose, and the sedimentary environment began to change from marine to terrestrial. During this period, volcanic and magmatic activities were frequent, and the early fractures were filled with hydrothermal material. After the Late Permian, the sedimentary environment of the outcrop area became completely terrestrial. The 250° and 210° fracture assemblages did not have marine-originating filling material, so they developed after this period.

[0122] One room area 87 Sr / 86 Sr values ​​range from 0.7085 to 0.7106. Figure 6 Overall, the cracks and cavern fillings 87 Sr / 86 The Sr values ​​were higher than those measured in the original rock. The high values ​​are mainly due to the large amount of Sr carried by the fluid as it passed through the siliceous debris deposits. 87 Sr. Typically, during the infiltration of freshwater into the atmosphere, the fluid carries a certain amount of... 87 Sr, leading to 87 Sr / 86 Sr value increases and δ 13 The C value is negative, while the test results show that δ 13 The C value is mostly similar to that of the original rock, and atmospheric freshwater infiltration is only part of the filling material. 87 Sr / 86The increase in Sr value is not the primary cause of this phenomenon in the test sample. The infill material used in the test was calcite and fluorite from a fluorite cave profile. Fluorite is formed by the combination of fluoride ions and calcium ions released from underground hydrothermal fluids. Therefore, the high Sr value of the fluorite infill material is not the main reason for this phenomenon. 87 Sr / 86 The Sr value likely originated from hydrothermal upwelling, carrying sediment from Precambrian clastic rocks or bedrock. 87 Sr caused cracks in the fluorite cave profile and the interior filling material, although δ 13 The C values ​​are generally similar to those of the original rock, but 87 The phenomenon of high Sr.

[0123] Based on the above analysis of the superposition sequence of fracture assemblages, the formation stages of the outcrop fractures were determined. During the Middle Caledonian period, there was relatively strong compression and uplift in the outcrop area, and the earliest 330° fracture assemblages developed during this period. During the Late Caledonian to Early Hercynian periods, even more intense compression and uplift occurred. At this time, the 270° fracture series formed and was uplifted to near the surface, resulting in intense karstification. The filling of the fracture cavities was influenced by atmospheric freshwater. During the Late Hercynian, volcanic activity and magma upwelling activated some pre-existing faults, and hydrothermal fluids intruded along faults and fractures, leading to the filling of some early fractures with hydrothermal materials such as fluorite and sulfur. During the Late Hercynian to Indosinian periods, the strata underwent another period of compression and uplift, ending the marine sedimentary environment. The filling of the resulting 250° and 210° fracture series occurred primarily in a freshwater environment; therefore, no significant hydrothermal filling material was observed. During the Himalayan orogeny, after a brief period of Paleogene sedimentation, the outcrop area experienced renewed tectonic activity, exposing the Ordovician strata to the surface. During this period, the 300° and 0° fracture series were formed. These fractures had a short formation time and were exposed to the surface for a long period. The karst water had low calcium carbonate saturation, resulting in weak filling and a low degree of fracture filling.

[0124] Different fracture assemblages correspond to paleostress fields of different periods. The direction of the maximum principal stress during the fracture formation period is inferred from the fracture tectonic traces and shear angle characteristics, and the magnitude of paleostress is then estimated by combining this with acoustic emission experiments. Acoustic emission technology, also known as the AE method, can record the stress that rocks have experienced. It is highly practical and applicable to different geological conditions, such as outcrops or downholes. For paleostress measurement, it is necessary to statistically analyze the response characteristics of the cumulative number of acoustic emissions, and to statistically analyze each significant point of steep slope increase on the response curve. To ensure the accuracy of the experiment, multiple repeated experiments are required at the same sampling point. According to the principle of the multi-stage nature of the Kaiser effect in stratigraphic rocks, older stratigraphic rocks record more tectonic movement stages. Previous studies have conducted multiple acoustic emission experiments on northwestern Tarim Basin and surrounding areas, including experiments on Cenozoic, Mesozoic, and Paleozoic stratigraphic rocks. The rock experimental data can be used as control data. Therefore, Middle and Upper Ordovician rocks in the outcrop area were selected as samples. In order to better analyze the paleotectonic stress of the Ordovician strata in the northwestern Tarim Basin, it is necessary to select fresh rocks without obvious cracks for sampling experiments.

[0125] Based on the slope variation characteristics of the cumulative acoustic emission curve, Kaiser effect points were identified and read. Most samples exhibited 3-4 Kaiser effect points during testing, indicating that the outcrop rock had undergone at least 3-4 phases of tectonic movement. The test results of rock samples from different horizontal directions were statistically analyzed to calculate the frequency of different memory stress intervals in several samples, and a histogram of memory stress frequency was plotted. Figure 7 Statistical results show that memory stress has four peak ranges, including 40-60 MPa, 60-80 MPa, 80-100 MPa and 140-160 MPa, which correspond to the early Himalayan orogeny, the Indosinian-Yanshanian tectonic movement and two Paleozoic tectonic movements in the study area, respectively.

[0126] Different peak ranges of memory stress represent different tectonic movements. Based on the memory stress in different directions, the magnitude and direction of the principal stresses for different tectonic movements are calculated. Subsurface rock masses are typically subjected to principal stresses in three directions. Based on the memory stresses in each direction statistically analyzed by acoustic emission experiments, the average value is taken, and the magnitude and direction of the principal stresses are calculated using the following formula:

[0127]

[0128] in, , and These are the stress values ​​(MPa) obtained from horizontal sampling tests at 0°, 45°, and 90°. The angle between the direction of the maximum principal stress and the due north direction is expressed in degrees. It is the maximum principal stress in the horizontal direction, in MPa; The minimum principal stress in the horizontal direction is , MPa.

[0129] Based on the formation law of conjugate shear cracks and the orientation characteristics of outcrop cracks, the direction of the maximum principal stress calculated from the acoustic emission experimental results was matched, and the stress field staging results of cracks at various stages were coupled with those from the acoustic emission experiments. The experimental analysis is consistent with the crack orientation characteristics, and the stress magnitude and direction characteristics are shown in Table 1. This result can provide a reference for setting the boundary conditions in the fourth step of paleostress simulation.

[0130]

[0131] Step 4: Paleostress field simulation during the main fracture formation period

[0132] Based on the digital outcrop model scanned by UAV aerial photography in step 1, a 3D TIN grid model was formed from point clouds. Point cloud coordinate information was extracted and imported into Petrel modeling software to construct a current outcrop surface model of the Yijianfang area. Combined with lithological identification through field investigation, the interface locations of various strata on the outcrop surface were determined. Relevant coordinate information of the stratigraphic interfaces was collected through digital outcrop observation. Based on the coordinate information of the stratigraphic interfaces, models of each stratigraphic interface were constructed in Petrel modeling software. Combined with the outcrop surface model, a 3D model reflecting the current real outcrop topography and stratigraphic information was constructed, restoring the current outcrop stratigraphic characteristics. Based on the stratigraphic faulting and valley morphology on the digital outcrop model, combined with the results of field investigation, faults and large-scale fractures were identified. Based on the stratigraphic undulations and fault faulting relationships, the stratigraphic model before weathering and erosion was restored and completed. The original stratigraphic model before weathering and erosion was restored in the modeling software. Based on the outcrop-reconstructed stratigraphic model and fault data, 3DMove was used to reconstruct the tectonic evolution, determine the fault tectonic evolution characteristics of different fracture formation stages, and combine the stratigraphic model with faults of different stages to construct a tectonic geological model corresponding to the four fracture formation stages of the Yijianfang outcrop area, which was used for paleostress field simulation modeling of different stages.

[0133] Paleostress field reconstruction is fundamental to fracture prediction. After determining the main fracture formation period, it is necessary to simulate the paleostress field at the time of fracture formation. Following the determination of the fracture formation period, four structural geological models were established based on outcrop digital models. Structural surfaces and fault data were extracted from these models, and corresponding paleostress field finite element models for the fracture formation period were constructed using ANSYS finite element simulation software. A combination of dense and sparse meshing was employed for mesh generation. The target strata in this study include the Tumushuke Formation, Yijianfang Formation, and the upper section of the Yingshan Formation exposed in the Yijianfang area. These strata are also the main underground reservoirs of this oilfield.

[0134] Based on lithological analysis and rock mechanics experiments of different strata in the Yijianfang area, the lithological characteristics and rock mechanics parameters of each stratum were determined. The rock mechanics parameters of each stratum in the model were then divided and assigned values, and combined with the finite element model to construct an outcrop geomechanical model. After the outcrop geological model was constructed, based on the paleostress characteristics analyzed by acoustic emission experiments, boundary conditions for stress field simulation were set: During the Middle Caledonian period, under the influence of the Kunlun Ocean Plate, the outcrop area was subjected to NE-SSW compressive stress. The model boundary conditions are as follows: Figure 8 -a; Late Caledonian to Early Hercynian period, inheriting the earlier kinetic background, the outcrop area is subjected to NW-SE compression, the model boundary conditions are as follows. Figure 8 -b; During the Late Hercynian-Indosinian period, the outcrop area was subjected to the combined influence of the Kunlun Plate and the Southern Tianshan Ocean, resulting in a stress field exhibiting NE-SW bilateral compression. The model boundary conditions are as follows: Figure 8 -c; During the Himalayan orogeny, the outcrop area was subjected to NNW-SSE compression due to thrusting from the southern Tianshan Mountains. The model boundary conditions are as follows: Figure 8 -d. Then, based on the finite element principle, paleostress field simulation is performed, and the final output is paleostress field data and stress distribution cloud map representing the four crack formation periods in the form of nodes or elements. Each node and element includes coordinate, stress and strain data.

[0135] Step 5: Quantitative characterization of multi-stage crack parameters

[0136] The characteristics of crack density distribution are important parameters for constraining discrete crack modeling. Rock fracturing under stress generally follows certain patterns. The formation and evolution of tectonic cracks refer to the process by which the internal structure of a rock gradually fractures under external stress loads due to shearing or tension. In nature, different stress fields with different characteristics, such as direction and magnitude, produce cracks with varying orientations and types. Previous experiments on numerous rocks have shown that the fracturing process in rocks under external compressive loads generally goes through several stages, typically including a compaction stage, an elastic deformation stage, and a plastic deformation stage. Carbonate rocks, having undergone long-term diagenesis, have low porosity and are stronger and more brittle than dense sandstone, resulting in slightly different fracturing evolution patterns. Based on comprehensive analysis of stress-strain curves, acoustic emission signals, and permeability variation curves, the fracturing evolution stages of carbonate rocks are divided into an initial compaction stage, an elastic strain stage, a plastic deformation stage, and a fracture propagation stage. The threshold ranges for each stage are determined through experimental results. Taking the Ordovician carbonate rocks of the Tarim Basin as an example, the initial compaction stage is generally when the loaded stress is less than 10% of the load limit. The existing cracks in the rock mass close, and the fractures have not yet developed. The elastic strain stage is when the loaded stress is between 10% and 75% of the load limit. The rock undergoes elastic strain, and a small number of microcracks develop, which have little effect on the transformation of oil and gas porosity and permeability. The plastic deformation stage is when the loaded stress exceeds 75% of the load limit. Microcracks begin to develop in large numbers, and the rock undergoes irreversible plastic deformation, and the permeability of the rock begins to increase. The fracture propagation stage is when the loaded stress exceeds a certain proportion of the load limit (83% under low confining pressure or 91% under high confining pressure). The number of microcracks reaches its peak, begins to merge and expand into larger-scale fractures, and finally forms macroscopic fractures, leading to rock mass fracturing and instability.

[0137] Generally, due to different stress field states, the mechanical properties of cracks also differ. The most common structural cracks are divided into tensile cracks and shear cracks. Therefore, in actual crack prediction, it is necessary to select an appropriate fracture criterion for prediction and judgment. Based on the basic principles of rock fracture criteria, a fracture criterion for carbonate rocks considering confining pressure was determined through multiple trials. Under triaxial compressive stress, the two-stage Coulomb-Mohr criterion is used as the judgment criterion for shear cracks, while the Griffith criterion is used when tensile cracks occur under extensional stress.

[0138] According to the Coulomb-Mohr criterion, when rock undergoes shear fracture, the shear stress on the fracture surface must reach the shear strength of the rock mass, and the friction between the fracture surfaces is also relevant. The magnitude of this friction arises from the compressive effect of the normal stress. Therefore, it can be said that under triaxial compression, when rock undergoes shear fracture at a certain angle, there is a certain functional relationship between the shear stress and the normal stress:

[0139]

[0140] In the formula, τ is the shear stress on the fracture surface, MPa; C is the cohesion of the rock, MPa; f is the internal friction coefficient of the rock, which can be defined as f=tan φ, where φ is the internal friction angle of the rock; σ is the normal stress on the fracture surface, MPa. Based on the point of tangency between the stress Mohr's circle and the fracture envelope, the formulas for calculating the normal stress and shear stress on the fracture surface when shear fracture occurs can be derived:

[0141]

[0142] In the formula, α is the angle between the fracture surface and the minimum principal stress, in °; σ1 is the maximum principal stress on the rock as a whole, in MPa; σ3 is the minimum principal stress on the rock as a whole, in MPa. There is a certain angular relationship between α and φ:

[0143]

[0144] Under triaxial compression conditions, rock fracture satisfies the two-stage Mohr-Coulomb criterion. Combining the results of triaxial mechanical experiments on carbonate rocks, and simplifying equations 5 and 6 into equation 4, we can obtain:

[0145]

[0146] Because carbonate rocks exhibit different properties under different confining pressures, stress Mohr's circles and their envelopes were plotted based on triaxial mechanical experiments. The envelope morphology was extracted and presented as a two-segment Mohr-Coulomb curve. Two internal friction angles, φ1 and φ2, can be obtained from the curve, with their slopes (i.e., internal friction coefficients) being k1 and k2, respectively. The confining pressure boundary value is σ0 = 30 MPa. Specifically, when 0 < σ3 < σ0, the internal friction angle is φ = φ1 = 65°; when σ3 ≥ σ0, the internal friction angle is φ = φ2 = 25°. The two-segment Mohr-Coulomb criterion can reasonably determine whether rock fracture has occurred during shear fracture and can indicate the fracture direction.

[0147] Griffith proposed the concept of the Griffith fracture, believing that rocks inherently possess weak points and micro-fractures, which he abstracted as elliptical fractures. When rocks are subjected to tensile stress, stress concentration occurs at the ends of the Griffith fracture. When the local stress exceeds the tensile strength of the rock, it begins to propagate at its ends. Under sustained tensile stress, the fracture gradually expands, leading to rock instability and ultimately the formation of a tensile crack. Therefore, according to Griffith's strength theory, the tensile stress state can be summarized into the following two cases:

[0148] when At that time, the rupture criterion is:

[0149]

[0150] when At that time, the rupture criterion is:

[0151]

[0152] In the formula: σ t Tensile strength of rock, expressed in Pa; This is the tensile fracture angle.

[0153] Based on the paleostress field simulation results, under compressive stress (i.e., when σ3 ≥ 0), the two-stage Coulomb-Mohr fracture criterion is used to determine the stress state and fracture status at the nodes in the finite element simulation. Under tensile stress (i.e., when σ3 < 0), the Griffith fracture criterion is used to determine the fracture status at the nodes. If the node is determined to be unfractured, the crack parameter is set to 0. If fracture is determined, the crack parameter is calculated using the following crack parameter calculation model considering confining pressure.

[0154] According to relevant theories of elasticity and the maximum strain energy density theory, brittle rocks accumulate elastic strain energy under external stress. When the rock fractures, part of this accumulated energy is used to generate the crack surface, while the remainder is released as elastic waves. The stress reaches the load limit σ. c Microcracks begin to develop extensively when the elastic strain energy density ω needs to be overcome before crack formation occurs, which corresponds to 85% of the total elastic strain energy density. e The total strain energy density ω of the rock can be expressed by the following equation:

[0155]

[0156] In the formula, ω is the total strain energy density of the rock, J / m 3 ;ω f The strain energy density required for crack initiation, J / m 3 σ1, σ2, and σ3 represent the maximum, intermediate, and minimum principal stresses experienced by the rock, respectively, in MPa; ε1, ε2, and ε3 represent the strain of the rock in the directions of the maximum, intermediate, and minimum principal stresses, respectively. Compared with dense sandstone and carbonate rocks, the Ordovician carbonate rocks in the Tarim Basin are more brittle and stronger, with more easily developed fractures. Their fracture characteristics and key stress values ​​differ from those of brittle sandstone. Based on the analysis of carbonate rock fracture evolution, after the elastic strain stage, the rock can be further subdivided into a plastic deformation stage and a fracture propagation stage. During the plastic deformation stage, microcracks begin to appear in large numbers, corresponding to the 0.85σ1 stress of brittle sandstone. c The stress critical point is the stage after the formation of microcracks, but this node appears earlier in carbonate rocks, with the stress critical point at approximately 0.75σ. c .

[0157] In summary, the deformation and fracture development characteristics of carbonate rocks differ under low and high confining pressure conditions. Furthermore, considering the two-segment Mohr-Coulomb curves obtained from rock mechanics experiments, the mechanical properties of the rock change significantly before and after the confining pressure reaches σ0. Therefore, the calculation of fracture parameters in carbonate rocks will be discussed in two separate cases:

[0158] (1) Low confining pressure conditions (σ3<30MPa)

[0159] At this stage, the strain energy generated by the plastic deformation is relatively small, and most of it is related to the final rock mass fracture. Therefore, this portion of the strain energy can be attributed to the surface energy of the fracture released during the final fracture, hence the following relationship:

[0160]

[0161] Under triaxial compressive stress, based on elastic theory and rock fracture evolution analysis, 0.75σ can be calculated. p ω at time e :

[0162]

[0163] In the formula, σ p It is the magnitude of the stress that causes rock fracture, in MPa; σ 0.75p , ε 0.75p The axial stress of the rock is at 0.75σ. p The magnitude of stress and strain.

[0164] Combining equations 11 and 12 above, we can simplify to:

[0165]

[0166] In the formula, E is the elastic modulus of the rock (GPa); μ is the Poisson's ratio of the rock.

[0167] Under tensile stress, rock fractures instantaneously after reaching peak strength, thus requiring an elastic strain energy density ω to overcome before cracks can form. e To achieve the tensile strength σ t The strain energy density corresponding to this time is:

[0168]

[0169] (2) High confining pressure conditions (σ3≥30MPa)

[0170] At this point, crack formation requires overcoming not only elastic strain energy but also plastic strain energy. According to rock fracture evolution analysis, the stress in the plastic strain stage is at 0.75σ.p ~0.91σ p Since this stage is short-lived, the formation of cracks requires overcoming the strain energy of this stage.

[0171]

[0172] In the formula, ω p The strain energy density consumed by the rock during the plastic strain stage, in J / m³. 3 Through stress-strain curve analysis, this curve can be approximated as a straight line, and its slope can be defined as E. b The strain energy density ω during the plastic strain stage p This can be represented by the following formula:

[0173]

[0174] In the formula, σ 0.91p , ε 0.91p The axial stress of the rock is at 0.91σ. p The magnitude of stress and strain.

[0175] Substituting equation 13 into equation 12, we get:

[0176]

[0177] Under triaxial compressive stress (σ3 ≥ 0), based on the paleostress field simulation results, the rock fracture is assessed, and the crack density and aperture are calculated using the following formulas:

[0178]

[0179] Under the condition of tensile stress (σ3<0), the calculation methods for parameters such as crack density and aperture are as follows:

[0180]

[0181] In the formula, D lf The calculated crack linear density parameter is given in terms of cracks / m; D vf For the calculated crack density parameters, m 2 / m 3 C0 represents the cohesion of the rock, in MPa; J represents the energy required to generate a crack per unit area, in J / m³. 2 J0 is the surface energy of the crack under zero confining pressure, in J / m³. 2 b is the crack opening, in meters; L1 and L3 are the lengths of the characterizing unit along the directions σ1 and σ3, respectively, in meters; ε is the maximum tensile strain under the current stress state, dimensionless; ε0 is the maximum elastic tensile strain, dimensionless; θ is the fracture angle, in degrees.

[0182] In Equation 19, there are two cases regarding the calculation of crack linear density parameters:

[0183] when When, then the fracture angle is

[0184]

[0185] when When, then the fracture angle is

[0186]

[0187] That is, the volume density and linear density parameters of the crack are the same.

[0188] Based on the above methods for calculating crack parameters, quantitative calculations of cracks are performed using restored paleostress field data to obtain the distribution characteristics of crack density in four phases, such as... Figure 9 Vertically, fractures of all phases are more developed in the Yijianfang and Yingshan Formations, and relatively less developed in the Tumushuke Formation. Horizontally, the first phase of fractures mainly developed during the Middle Caledonian period, primarily distributed in the eastern part of the outcrop area; the second phase developed during the Late Caledonian to Early Hercynian period, mainly distributed in the southern and southwestern parts of the outcrop area; the third phase developed during the Late Hercynian to Indosinian period, with the largest number of fractures, distributed in the western and southern parts of the outcrop area; and the fourth phase developed during the Himalayan period, mainly in the southeastern part of the outcrop area. The distribution characteristics of fractures in each phase are basically consistent with the outcrop development, and can be used for phased fracture modeling. The calculated fracture parameters can be imported into the outcrop stratigraphic model in the form of node distributions to construct fracture density and fracture aperture models, thereby displaying fracture development characteristics and providing constraint properties for the construction of discrete fracture models.

[0189] Step 6: Multi-stage, multi-scale crack discretization modeling

[0190] While fracture prediction based on tectonic stress fields can characterize the distribution features of fracture parameters, it lacks intuitiveness and cannot meet the quantitative and intuitive needs of prototype models. In contrast, the Discrete Fracture Network (DFN) model is more intuitive, combining fracture parameters from geomechanical simulations as constraints for fracture modeling, and can better represent fracture distribution and development characteristics. Different fracture modeling methods were adopted to model fracture grids of different phases within the outcrop area, considering the development characteristics of fractures at different scales. Large fractures exceeding 100 meters in size and penetrating the entire stratum are easily identified in the outcrop. Therefore, deterministic modeling techniques were used for large fracture modeling. Through image processing and 3D reconstruction techniques, high-precision coordinate data of large-scale fractures at different phases were systematically extracted from the digital outcrop model, including parameters such as fracture development location, attitude, and length. Taking the third-phase fractures as an example, large-scale fractures trending 250° and 210° can be identified using a digital outcrop model, and the coordinate information of the structural traces can be extracted. Within the framework of the current outcrop stratigraphic feature model, the coordinates and geometric parameters of the large-scale fractures are inserted, directly transforming them into corresponding large-scale fracture sections, thus completing the modeling of large-scale fractures in different phases. Since small and medium-sized fractures have a relatively small spatial scale, a stochastic modeling method can be used to characterize their geometry and distribution. This method requires constraints through fracture attribute distribution patterns to randomly generate the center locations of small and medium-sized fractures, and then assigning fracture attributes to construct the corresponding fracture model. Therefore, it is necessary to integrate constraints such as fracture intensity (density), size distribution, dip / strike, and spatial clustering patterns through geostatistical understanding to generate a statistically representative fracture network. The fracture density distribution characteristics use the third-phase fracture density prediction results obtained in step 5 as attribute constraints, but the fracture density prediction results include the sum of small and medium-scale fractures, requiring the setting of the fracture distribution number according to the fracture quantity ratio. Based on statistical analysis of outcrop cracks, the ratio of small-scale to medium-scale cracks was set to 7:3 for modeling small- and medium-scale cracks. The crack orientation distribution was modeled using the Fisher model, employing the distribution range of the dip and dip angle parameters from the outcrop statistics. A power-law function was used to characterize the crack length distribution in the crack scale settings. Based on these parameter settings, models for medium-scale and small-scale cracks were obtained.

[0191] By comprehensively characterizing fracture networks at different scales and directions using multiple factors and models, a discrete fracture model at the small and medium scales, incorporating multivariate data, is established. Within the grid framework of the outcrop geological model, large, medium, and small fractures are superimposed to construct a multi-scale discrete fracture model. Following this approach, fractures at different stages can be modeled separately using geomechanical simulation results, and then superimposed to obtain the current outcrop fracture network model, such as... Figure 10 .

[0192] Step 7: Fluid simulation based on crack discrete network

[0193] Fluid simulation based on the DFN model yields the permeability distribution characteristics of the fracture network. Combining this with the permeability characteristics of the fracture network at different stages allows for further analysis of the control effect of fractures on karst development at different times. The cubic law is a fundamental model in subsurface fluid mechanics, idealizing fractures as parallel plates and deriving flow equations under the assumption of laminar flow through planar pores. The basis for permeability calculations primarily stems from the geometric characteristics of fractures and fluid flow simulations. Based on the Oda method and the cubic law, a permeability tensor is constructed. The calculations describe the laminar flow of fluid along a smooth crack system:

[0194]

[0195] in, The symbol for Kronecker; , , , represent the projections of the crack normal vector in the i and j directions, respectively. By coarsening the discrete crack network (DFN) model, key crack properties, including crack line density and connectivity, can be systematically extracted from the scaled-up mesh. Subsequently, single-phase flow simulations using the finite element method (FEM) or finite volume method (FVM) can calculate the equivalent permeability tensor in different coordinate directions, providing crucial input for studying the heterogeneity of fluid flow.

[0196] Outcrop observations show that the ancient cave system extends in the same direction as the 270° and 330° trending fractures. This geometric consistency suggests that these fractures are genetically related to the ancient karst development process. In this study, fluid flow simulations were conducted for two specific directions (270° and 330°) of the karst caves, using permeability distribution characteristics to characterize the trend direction of dissolution development. By combining the regional tectonic background and ancient karst development characteristics, two different fracture network models were used for fluid flow simulation: (1) a 330° fracture network model from the Middle Caledonian period; (2) a model from the Early Hercynian period that includes both 270° and 330° fracture networks. Each model corresponds to a specific geological period characterized by different stress states. This method can compare and analyze the fluid transport dynamics of fracture networks in different periods. The simulation results are as follows: Figure 11 .

[0197] During the late Ordovician depositional period, tectonic uplift in the southern part of the Bachu region triggered southward-to-northward paleofluid seepage activity in the Yijianfang outcrop area. Under this tectonic context, seepage simulations were performed on a model of the 330° fracture network that had developed at that time. The results show that the high permeability areas coincide with the fault and fracture development areas in the western part of the study area. Figure 11-a). The potential for karstification and cave formation is greatest along these faults, consistent with observed dissolution characteristics. In the Yifang Formation, faults typically generate a 10-20 meter dissolution control zone within which the permeability of the fault-affected fracture network is 6-16 mD. This area often develops large-scale dissolution cavities. The combined effect of the fault's influence range and hydraulic parameters determines the spatial distribution and preferred development direction of the karst system, highlighting their crucial role in the underground dissolution process. Figure 11 -c).

[0198] Analysis of the residual features of the ancient karst cave indicates that the east-west trending cave was formed by west-to-east dissolution, suggesting a genetic link to the 270° fracture. To investigate this connection, east-west seepage simulations were performed on a discrete fracture network (DFN) model incorporating 270° and 330° fracture networks. The simulation results show ( Figure 11 -b), high permeability zones exist in the southeastern and central-eastern regions with east-west trending fault systems. Under these conditions, the seepage efficiency within a fracture network is significantly better than that of a single 330° fracture network. A strong correlation exists between permeability distribution and fracture density models ( Figure 11 -d).

[0199] Step 8: Verification and correction of the prototype model

[0200] Based on the existing outcrop geological model mesh, the fracture discretization model is coarsened to obtain the fracture-matrix coupling coefficient and fracture density attributes. The fracture discretization model and matrix model can serve as prototype models for outcrop fractures. The fracture mesh and attribute model data are projected onto the existing outcrop surface. Typical outcrop profiles are selected for statistical analysis of fracture density and orientation, which is then compared and verified with the fracture density and orientation data in the prototype model. The outcrop prototype model is corrected by adjusting the boundary conditions of the geomechanical simulation to maintain basic consistency with the outcrop fracture characteristics. Comparison of the fracture orientation and length in the statistical model with the outcrop statistical results shows that the simulated fracture orientation and length characteristics are basically consistent with the outcrop statistical results. Through analysis of fracture and structural surface truncation traces, typical areas in the outcrop are selected for fracture morphology comparison. Figure 12 Comparison of -b) and crack line density ( Figure 12 -c). The comparison results show that the model's fracture morphology matches the actual fracture characteristics in the outcrop by more than 80%, and the average fracture line density similarity is above 85%. The fracture model has a high similarity to the outcrop and can be used as a prototype geological model for outcrop fractures.

[0201] Comparison of karst development zones observed at outcrops reveals a high correlation between high-value areas in fracture modeling and fluid simulation results and karst development zones. A prototype fracture model of the Yijianfang outcrop area reproduced the development characteristics of outcrop fractures, extracting fracture data from different periods and scales. Quantitative analysis of the relationship between outcrop fracture development characteristics and geological features such as distance from faults, fault strike, and tectonic curvature provides constraints for fracture-cavity evolution analysis and reservoir fracture modeling. Comparison of fracture network permeability distribution and cave development range shows that fracture development intensity dominates the spatial boundary of karst expansion, highlighting the crucial role of the fracture network in guiding fluid flow and promoting local dissolution. During the first phase of karst activity, caves mainly developed in areas with a fracture density of 0.15 fractures / m and a fracture permeability greater than 5 mD. In the second phase of karst activity, bedding-parallel karstification intensified, with faults and fractures being the main inducing factors; caves mainly developed in areas with a fracture density of 0.25 fractures / m and a fracture permeability greater than 15 mD. These data can provide crucial references for karst fissure and cavern prediction and underground reservoir evaluation.

[0202] Overall, this invention provides an effective method for constructing a prototype model of outcrop fractures in carbonate rocks, which has high practical value for predicting and modeling fractures in carbonate rocks and for digital characterizing outcrops.

[0203] The present invention has been described above by way of example, but the present invention is not limited to the specific embodiments described above. Any modifications or variations made based on the present invention shall fall within the scope of protection claimed by the present invention.

[0204] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for modeling carbonate rock fractures based on digital outcrops and geomechanical simulation, characterized in that, Includes the following steps: S1. Constructing a digital outcrop model of a typical outcrop: Based on UAV oblique photogrammetry technology, the selected outcrop area is scanned and three-dimensionally reconstructed to obtain a digital outcrop model, and geometric information of strata, faults and outcrop surface cracks is extracted from it. S2. Statistical analysis of outcrop crack parameters: Combining digital outcrop models with field identification, the orientation, length and aperture parameters of outcrop surface cracks are statistically analyzed, and different crack groups are classified according to their orientation and mechanical characteristics. S3. Determine the main fracture formation period of the outcrop crack: Determine the filling period and the formation sequence of the crack by isotope testing and structural trace analysis of the crack filling material, and obtain the direction and magnitude of paleostress in each period by combining acoustic emission experiments. S4. Paleostress field simulation during major fracture formation periods: A geomechanical model is established based on the restored paleotectonic model and rock mechanical parameters. The paleostress data obtained in S3 is used as boundary conditions to simulate the distribution of paleostress fields during each major fracture formation period. S5. Quantitative characterization of multi-stage fracture parameters: Based on the stress field results of S4, the rock fracture criterion is used to identify fracture, and the density and aperture of each stage of fracture are calculated according to the quantitative relationship model between stress-strain and fracture parameters. S6. Multi-stage, multi-scale crack discrete modeling: Deterministic modeling based on digital outcrop data is used for large-scale cracks, and stochastic modeling with crack density obtained from S5 as constraint is used for medium and small-scale cracks. Finally, the models are superimposed to form a multi-stage crack discrete network model. S7. Fluid simulation based on discrete fracture network: Fluid flow simulation is performed based on the discrete fracture network model in S6 to calculate fracture permeability and analyze the control effect of multi-stage fractures on the development of paleokarst. S8. Verification and correction of the outcrop prototype model: Compare and correct the predicted attributes of the fracture model with the measured outcrop data, and output the verified fracture-matrix coupled prototype geological model.

2. The method according to claim 1, characterized in that, In step S1, the UAV oblique photography technology includes fixed-altitude aerial photography and close-up aerial photography; wherein, the fixed-altitude aerial photography constructs a model with a resolution of 0.1~0.3 meters, which is used to analyze regional structural features; in close-up aerial photography, the distance between the UAV and the geological body is less than 15 meters, and the resolution of the constructed model is 3~6 millimeters, which is used to statistically analyze crack parameters.

3. The method according to claim 1, characterized in that, In step S2, the superposition sequence of multi-stage cracks is determined by analyzing the truncation, faulting, confinement and sliding relationships between different crack groups.

4. The method according to claim 1, characterized in that, In step S3, the acoustic emission experiment specifically involves: applying triaxial loading to directional rock plunger samples from multiple directions, identifying Kaiser effect points on the cumulative acoustic emission curve to obtain memory stress; and calculating the maximum horizontal principal stress using a formula based on the memory stress values ​​obtained from sampling tests at 0°, 45°, and 90° in the horizontal direction. Minimum horizontal principal stress and the direction angle of maximum principal stress .

5. The method according to claim 1, characterized in that, In step S5, the rock fracture criterion is as follows: when the minimum principal stress σ3 ≥ 0, the two-stage Coulomb-Mohr criterion is used to determine shear fracture, and the fracture criterion is as follows: Where σ1 is the maximum principal stress, σ3 is the minimum principal stress, C is the cohesion of the rock, and φ is the internal friction angle of the rock; when σ3<0, the Griffith criterion is used to determine tensile fracture.

6. The method according to claim 5, characterized in that, In the two-stage Coulomb-Mohr criterion, the value of the internal friction angle φ is determined according to the relationship between the confining pressure σ3 and the critical confining pressure σ0: when 0 < σ3 < σ0, φ is the first internal friction angle φ1; when σ3 ≥ σ0, φ is the second internal friction angle φ2.

7. The method according to claim 1, characterized in that, In step S5, the quantitative relationship model is constructed based on the strain energy density theory. Specifically, the fracture density is determined by calculating the ratio of the energy released by rock fracture to form the fracture surface to the energy required to generate a unit area of ​​fracture. Then, the fracture angle, the size of the characterizing unit body, and the strain state are combined to calculate the fracture linear density and aperture.

8. The method according to claim 1, characterized in that, In step S7, the calculation of fracture permeability is based on the Oda method and the cubic law, specifically the permeability tensor is calculated based on the fracture aperture, linear density and normal vector information of the fracture surface.

9. A storage device storing a computer program, characterized in that, When the computer program is executed by a processor, it implements the method as described in any one of claims 1 to 8.

10. A computing device, characterized in that, include: Memory, used to store computer programs; A processor for executing a computer program in the memory to implement the method as described in any one of claims 1 to 8.

Images