A three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic bidirectional coupling

Abstract

The present invention relates to the technical field of numerical simulation of crustal dynamics, and provides a three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic bidirectional coupling, including constructing a three-dimensional geometric model of a target area and configuring deep rock elastoplastic rheological parameters and a shallow surface erosion rate equation; generating deep unstructured grids and shallow surface irregular grids by using a differential grid strategy, and setting an adaptive time step control; configuring a nonlinear solver to achieve staged solution of thermo-mechanical coupling equations, and completing cross-scale data transfer through bilinear interpolation; running the simulation to output a material property field, a pressure field and a strain rate field; converting the results into a standardized format and visually marking stress concentration areas. Through the deep and shallow dynamic bidirectional coupling mechanism, the present invention realizes the synchronous simulation of tectonic driving and surface response, and significantly improves the simulation accuracy and engineering applicability of crustal dynamics evolution.

Description

Technical Field

[0001] The present invention relates to the technical field of numerical simulation of crustal dynamics, and provides a three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic bidirectional coupling. Background Technique

[0002] In-situ stress is the internal stress existing in natural rock masses, and its formation causes cover multiple factors such as plate extrusion, geothermal gradient change, mantle convection, gravity, and tectonic movement. Local geological features such as faults, fractures, and rock heterogeneity further lead to the dynamic differentiation and complexity of in-situ stress. The distribution of the in-situ stress field directly affects the stability of rock masses, the activation of buried faults, the evolution of oil and gas reservoirs, and the occurrence of geological disasters.

[0003] Traditional in-situ stress measurement methods, such as borehole testing and acoustic wave detection, rely on high-precision equipment, have limitations of high cost and sparse measurement points, and it is difficult to comprehensively obtain regional stress field data. Especially in areas with developed buried faults, it is difficult for traditional methods to accurately capture stress anomaly characteristics, seriously restricting the efficiency of engineering risk assessment and resource exploration. In-situ stress numerical simulation can quantitatively analyze the spatio-temporal evolution law of the stress field by constructing a three-dimensional dynamic model, reveal the interaction mechanism between deep structures and shallow processes, and provide an efficient and economical solution for buried fault location, resource target area delineation, and disaster warning. However, the existing in-situ stress numerical simulation technology still has the following key defects:

[0004] Traditional methods treat the deep structural evolution, such as mantle convection and lithospheric extension, separately from the shallow processes. For example, deep models like ASPECT focus on thermo-mechanical evolution, while shallow models like Badlands are simplified to low-dimensional simulations. There is a lack of real-time data interaction between the two, which results in the inability to feedback deep stress driving to shallow topographic evolution in real time, and the regulatory effect of shallow sediment load on the deep stress field is ignored. It is difficult to accurately simulate complex processes such as the formation of fault basins and tectonic-sedimentary co-evolution. In addition, existing technologies have not effectively solved the problem of resolution differences between the deep (kilometer-scale grid) and the shallow (hundred-meter-scale grid). Traditional interpolation algorithms have significant errors when transmitting the stress field, especially in key areas such as fault zones and topographic abrupt change regions. The grid mismatch causes the distortion of stress concentration characteristics and affects the positioning accuracy of buried faults. On the other hand, there is a time-scale difference between deep structural evolution (million-year scale) and shallow processes (thousand-year scale). Existing models mostly use fixed time steps or one-way time series synchronization, making it difficult to balance numerical stability and computational efficiency. For example, if the shallow model uses the deep adaptive time step, it is easy to cause topographic evolution oscillations due to non-satisfaction of the CFL condition; while using explicit iteration in the deep is difficult to support long-term tectonic simulations, restricting engineering applications. Finally, traditional methods rely on static parameter settings and do not fully integrate measured data such as borehole in-situ stress and LiDAR topography to dynamically optimize the model. Parameter errors such as viscosity and friction angle will accumulate over the simulation time, resulting in the deviation of prediction results from actual geological conditions, especially in high-risk areas where the errors are amplified, reducing the reliability of risk assessment.

[0005] The above problems limit the accuracy and applicability of in-situ stress numerical simulation in engineering buried fault location, resource exploration, and disaster prevention and control, and there is an urgent need for a new solution. Summary of the Invention

[0006] To solve the problems in the background technology, the present invention provides a three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic bidirectional coupling, which includes the following steps:

[0007] S1. Model configuration: Based on the geological structure characteristics of the target area, construct a three-dimensional geometric model, define the elastoplastic rheological parameters of deep rocks, establish a shallow erosion rate equation, and set the gravity field parameters, radioactive heat generation parameters, and layered temperature field configuration.

[0008] S2. Mesh setting: Use a differential mesh strategy to generate deep unstructured meshes and shallow irregular meshes, and set boundary conditions, adaptive time step control strategies, and mesh adaptive optimization strategies.

[0009] S3. Solver configuration: Configure a non-linear solver to handle the thermo-mechanical coupling equation, establish a staged solution strategy, configure a deep and shallow dynamic bidirectional coupling mechanism, and achieve cross-scale data transfer through bilinear interpolation algorithms.

[0010] S4. Run simulation: Output physical field data including viscosity, density, lithostatic pressure, tectonic stress, and strain rate, and record the topographic elevation evolution data;

[0011] S5. Data output: Generate a superimposed visualization map of the three-dimensional stress isosurface and topographic contour lines, mark the principal stress distribution and stress concentration areas, and export the key parameter table;

[0012] S6. Post-processing: Calculate the error between the simulation results and geological observation data, and trigger the parameter optimization loop when the error exceeds the threshold; Use a filtering algorithm to correct the viscosity, friction coefficient, and calibrate the erosion equation coefficient; Generate a fault location and activation risk assessment report.

[0013] Furthermore, in S1, the construction of the three-dimensional geometric model includes: Based on the actual geological background and geographical environment of the target area, define the length L, width W, and depth H of the model as the three-dimensional spatial range of the target research area from the deep crust to the surface; Among them, in the elastoplastic rheological parameters of deep rocks, the unit of viscosity η is Pa·s, which characterizes the flow resistance of rocks under high stress; the unit of the internal friction angle φ is degree, which is determined by rock mechanics experiments and is used to describe the non-linear relationship between the shear strength of rocks and the normal stress; the unit of cohesion C is MPa, which characterizes the binding force between rock particles and is calibrated by triaxial compression experiments; the erosion rate equation of the shallow geodynamic model is , where Q is the surface water flow rate (m³ / s), S is the slope, k is an empirical coefficient ( ), x and y are exponential parameters fitted based on regional hydrological observation data; the gravitational acceleration g is set to 9.81 m / s², and the direction is vertically downward; the unit of the radioactive heat generation rate A is μW / m³, which is calculated and determined according to the abundance of radioactive elements and the decay heat generation rate in the target area; among them:

[0014] The elastoplastic rheological parameters adopt the Drucker-Prager model, and its yield criterion is:

[0015]

[0016] Among them, is the second invariant of the deviatoric stress, is the first invariant of the stress tensor, , ; is the internal friction angle, which is an angular parameter characterizing the shear strength of rocks and reflects the internal friction characteristics of materials during shear failure; this model introduces a temperature-dependent viscosity η(T) to characterize the creep behavior of deep rocks in a high-temperature environment, and the viscosity expression is:

[0017]

[0018] Among them, is the reference viscosity, is the activation energy, is the gas constant, is the absolute temperature (K);

[0019] Furthermore, in S1, the temperature field configuration includes: the upper crust temperature is jointly determined by heat flow and radioactive heat generation, and the calculation formula is: where T is the absolute temperature (K), z is the depth (km), 0.065 is the heat flow density (W / m²), 2.5 is the thermal conductivity of the rock (W / (m·K)), 1.5× is the radioactive heat generation rate (μW / m³); the lower crust temperature linearly decreases to 1573 K at the bottom of the lithosphere.

[0020] Furthermore, in S2, the generation of the deep unstructured grid includes: initially discretizing the deep model using tetrahedral elements, and the element size is dynamically adjusted according to the fracture zone density. The grid resolution in the fracture zone area is encrypted to the 100-meter level; the shallow irregular grid is generated by the Delaunay triangulation algorithm, the horizontal resolution is set to 50 - 100 meters, and according to the terrain curvature, the curvature threshold is set to to trigger local encryption to ensure grid refinement in the terrain mutation area; the horizontal tectonic movement rate of the lateral boundary of the deep model is set according to the regional plate movement data, with a range of 0.1 - 10 mm / yr to simulate the crustal extension or extrusion process.

[0021] Furthermore, in S2, the adaptive time step control strategy includes: the deep dynamic step size is adaptively adjusted according to the stress change rate When MPa / yr, the time step is shortened to within 100 years; when MPa / yr, the time step is extended to 1000 years; the shallow fixed time step follows the CFL condition, that is:

[0022]

[0023] where is the minimum grid size (m), is the maximum surface water flow velocity (m / s) to ensure the numerical stability of terrain evolution;

[0024] The grid adaptive optimization strategy includes: when the deep local strain rate exceeds , grid encryption is triggered, and the grid resolution in the fracture zone area is increased to 200 meters; the shallow model triggers local refinement of the Delaunay triangulation according to the terrain curvature threshold to ensure that the grid size in the area where the slope change rate ≥ 5° / 100 meters is ≤ 50 meters.

[0025] Furthermore, in S3, the staged solution strategy includes: in the first stage, the algebraic multigrid (AMG) preconditioning technique is adopted to reduce the residual of the Stokes equation to 10% of the initial value; in the second stage, the Newton-Krylov iteration method is used with a tolerance of 1× , and the maximum number of iterations is 5000; for cross-scale data transfer, the bilinear interpolation algorithm is used to map the stress field of the deep kilometer-scale grid to the shallow hundred-meter-scale grid. The interpolation weight is calculated according to the element shape function, and the error is controlled within 5%.

[0026] In S3, the process of configuring the shallow-deep coupling mechanism includes: the shallow sediment thickness d sed dynamically feedbacks the deep viscosity field through the ballast effect, and the viscosity correction formula is:

[0027]

[0028] where d sed is in meters, is the initial viscosity (Pa·s); the deep horizontal tectonic stress σ xy is transferred to the shallow model through bilinear interpolation to drive the topographic elevation h to calculate the settlement amount according to the formula Δh = σ xy / (ρg), where ρ is the rock density (kg / m³) and g is the acceleration due to gravity (m / s²).

[0029] Furthermore, in S4, the specific process of running the simulation includes:

[0030] The lithostatic pressure is calculated by the formula , where:

[0031] is the rock density, characterizing the mass distribution characteristics of the rock mass; is the acceleration due to gravity, with the direction vertically downward; is the current depth, representing the vertical distance from the target point to the surface;

[0032] The tectonic stress , representing the dynamic stress field generated by crustal tectonic movements, is obtained by solving the Stokes equation and is superimposed with the lithostatic pressure to obtain the total pressure field ;

[0033] The displacement of the buried fault is calculated by the time integration of the velocity field , that is , with an accuracy of centimeter level; where: is the starting time of the simulation; is the current simulation time.

[0034] Further, in S5, the generation of the three-dimensional stress isosurface includes: the maximum principal stress and the minimum principal stress obtained by eigenvalue decomposition, and the area with the principal stress ratio is marked as the stress concentration area; the rockburst risk area is determined by the Griffith criterion, that is:

[0035]

[0036] wherein, is the tensile strength of the rock (MPa), and the threshold is set to 5 MPa; the fault slip amount table records the displacement , the slip rate and the activation probability , wherein is calculated based on the ratio of the shear stress to the normal stress .

[0037] Further, in S6, the parameter optimization loop includes: using the Kalman filter algorithm to correct the viscosity and the internal friction angle , and the observed data is the borehole in-situ stress measurement value; the erosion equation coefficients , , are calibrated by the Markov chain Monte Carlo MCMC method, the objective function is the root mean square error RMSE between the simulated terrain elevation and the LiDAR measured data, and the threshold is set to ±2 m; when the error exceeds the limit, 10 iterations of optimization are triggered until RMSE < 1 m.

[0038] In S6, the fault location report includes determining the spatial coordinates, strike angle , dip angle and slip type of the hidden fault; the activation risk assessment is based on the probability distribution of the historical earthquake catalog and the simulated slip rate, and the risk level is divided into low, medium, and high, and a risk heat map and emergency plan suggestions are output.

[0039] The present invention also proposes the use of the three-dimensional in-situ stress numerical simulation method based on the deep and shallow dynamic bidirectional coupling in the hidden fault location and resource exploration.

[0040] The beneficial effects achieved by the present invention are as follows:

[0041] First, the present invention integrates the deep dynamics module ASPECT and the shallow geodynamics module Badlands, and constructs a two-way dynamic interaction mechanism for the deep and shallow models. In each time step, the deep stress field drives the evolution of the shallow surface topography in real time, while the shallow sediment load is feedback to the deep viscosity field through the ballast effect, forming a closed-loop coupling, breaking through the limitation of the traditional model in dealing with the deep and shallow processes in isolation, realizing the synchronous simulation of tectonic driving and surface response, significantly improving the physical consistency of the crustal dynamics evolution, and providing a more realistic dynamic environment for the simulation of complex processes such as the activation of buried faults and the formation of basins.

[0042] Second, the present invention solves the difference in the grid resolution of the deep and shallow models, uses the bilinear interpolation algorithm to realize the cross-scale stress field mapping, and optimizes the key areas through the dynamic encryption strategy of the unstructured grid. The grid of the deep fracture zone is encrypted to the 200-meter level according to the strain rate threshold, and the local grid is refined to 50 meters in the area with sudden change of the shallow surface topography curvature, effectively reducing the numerical error caused by the grid difference, while improving the ability to capture details such as stress concentration in the fracture zone and terrain mutation, and enhancing the analysis accuracy of the model for complex geological structures.

[0043] Third, in view of the difference in the time scales of the deep and shallow processes, the present invention designs a differential time step control strategy. The deep part adopts a dynamic step size (100 - 1000 years) adaptive to the stress change rate, and the shallow part sets a fixed step size based on the CFL condition to ensure the numerical stability of the topography evolution. At the same time, the solver adopts a staged processing, preprocessing + high-precision iteration and hybrid iteration algorithms to solve the nonlinear problem of the thermo-mechanical coupling equation, taking into account the simulation requirements of the long-term tectonic evolution in the deep part and the instantaneous process in the shallow part. On the premise of ensuring the accuracy, the calculation efficiency is greatly improved, supporting the engineering application of large-scale three-dimensional simulation.

[0044] Fourth, the present invention introduces the Kalman filter and the Markov chain Monte Carlo MCMC method to realize the dynamic optimization of the model parameters. When the simulation error exceeds the threshold, an iterative calibration loop is triggered. Combining the measured data such as borehole in-situ stress and LiDAR topography, the rheological parameters and boundary conditions are corrected. Through data-driven optimization, the root mean square error RMSE between the simulation result and the geological observation is reduced to less than 1 meter, significantly improving the adaptability of the model to the actual geological environment, and providing high-confidence data support for the location of buried faults and the assessment of rockburst risk.

[0045] Fifth, the present invention realizes multi-dimensional analysis of data such as the in-situ stress field and topographic evolution. The three-dimensional stress isosurface is superimposed on the topographic contour lines to mark the stress concentration areas and rock burst risk areas, and a parametric report such as the fault slip amount and activation probability is generated, converting the abstract numerical results into intuitive engineering decision-making basis, greatly improving the positioning accuracy of hidden faults. At the same time, it supports actual scenarios such as resource exploration target area delineation and engineering disaster warning, reducing the exploration cost. Description of the Drawings

[0046] Figure 1 is a flow chart of a three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic bidirectional coupling;

[0047] Figure 2 is the initial three-dimensional mesh division diagram in Embodiment 2;

[0048] Figure 3 is the three-dimensional initial geothermal cloud map of the calculation area in Embodiment 3.

[0049] Figure 4 is the three-dimensional in-situ stress and strain cloud map of a certain project (oblique) longitudinal section in Embodiment 4.

[0050] Figure 5 The in-situ stress and strain diagram of the calculation area after adaptive mesh refinement near the fault in Embodiment 5.

[0051] Figure 6 is the detailed flow chart of the three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic bidirectional coupling of the present invention. Detailed Embodiments

[0052] Next, the technical solutions in the present invention will be clearly and completely described in conjunction with the accompanying drawings in the present invention. In addition, the forms of the structures described in the following embodiments are only examples, and the present invention is not limited to the structures described in the following embodiments. All other embodiments obtained by those of ordinary skill in the art without creative efforts belong to the scope of protection of the present invention.

[0053] The parameter symbols, units and explanations involved in the present invention:

[0054] L, W, H (km or m): Model space parameters, representing the length, width and depth of the three-dimensional geometric model of the target area respectively;

[0055] η (Pa·s): Viscosity of deep rock, characterizing the rock flow resistance, varying with temperature (formula: η(T) = ·exp( / (RT)), where is the reference viscosity, is the activation energy, R is the gas constant, and T is the absolute temperature);

[0056] φ (°): The angle of internal friction, determined through rock mechanics experiments, describes the non-linear relationship between the shear strength of rock and the normal stress;

[0057] C (MPa): Cohesion, characterizing the bonding force between rock particles, calibrated by triaxial compression experiments;

[0058] : Superficial erosion rate equation, Q is the surface water flow rate (m³ / s), S is the slope (dimensionless), k is an empirical coefficient ( ), x and y are exponential parameters fitted based on hydrological data;

[0059] g (9.81 m / s²): Acceleration due to gravity, with the direction vertically downward;

[0060] A (μW / m³): Radioactive heat generation rate, determined by calculating the abundances of radioactive elements and the heat generation rate of decay in the target area;

[0061] (dimensionless): The second invariant of the deviatoric stress, characterizing the shear stress component;

[0062] (MPa): The first invariant of the stress tensor, i.e., the hydrostatic pressure;

[0063] α, β (dimensionless): Drucker-Prager yield criterion parameters, α = 2sinφ / (√3(3 - sinφ)), β = 6Ccosφ / (√3(3 - sinφ));

[0064] T (K): Absolute temperature, the temperature of the upper crust is jointly determined by heat flow and radioactive heat generation (formula: T = 273 + (0.065 / 2.5)(20 - z) - (1.5× (20 - z)²) / (2×2.5)), z is the depth (km), 0.065 W / m² is the heat flux density, and 2.5 W / (m·K) is the thermal conductivity of the rock;

[0065] Δσ / Δt (MPa / yr): Stress change rate, used for controlling the dynamic time step of the deep model (when Δσ / Δt > 1 MPa / yr, the time step is shortened to 100 years, and when < 0.1 MPa / yr, it is extended to 1000 years);

[0066] Δx (m): The minimum grid size of the superficial model, combined with the CFL condition (Δt ≤ Δx / v max ) to set a fixed time step, v max is the maximum surface water flow velocity (m / s);

[0067] : Strain rate threshold (triggers grid encryption to 200 m at depth );

[0068] Curvature ( ): Superficial terrain curvature, threshold 0.01 triggers local grid refinement to 50 m;

[0069] , (MPa): Maximum and minimum principal stresses, obtained through eigenvalue decomposition, the area where the principal stress ratio / > 3 is the stress concentration area;

[0070] (5 MPa): Rock tensile strength threshold, the condition for rockburst risk determination is - > 2T0;

[0071] D (mm): Displacement of the buried fault, calculated by integrating the change of velocity field v over time (D = ∫v dt);

[0072] v slip (mm / yr): Fault slip rate;

[0073] P activate (%): Fault activation probability, calculated based on the ratio of shear stress τ to normal stress σ n ;

[0074] RMSE (m): Root mean square error, used to calibrate the coefficients of the erosion equation (threshold ±2 m, target RMSE <1 m after optimization);

[0075] θ (°): Fault strike angle, 0° is due north;

[0076] δ (°): Fault dip angle, range 0° - 90°;

[0077] ρ (kg / m³): Rock density, used for calculating the lithostatic pressure;

[0078] k, x, y: Coefficients of the erosion equation, calibrated by the MCMC method, the objective function is to match the LiDAR terrain data;

[0079] d sed (m): Superficial sediment thickness, used to correct the deep viscosity through the ballast effect feedback (formula: η new = ·(1 + 0.12d sed ));

[0080] σxy (MPa): Deep horizontal tectonic stress, driving the topographic settlement Δh = σ xy / (ρg).

[0081] English abbreviations and their definitions involved in the present invention:

[0082] AMG (Algebraic Multigrid): Algebraic multigrid preconditioning technology, used to reduce the stiffness of the Stokes equations;

[0083] CFL (Courant-Friedrichs-Lewy): Numerical stability criterion, restricting the time step of the shallow model;

[0084] MCMC (Markov Chain Monte Carlo): Markov chain Monte Carlo method, used to calibrate the coefficients of the erosion equation;

[0085] LiDAR (Light Detection and Ranging): Light detection and ranging technology, providing high-precision topographic observation data;

[0086] HDF5 (Hierarchical Data Format 5): Scientific data storage format, supporting the reading and writing of large-scale physical field data;

[0087] VTK (Visualization Toolkit): 3D visualization toolkit, used to generate the overlay map of stress isosurfaces and topography;

[0088] DEM (Digital Elevation Model): Digital elevation model, used for terrain evolution verification;

[0089] RMSE (Root Mean Square Error): Root mean square error, measuring the deviation between the simulated terrain and the measured data;

[0090] GIS (Geographic Information System): Geographic information system, generating risk heat maps and spatial analysis;

[0091] PETSc (Portable, Extensible Toolkit for Scientific Computation): Parallel computing library, supporting the solution of the deep-shallow coupling model;

[0092] NetCDF (Network Common Data Form): Scientific data interface standard, used for cross-platform data interaction;

[0093] Drucker - Prager: An elastoplastic mechanics model, a criterion for describing the yield behavior of rocks;

[0094] Boussinesq approximation: An assumption for simplifying density changes in thermal convection simulations;

[0095] Paraview: An open - source 3D visualization software for rendering stress fields and temperature fields;

[0096] GMT (Generic Mapping Tools): A geographical mapping tool for generating two - dimensional dynamic evolution maps.

[0097] Refer to Figure 1 - Figure 6 , the present invention provides a three - dimensional in - situ stress numerical simulation method based on deep - shallow dynamic bidirectional coupling, which includes the following steps:

[0098] S1. Model configuration: Construct a three - dimensional geometric model based on the geological structure characteristics of the target area, and set the spatial parameters of the model including length L, width W, and depth H; Define the elastoplastic rheological parameters of deep rocks, including viscosity η, a physical quantity characterizing the flow resistance of rocks, internal friction angle φ, an angular parameter describing the shear strength of rocks, and cohesion C, a stress parameter characterizing the bonding force between rock particles; Establish an erosion rate equation for the shallow geodynamic model, and the equation is a function including flow rate Q, slope S, and empirical coefficients k, x, y; Set the gravity field parameters, including gravitational acceleration g, with a direction perpendicular to the downward direction, and the heat generation rate A of the radioactive heat - generating layer, characterizing the heat generation power of radioactive elements;

[0099] S2. Mesh setting: Generate deep unstructured meshes using tetrahedral elements and shallow irregular meshes using triangular elements; Set the boundary conditions of the deep model: The bottom is a free - slip boundary allowing tangential movement and restricting normal displacement, and a horizontal tectonic movement rate is applied laterally; Set the open - boundary conditions of the shallow model to allow sediment migration and the bottom coupling interface to receive deep - stress drive; Set a differential configuration time - step control strategy: The deep part adopts a dynamic step size adaptive to the stress change rate, and the shallow part adopts a fixed step size based on the CFL condition stability criterion;

[0100] S3. Solver configuration: Configure a non - linear solver to handle the thermo - mechanical coupling equation, and set the iteration tolerance and the maximum number of iterations; Adopt a phased solution strategy to handle the Stokes equation: In the first stage, reduce the equation stiffness through pre - processing, and in the second stage, perform high - precision iterative correction; Implement cross - scale data transfer between deep and shallow meshes through an interpolation algorithm, including stress - field mapping and terrain feedback;

[0101] S4. Run Simulation: Based on the deep-shallow dynamic bidirectional coupling mechanism, solve the deep thermal-mechanical equations and the shallow erosion dynamics equations in real-time and synchronously, and output the fully coupled physical field data including the viscosity field (dynamically updated according to the temperature-dependent viscosity formula in S1), density field, lithostatic pressure field (calculated by the formula in S1), tectonic stress field (obtained by solving the Stokes equation with the non-linear solver configured in S3), and strain rate field; synchronously record the topographic elevation evolution data (topographic changes driven by the shallow erosion rate equation in S2), and calculate the displacement of the buried fault using the implicit time integration algorithm with centimeter-level accuracy; perform bidirectional data synchronization every 100 shallow time steps and 10 deep time steps, and implement cross-scale stress field mapping through the bilinear interpolation algorithm in S3 to ensure the physical consistency between the deep tectonic stress and the shallow topographic evolution;

[0102] S5. Data Output: Convert the simulation results into a standardized format file, including the overlay visualization of the three-dimensional stress isosurfaces and topographic contour lines, is the maximum principal stress, is the minimum principal stress; based on the principal stress ratio determination, mark the stress concentration areas and rockburst risk areas; export the key parameter tables of the fault slip amount and sedimentation rate;

[0103] S6. Post-processing and Verification: Calculate the error between the simulation results and the geological observation data, and trigger the parameter optimization loop when the error exceeds the threshold; use the filtering algorithm to correct the viscosity, friction coefficient, and calibrate the erosion equation coefficient; generate a fault location and activation risk assessment report.

[0104] In S1, the construction of the three-dimensional geometric model further includes: Based on the actual geological background and geographical environment of the target area, define the length L, width W, and depth H of the model as the three-dimensional spatial range of the target research area from the deep crust to the surface; among them, in the elastoplastic rheological parameters of the deep rock, the unit of viscosity η is Pa·s, which characterizes the flow resistance of the rock under high stress; the unit of the internal friction angle φ is degree, which is determined by rock mechanics experiments and is used to describe the non-linear relationship between the shear strength of the rock and the normal stress; the unit of cohesion C is MPa, which characterizes the bonding force between rock particles and is calibrated by triaxial compression experiments; the erosion rate equation of the shallow geodynamics model is , where Q is the surface water flow rate (m³ / s), S is the slope, k is the empirical coefficient ( ·s²), x and y are the exponential parameters fitted based on the regional hydrological observation data; the gravitational acceleration g is set to 9.81 m / s², and the direction is vertically downward; the unit of the radioactive heat generation rate A is μW / m³, which is calculated and determined according to the abundance of radioactive elements and the decay heat generation rate in the target area.

[0105] In S1, the elastoplastic rheological parameters further adopt the Drucker-Prager model, and its yield criterion is:

[0106]

[0107] where is the second invariant of the deviatoric stress, is the first invariant of the stress tensor, , ; is the internal friction angle, an angular parameter characterizing the shear strength of the rock, reflecting the internal friction characteristics of the material during shear failure; by introducing the temperature-dependent viscosity η(T), this model characterizes the creep behavior of deep rocks in a high-temperature environment, and the viscosity expression is:

[0108]

[0109] where is the reference viscosity, is the activation energy, is the gas constant, is the absolute temperature (K);

[0110] In S1, the temperature field configuration further includes: The temperature of the upper crust is jointly determined by the heat flow and radioactive heat generation, and the calculation formula is: where T is the absolute temperature (K), z is the depth (km), 0.065 is the heat flow density (W / m²), 2.5 is the thermal conductivity of the rock (W / (m·K)), 1.5× is the radioactive heat generation rate (μW / m³); The temperature of the lower crust linearly decreases to 1573 K at the bottom of the mantle.

[0111] In S2, the generation of the deep unstructured grid further includes: The deep model is initially discretized using tetrahedral elements, and the element size is dynamically adjusted according to the fracture zone density. The grid resolution in the fracture zone area is encrypted to the hundred-meter level; The shallow irregular grid is generated by the Delaunay triangulation algorithm, the horizontal resolution is set to 50 - 100 meters, and according to the terrain curvature, the curvature threshold is set to 0.01 , triggering local encryption to ensure grid refinement in the terrain mutation area; The horizontal tectonic movement rate of the lateral boundary of the deep model is set according to the regional plate movement data, ranging from 0.1 - 10 mm / yr, simulating the crustal extension or extrusion process.

[0112] In S2, the time step control strategy further includes: The deep dynamic step size is adaptively adjusted according to the stress change rate , when MPa / yr, the time step is shortened to within 100 years; when When the step size is extended to 1000 years at MPa / yr; the shallow fixed step size follows the CFL condition, that is:

[0113]

[0114] Among them, is the minimum grid size (m), is the maximum surface water flow velocity (m / s), ensuring the numerical stability of terrain evolution;

[0115] In S2, the grid adaptive optimization strategy further includes: when the deep local strain rate exceeds grid encryption is triggered, and the grid resolution in the fracture zone area is increased to 200 meters; the shallow model triggers local refinement of Delaunay triangulation according to the topographic curvature threshold of 0.01 to ensure that the grid size in the area where the slope change rate ≥ 5° / 100 m is ≤ 50 m.

[0116] In S3, the staged solution strategy further includes: in the first stage, the algebraic multigrid AMG preconditioning technology is used to reduce the residual of the Stokes equation to 10% of the initial value; in the second stage, the Newton-Krylov iteration method is used, with a tolerance of 1× , and the maximum number of iterations is 5000 times; the cross-scale data transfer uses the bilinear interpolation algorithm to map the stress field of the deep kilometer-level grid to the shallow hundred-meter-level grid, and the interpolation weight is calculated according to the element shape function, with the error controlled within 5%;

[0117] In S3, the realization of the deep-shallow coupling mechanism further includes: the shallow sediment thickness d sed dynamically feedbacks the deep viscosity field through the ballast effect, and the viscosity correction formula is:

[0118]

[0119] Among them, d sed is in meters, is the initial viscosity (Pa·s); the deep horizontal tectonic stress σ xy is transferred to the shallow model through bilinear interpolation, and drives the topographic elevation h to calculate the settlement amount according to the formula Δh = σ xy / (ρg), where ρ is the rock density (kg / m³) and g is the acceleration due to gravity (m / s²).

[0120] In S4, the lithostatic pressure is calculated by the formula , where:

[0121] is the rock density, characterizing the mass distribution characteristics of the rock mass; is the acceleration due to gravity, with the direction perpendicular to the downward direction; is the current depth, representing the vertical distance from the target point to the ground surface;

[0122] Tectonic stress , representing the dynamic stress field generated by crustal tectonic movements, is obtained by solving the Stokes equation and is superimposed with the lithostatic pressure to obtain the total pressure field ;

[0123] Displacement of the buried fault is calculated by time integration of the velocity field , that is , with an accuracy of centimeter level; where: is the starting time of the simulation; is the current simulation time.

[0124] In S5, the generation of the three-dimensional stress isosurface further includes: the maximum principal stress and the minimum principal stress are obtained through eigenvalue decomposition, and the area with the principal stress ratio is marked as the stress concentration area; the rockburst risk area is determined by the Griffith criterion, that is:

[0125]

[0126] where is the tensile strength of the rock (MPa), and the threshold is set to 5 MPa; the displacement of the fault slip is recorded in the table , the slip rate and the activation probability , where is calculated based on the ratio of the shear stress to the normal stress .

[0127] In S6, the Kalman filter algorithm is used to correct the viscosity and the internal friction angle , and the observed data is the measured value of the in-situ stress of the borehole; the coefficients of the erosion equation , , are calibrated by the Markov chain Monte Carlo MCMC method, and the objective function is the root mean square error RMSE between the simulated terrain elevation and the LiDAR measured data, with the threshold set to ±2 m; when the error exceeds the limit, 10 iterations of optimization are triggered until RMSE < 1 m. The fault location report further includes: determining the spatial coordinates of the buried fault, the strike angle , the dip angle and slip types; the activation risk assessment is based on the probability distribution of historical earthquake catalogs and simulated slip rates. The risk levels are divided into low, medium, and high, and a risk heat map and emergency plan suggestions are output.

[0128] In S6, η (viscosity): a measure of the fluid flow resistance, with the unit of Pascal-second (Pa·s). It reflects the flow characteristics of rock formations or magma in the geological model. φ (internal friction angle): an index of the shear strength of materials, with the unit of degree (°). It is used to describe the stability of faults or rock formations under shear stress. Observation data: borehole in-situ stress measurement values, including vertical stress (σ v ) and horizontal stress (σ h ), with the unit of megapascal (MPa). It is directly obtained through borehole sensors and is used to invert formation mechanical parameters.

[0129] The Kalman filter process includes: prediction, predicting the current state of η and φ based on the geodynamic model. Update, comparing the measured borehole in-situ stress values with the predicted values and adjusting the parameters by minimizing the covariance matrix. Iteration, dynamically correcting the parameters to adapt to the formation heterogeneity.

[0130] Calibrate the erosion equation coefficients (k, x, y) by MCMC. The erosion equation is in the form of , where:

[0131] : erosion rate (m / yr); : runoff (m³ / s); : slope (°); : empirical coefficients to be calibrated. Objective function: root mean square error (RMSE);

[0132] Formula: simulate , with the unit of meter (m).

[0133] The goal is to make the RMSE between the simulated terrain elevation and the LiDAR measured data < 1 meter.

[0134] The MCMC method uses the Metropolis-Hastings algorithm to generate the posterior distribution of parameters and explores high-probability parameter combinations through the acceptance-rejection criterion. If RMSE > 2 meters, start 10 iterations of optimization. Each iteration updates k, x, y and re-evaluates RMSE until the standard is met.

[0135] The process of fault location and activation risk assessment includes determining the parameters of buried faults.

[0136] Spatial coordinates: The three-dimensional position of the fault, including longitude (λ), latitude (φ), and depth (z, unit: km). Strike angle θ: The extension direction of the fault line on the horizontal plane, with due north being 0°, measured clockwise (range 0 - 360°). Dip angle δ: The angle between the fault plane and the horizontal plane (range 0 - 90°).

[0137] Slip types include strike-slip: horizontal displacement of plates (such as the San Andreas Fault). Thrust: vertical compression of plates (such as the Himalayan Fault). Normal slip: vertical stretching of plates (such as the East African Rift).

[0138] Activation risk assessment process: Data sources include historical earthquake catalogs: recording the time, location, magnitude (Mw), and focal mechanism of historical earthquakes. Simulated slip rate: Calculating the annual slip of the fault (unit: mm / yr) through a dynamic model to generate a probability density function (PDF).

[0139] The risk level classification process is as follows: Low risk: The peak of the slip rate PDF < 1 mm / yr, and the recurrence period of historical earthquakes > 1000 years. Medium risk: Slip rate 1 - 5 mm / yr, recurrence period 100 - 1000 years. High risk: Slip rate > 5 mm / yr, recurrence period < 100 years.

[0140] Output results include: Risk heat map: The spatial distribution of risk levels represented by color gradients (green - yellow - red) in a GIS map. Emergency plan suggestions: High-risk areas: Deploy real-time monitoring, formulate evacuation routes, and reinforce infrastructure. Medium-risk areas: Conduct regular earthquake drills and improve building seismic standards. Low-risk areas: Public science popularization and long-term monitoring plans.

[0141] Example 1, the specific process of a three-dimensional in-situ stress numerical simulation method based on deep-shallow dynamic bidirectional coupling is described through the specific steps of this example.

[0142] For the high-concurrency requirements of deep - shallow coupling simulation, in this example, before starting the operation module, we need to configure high-performance computing resources and software environments to ensure the efficient operation of large-scale parallel tasks.

[0143] Submit parallel task scripts through the job scheduling system to apply for multiple computing nodes, fully meeting the high concurrency requirements of three-dimensional crustal dynamics simulation, saving time, and improving the computing speed and efficiency. To support the parallel computing requirements of large-scale three-dimensional simulations, the PETSc parallel computing library is loaded synchronously to achieve efficient scheduling of multi-node tasks; at the same time, the HDF5 and NetCDF data interface libraries are integrated to ensure the standardized reading and writing of deep stress fields, shallow surface topography, and coupled intermediate files. The coupling controller is embedded in the system as the core hub, responsible for managing data exchange and timing synchronization between the deep and shallow modules. In addition, the input / output paths are standardized to simplify the data management process and avoid file conflicts. At the same time, a dedicated Python environment is activated, and a geographic information processing toolchain is integrated for spatial interpolation, dynamic visualization of terrain data, and generation of result maps.

[0144] Through the above steps, the system has completed the full-process initialization from hardware resource allocation to software function integration, laying a reliable technical foundation for subsequent deep-shallow coupling simulation.

[0145] Construct a three-dimensional geometric model based on the geological background, define material properties and physical boundary conditions to ensure the accuracy of the simulated physical processes.

[0146] First, we construct the geometric model of this area based on the actual regional geological background and geographical environment of the basin. For the deep model, according to the actual natural geographical environment of the basin, the length, width, and depth are set to construct a rectangular area and perform initial meshing on this area. The initial mesh uses unstructured tetrahedral elements, and then the adaptive mesh function is loaded on the fault zone area of this structural model to simulate the tectonic activities in the deep crust. The shallow surface dynamics model coupled with the deep model generates irregular triangular meshes (with a horizontal resolution reaching the order of 100 meters) based on the upper surface of the deep model and uses the Delaunay triangulation algorithm to optimize the mesh quality. This is because irregular triangular meshes fit the geographical state of the shallow surface topography better, thereby improving the accuracy of model fitting.

[0147] Material properties: To accurately simulate the deep-shallow coupling dynamics process of the crust, the material properties and boundary conditions are refined. The deep rock adopts the Drucker-Prager elastoplastic rheological model, and its parameters are based on rock mechanics experiments and geological exploration data. This model can accurately characterize the fracture behavior and long-term creep characteristics of rocks under high stress. Combined with temperature-dependent viscosity parameters, it effectively simulates the driving effect of mantle convection on deep tectonic activities.

[0148] The dynamic behavior of shallow surface sediments is described by the erosion rate equation ( )(Quantification is carried out, where Q is the flow rate and S is the slope. This equation can reflect the dynamic feedback of surface material migration on tectonic activities by integrating hydrological observation data and terrain evolution laws. The erosion rate in areas with high flow rates and steep slopes increases significantly, resulting in rapid topographic incision and affecting the local stress distribution.)

[0149] Boundary conditions: The boundary conditions are set considering both physical rationality and numerical stability:

[0150] Deep model: A free-slip boundary is adopted at the bottom, and horizontal stretching rates are applied on the east and west sides to simulate the accumulation of tectonic stress during the separation of terranes; the north and south sides are fixed to restrict the lateral expansion of the simulation area.)

[0151] Shallow model: The bottom is coupled with the deep model, receiving deep dynamic driving through stress field transfer, and at the same time feeding back sediment loads to the deep; the open boundary allows sediment input and output to simulate the long-term impact of surface material migration on terrain evolution.)

[0152] The velocity field is defined by a piecewise function, and a reverse stretching rate is applied in a specific spatial area to simulate the differential movement of local fault zones.)

[0153] Temperature field configuration: The temperature field configuration is set by stratification according to the heat flow distribution and radioactive heat generation characteristics:

[0154] The temperature in the upper crust (0 - 10 km) is jointly determined by heat flow and radioactive heat generation, and the calculation formula is:

[0155]

[0156] The formula starts from the reference temperature of 273 K (0 °C) and superimposes two parts of effects:

[0157] Heat flow contribution: , where 0.065 represents the heat flow density (unit: W / m²), 2.5 is the thermal conductivity (unit: W / (m·K)), and z is the current depth (kilometer). This part indicates that as the depth increases (z increases), the temperature increment caused by the downward transfer of heat flow gradually decreases.)

[0158] Radioactive heat generation correction: , which reflects the non-linear effect of the heat generated by the decay of radioactive elements on temperature through a quadratic term, and the negative sign indicates the attenuation of the heat generation effect with increasing depth.)

[0159] This formula comprehensively considers the heat conduction and heat generation effects of radioactive elements, and accurately depicts the thermodynamic state of the strata.)

[0160] The temperature in the lower crust (10 - 40 km) decreases linearly to 1573 K at the bottom of the lithosphere, reflecting the thermal gradient change during the cooling of the lithosphere.)

[0161] Implementation of the deep - shallow coupling mechanism: In the deep - shallow coupling dynamics framework, the cyclic co - evolution of deep and shallow processes is realized through a dynamic two - way interaction mechanism. The horizontal stress calculated by the deep module transfers parameter changes to the shallow module in real - time through the coupling controller, driving topographic evolution. For example, when the tensile stress at a certain deep location exceeds the threshold, the stress gradient will trigger basin fault - subsidence, and the shallow module adjusts the topographic elevation h and sediment distribution d accordingly. sed Meanwhile, the sediment load simulated by the shallow module is dynamically fed back to the deep viscosity field through the ballast effect, forming a closed - loop coupling. Specifically, an increase in sediment thickness will cause changes in the deep viscosity, which in turn affects lithospheric deformation.

[0162] Mesh generation and adaptive optimization: To achieve a balance between computational efficiency and accuracy, this scheme adopts a phased mesh construction and dynamic optimization strategy. The initial mesh is constructed based on geological structure characteristics and simulation requirements:

[0163] The deep model uses unstructured tetrahedral meshes and sets the initial element size. This design can effectively reduce computational redundancy and provide a basis for subsequent adaptive refinement. For key structural regions such as fault zones and basin centers, grid encryption is preferentially implemented to improve the resolution for accurately capturing complex dynamic behaviors such as stress concentration, rock deformation, and heat convection.

[0164] The shallow model generates irregular triangular meshes based on the upper surface of the deep model and sets the initial horizontal resolution. Such a mesh structure can better fit the complex surface topography (including mountains and valleys), and local refinement is performed on regions with significant topographic curvature through a dynamic encryption strategy to improve the resolution and ensure high - precision characterization of micro - geomorphic features such as erosion gullies and alluvial fans.

[0165] The adaptive optimization strategy resolves the contradiction between limited computational resources and the multi - scale characteristics of physical processes by dynamically adjusting the mesh resolution:

[0166] Trigger conditions: In the deep model, when the local strain rate exceeds the threshold, it indicates that there is significant rock plastic deformation or stress accumulation in this region, triggering dynamic grid encryption to improve numerical accuracy; in the shallow model, according to sudden changes in topographic curvature or sediment thickness (including landslide and alluvial fan development areas), the grid resolution is automatically increased to capture the details of instantaneous topographic evolution.

[0167] Parameter configuration: Each time an adjustment is made, high - strain - rate elements are refined to focus on key regions, while low - change regions are coarsened to reduce redundant calculations. The optimization operation is performed every 10 time steps, which can not only avoid performance losses caused by frequent adjustments but also ensure the reasonable allocation of computational resources.

[0168] Through phased construction and dynamic optimization, while ensuring the accuracy of high-concern areas such as fault zones and basin centers, the overall computational cost is significantly reduced. The flexibility of the unstructured grid adapts to complex geometries, and the adaptive mechanism realizes the optimal allocation of computational resources through the "allocate as needed" principle, providing a feasible technical path for three-dimensional coupled simulations on tectonic time scales.

[0169] Solver configuration and time step control: To address the nonlinear characteristics and multi-scale time evolution requirements of deep-shallow coupled simulations, the solver parameters are finely configured and the time step coordination mechanism is optimized.

[0170] The nonlinear solver uses an iterative method to handle the thermo-mechanical coupling equations, setting the tolerance and the maximum number of iterations. This parameter design aims to balance computational accuracy and efficiency: the tolerance threshold ensures that the equation residuals converge to a reasonable range, avoiding resource waste caused by excessive iterations; the maximum iteration limit prevents infinite loops in case of non-convergence. For the high complexity of the Stokes equation, a two-stage solution strategy is adopted: first, the equation stiffness is reduced through 2000 inexpensive iterations of preprocessing to quickly approximate the approximate solution; then, 5000 high-precision iterations are used to correct the details to ensure the exact matching of the stress field and the temperature field.

[0171] The time step size is controlled differently according to the time scales of physical processes: the deep model uses an adaptive time step size, which is dynamically adjusted according to the stress change rate. When tectonic activities are intense, the time step is automatically shortened to capture transient responses; during periods of tectonic calm, the time step is extended to improve the efficiency of million-year-scale simulations.

[0172] Based on the requirements of hydrodynamic stability, the shallow model sets a fixed time step size and strictly follows the CFL condition. This limit can effectively suppress numerical oscillations caused by rapid topographic evolution (including flood erosion, landslides), ensuring the reliability of surface process simulations.

[0173] The coupling synchronization mechanism realizes the dynamic interaction between the deep and shallow modules through time series alignment and data interpolation. Every 10 shallow time steps are synchronized with 1 deep time step. The bilinear interpolation algorithm is used to map the shallow topographic data to the deep grid, and at the same time, the deep stress field is transferred to the shallow model. This design not only avoids data mismatch caused by time scale differences but also reduces numerical errors introduced by inconsistent grid resolutions through interpolation smoothing, ensuring the continuity and consistency of cross-module physical fields (including stress, temperature, topography).

[0174] Through the collaborative control of solver parameter optimization and time steps, this solution significantly improves the computational efficiency of large-scale coupled simulations while ensuring numerical stability. The adaptive time step mechanism takes into account the simulation requirements of the long-term evolution of deep structures and the transient processes near the surface, while the bilinear interpolation technique effectively solves the problem of multi-scale data transfer, providing key technical support for the fully coupled analysis of complex geological systems.

[0175] Data Output and Visualization: To support the quantitative analysis of geological processes and model verification, this solution designs a systematic data output and visualization process. The data output adopts a differential strategy to adapt to the multi-time scale characteristics of deep-shallow coupled simulations: for the deep model, the stress field (σ ij ), temperature field (T), and velocity field (v i ) are output once every deep time step, and the data is stored in VTK format. This format supports efficient compression, parallel reading and writing, and hierarchical management of large-scale 3D data, especially suitable for the storage requirements of massive physical field data in long-term tectonic evolution simulations; for the shallow model, the topographic elevation (h), sediment thickness (d), and erosion rate (E) are recorded every ten shallow time steps and stored in HDF5 format. Its data organization based on the time dimension facilitates the rapid retrieval of spatio-temporal sequences and cross-platform interaction. By setting standardized data formats, both the integrity and traceability of the data can be ensured, and the complexity of later analysis can be significantly reduced. The output data is in HDF5 and VTK formats, compatible with mainstream geological analysis software (including Paraview, GMT), and provides API interface support for seamless docking with engineering platforms (including oil and gas exploration systems, disaster warning platforms), significantly reducing the technical migration cost.

[0176] The visualization process integrates a multi-dimensional toolchain to achieve a comprehensive analysis of simulation results: 3D Visualization: Generate stress isosurface maps through Paraview, overlay the contour lines of the shallow surface topography, and mark the stress concentration areas and rockburst risk areas. Such visualization can intuitively reveal the correlation between deep tectonic activities (including strike-slip movement, stress accumulation in fault zones) and the surface topography, providing direct evidence for the location and prediction of hidden faults in engineering and disaster risk warning.

[0177] 2D Dynamic Analysis: Use Matplotlib to draw topographic evolution maps to dynamically display the basin subsidence rate, river network migration, and sediment distribution patterns. Through time series animations, the impact of river incision and sediment transport on topographic reshaping within a million-year scale can be clearly presented, assisting in the inversion of the interaction mechanism between surface processes and tectonic activities.

[0178] Model validation is a core step to ensure the scientific reliability of simulation results. This solution conducts quantitative validation by comparing the subsidence rate of the simulated fault basin with geological records (including stratigraphic thickness and radiometric dating data), and strictly controls the error threshold. If the deviation exceeds the threshold, the parameter optimization process (including adjusting rheological model parameters or boundary conditions) will be triggered until the simulation results are consistent with the observed data. This validation mechanism can not only calibrate the rationality of the initial assumptions of the model, but also provide dynamic feedback for the numerical reconstruction of complex geological processes, significantly improving the engineering applicability of the prediction results.

[0179] Through efficient data management, multi-dimensional visualization and strict verification processes, this solution realizes seamless connection from numerical simulation to geological interpretation. The standardized output format reduces the barriers of multi-disciplinary data integration, while the visualization tool chain transforms the abstract physical field data into actionable geological information through an intuitive graphical interface, providing strong technical support for research in directions such as buried fault prediction, resource exploration, and disaster assessment.

[0180] Example 2: In this example, for the fault depression zone on the western margin of a certain basin, the three-dimensional tectonic evolution process under crustal extension is simulated. The focus is on analyzing the influence of dynamic densification of the fault zone grid on the simulation accuracy of the in-situ stress field, and verifying the applicability of the shallow-deep coupling model in complex tectonic areas. Implementation process:

[0181] S1 Model configuration: Geometric model construction: A three-dimensional rectangular area with a length of 200 km × a width of 150 km × a depth of 40 km is established, including 3 main control faults (strike NNE, dip angle 60°). Spatial parameters: L = 200 km, W = 150 km, H = 40 km, and the bottom boundary is set as the Moho surface. Physical parameter setting: Deep rheological parameters: The viscosity of the upper crust η = 1×10²² Pa·s (quartzite), the lower crust η = 5× Pa·s (granulite); the internal friction angle φ = 35° (the fault zone is encrypted to 28°), and the cohesion C = 25 MPa. Shallow surface erosion equation: , and the empirical coefficient k is calibrated through data from hydrological stations on the Loess Plateau. Temperature field initialization: The surface heat flow value is set to 65 mW / m², and the radioactive heat generation rate A = 1.2 μW / m³ (uranium-rich granite layer).

[0182] S2 Grid setting, as Figure 2 shown, Initial grid generation: Unstructured tetrahedral grids are used in the deep part ( Figure 2 ), with a basic resolution of 2 km, and a pre-set encryption area (red area) for the fault zone. Delaunay triangulation is implemented on the shallow surface, with a horizontal resolution of 500 m, and areas with a terrain curvature > 0.015 are automatically encrypted to 100 m. Dynamic adjustment strategy: Strain rate threshold: When the deep units When triggered, local encryption is up to 500 m; when the surface slope change rate > 5° / 100 m, the grid size ≤ 50 m. Time step control: the initial step size in the deep part is 1000 years and is shortened to 200 years when Δσ / Δt > 0.5 MPa / yr; the fixed step size on the surface is 50 years and the CFL number is maintained at 0.8.

[0183] S3 Solver configuration, Nonlinear solver: The Newton-Krylov iteration method is adopted, with a relative tolerance of 1e-6 and a maximum number of iterations of 3000. For the preprocessing of the thermo-mechanical coupling equation, the algebraic multigrid AMG is used to reduce the condition number of the stiffness matrix. Cross-scale data transfer: Bilinear interpolation is used for the mapping of the deep stress field to the surface, the weight matrix is calculated through the element shape function, and the interpolation error is controlled within 3%. Sediment load feedback formula: η new = (1 + 0.15 d sed ,d sed When d > 20 m, viscosity correction is triggered.

[0184] S4 Run simulation, Boundary condition loading: A stretching rate of 3 mm / yr is applied to the east and west boundaries to simulate the Cenozoic regional extensional environment. The free-slip boundary is adopted at the bottom, and the initial value of the thermal convection velocity field is 1 cm / yr. Dynamic process monitoring: The real-time output shows that the fracture zone grid is encrypted 3 times at 50,000 years ( Figure 2 yellow transition zone), and the resolution is increased to 200 m. Record that the main fracture displacement rate increases from 0.8 mm / yr to 2.1 mm / yr, reflecting the tectonic activation process.

[0185] S5 Data output, 3D visualization: Generate / a stress concentration area distribution map with a value > 3.5, and mark the high-risk zone at the edge of the fault basin. Export the fracture zone grid encryption log, showing that the calculation time increases by 18% but the stress accuracy is improved by 42%. Key parameter table: The fault slip amount table contains the displacement vector (strike-slip 1.2 m, dip-slip 0.3 m), and the activation probability P = 0.68.

[0186] S6 Post-processing verification, Error analysis: Comparing with the borehole in-situ stress data, the deviation between the simulated value of the maximum principal stress of 38.5 MPa and the measured value of 42 MPa is 8.3%. Parameter optimization: The viscosity field is corrected using the Kalman filter, and after 5 iterations, the deviation is reduced to 3.7%. The MCMC calibration of the erosion coefficient x changes from 0.8 → 0.83, improving the matching degree between the simulated gully incision rate and the DEM data by 29%.

[0187] This embodiment passes through Figure 2The dynamic grid encryption strategy shown successfully captures the stress concentration phenomenon on the western margin of the fault depression basin. The simulation results show that grid encryption improves the stress calculation accuracy of the fault zone and reduces the simulation error of the migration rate of the basin subsidence center. The deep-shallow coupling mechanism accurately reproduces the co-evolution process of deep fault activation and shallow graben development under the extensional tectonic background, providing a high-precision structural model for the prediction of subtle hydrocarbon reservoirs.

[0188] Example 3: For the high-temperature geothermal field on the eastern margin of the Qinghai-Tibet Plateau, a three-dimensional geothermal-stress coupling model is constructed to analyze the influence mechanism of radioactive element heat generation on the geothermal field distribution and guide the positioning of geothermal resource target areas. Implementation process:

[0189] S1 Model configuration, geometric model construction: A 120km×80km×30km three-dimensional model is established, including granite bedrock (uranium content 8ppm) and sedimentary cover (uranium content 2ppm). Spatial parameters: L = 120km, W = 80km, H = 30km, and the surface elevation is initialized according to SRTM data. Temperature field configuration:

[0190] Upper crust temperature formula: T = 273 + 0.065 / 2.5*(20 - z) - 1.5e-6*(20 - z)² / (2*2.5), where z is in km. Radioactive heat generation layer setting: The uranium-rich granite layer (A = 3.5 μW / m³) has a thickness of 8km and a burial depth of 2 - 10km. Physical parameters: The rock thermal conductivity is set in layers: 2.0 W / (m·K) for the sedimentary layer, 2.8 W / (m·K) for granite, and 3.5 W / (m·K) for the mantle. The thermal expansion coefficient , which affects the calculation of thermal stress.

[0191] S2 Grid setting, thermodynamics coupling grid: In the deep part, a hexahedron-dominated hybrid grid is used, and the vertical resolution in the thermal anomaly area is 100m. Near the surface, a terrain-following grid is implemented, and it is stratified and encrypted to 20m within 1km of the surface to accurately depict the geothermal gradient. Adaptive strategy: Temperature gradient threshold: When ΔT / Δz > 45°C / km, the horizontal grid is automatically encrypted to 500m. Time step control: Implicit time integration is used for heat conduction, with a maximum time step of 500 years; the convective term is treated explicitly, and the CFL number is 0.6.

[0192] S3 Solver configuration, thermal-mechanical coupling solution: A segregated solution strategy is adopted. First, the temperature field is iterated until convergence (residual < 1e-4), and then the thermal stress field is calculated. The Boussinesq approximation is enabled, and the density change , and the buoyancy term is added to the momentum equation. Radioactive heat source treatment: Heat generation rate spatial distribution function: , (attenuation depth of the uranium-rich layer). Heat flow boundary: The bottom heat flow value is set to 80 mW / m², and the lateral sides use adiabatic boundaries.

[0193] S4 runs the simulation with initial conditions loaded: the geothermal temperature field is initialized, the surface temperature is 10°C, and the Moho temperature is set to 600°C. The thermal convection velocity field is initialized to 1 cm / yr to simulate the influence of mantle upwelling. Dynamic evolution monitoring: Running for 100,000 years shows that the geothermal gradient in the thermal uplift area reaches 58°C / km, with a deviation of <5% from the measured geothermal well data. The output shows that a heat flux anomaly area (>100 mW / m²) is formed at the top of the granite body, with an area of 350 km².

[0194] S5 Data output and 3D visualization: Generate a superimposed map of geothermal isosurfaces and topography, and mark the geothermal anomaly target area (the temperature in the red area > 180°C). Output the distribution of the thermal stress tensor, showing that the maximum thermal stress of 18 MPa is located at the intersection of the buried faults. Key parameter derivation: The geothermal resource inventory shows that the heat stored in the shallow layer within 3 km reaches 5.6× J, which is equivalent to the heat generated by a 1500 MW power station operating for a hundred years.

[0195] S6 Post-processing verification and error analysis: Comparing the geothermal well temperature measurement curves, the simulated temperature at a depth of 3000 m is 218°C, with a deviation of 3.1% from the measured value of 225°C. Parameter optimization: The MCMC method is used to calibrate the radiogenic heat production rate. After 7 iterations, A0 changes from 3.2 → 3.5 μW / m³. The Kalman filter corrects the thermal conductivity field, and the thermal conductivity of granite changes from 2.8 → 2.7 W / (m·K), reducing the error to 2.3%.

[0196] This embodiment characterizes the basic characteristics of the thermodynamic structure of the target area through 3D geothermal temperature field simulation. The model reproduces well the vertical variation characteristics of the measured geothermal gradient, reflecting the thermodynamic differences between the rigid layer of the upper crust and the viscoplastic layer of the lower crust. The delineation results of the shallow geothermal target area are in good agreement with the drilling verification data, confirming the basic characterization ability of the model for the in-crust thermal-mechanical coupling process. This simulation provides a new visualization analysis method for understanding the spatial coupling relationship between crustal layered deformation and thermal anomalies.

[0197] Example 4: This example aims at the problem of instability of the surrounding rock in the deep roadway of a certain metal mine. The method of the present invention is used to construct a 3D in-situ stress model for a mining area of one kilometer, and analyze the superposition effect of tectonic stress and mining-induced stress. The simulation area contains 3 buried thrust faults, with a maximum burial depth of 1200 m, and focuses on analyzing Figure 4 the stress distribution characteristics of the longitudinal section shown. The implementation process is as follows:

[0198] S1 Model Configuration; Geometric Model: A three-dimensional area with a length of 3 km × width of 2 km × depth of 1.2 km is constructed. The strike of the main ore body is N30°E and the dip angle is 55°. Fault parameters are defined: The strike of Fault F1 is 110°, the dip angle is 65°, and the fault throw is 80 m. Physical parameters: The viscosity of deep granite η = 8×10²¹ Pa·s, the internal friction angle φ = 38°, and the cohesion C = 45 MPa; The viscosity of the shallow altered rock mass is reduced to . The temperature field adopts the measured geothermal gradient of 2.8℃ / 100m, and the maximum heat production rate A = 2.1 μW / m³.

[0199] S2 Mesh Setting; Initial Mesh: Unstructured tetrahedral meshes are used in the deep part, with a basic resolution of 50 m, and the area around the roadway is encrypted to 10 m ( Figure 4 yellow area). The shallow surface topography adopts Delaunay triangulation, with a horizontal resolution of 20 m. Dynamic Optimization: Set the strain rate threshold , and when the strain of the surrounding rock exceeds the limit due to mining disturbance, it will be automatically encrypted to a 5 m mesh.

[0200] S3 Solver Configuration; The hybrid Newton - quasi - Newton algorithm is adopted, with the maximum number of iterations set to 3000 times and the tolerance to 1e - 5. Cubic spline interpolation is used for cross - scale data transfer, and the stress field mapping error < 3%. The equivalent stiffness matrix of the support structure is introduced, and the bolt support is equivalent to anisotropic strengthening parameters implanted into the constitutive equation.

[0201] S4 Running Simulation; Boundary Conditions: Horizontally applied tectonic stress σh = 1.5σv is applied laterally to simulate the squeezing environment of the simulation area. Real - time output display shows that when mining reaches the - 800 m level, a stress concentration area ( Figure 4 red area) appears on the hanging wall of Fault F1, and the maximum principal stress reaches 58 MPa. A three - dimensional stress nephogram ( Figure 4 ) is generated, clearly showing the gradient distribution along the fault on the longitudinal section, and a high - strain area (ε > 5%) with a range of 20 m is formed at the fault intersection.

[0202] S5 Data Output; Export the distribution of the stress ratio of the key section ( / ), and mark the potential rockburst areas where / > 4. Generate an optimized support parameter table, and recommend that high - pre - tension bolts (300 kN) be densely arranged (spacing 0.8 m) in the stress concentration area.

[0203] S6 Post - processing Verification; Comparing with the microseismic monitoring data, the spatial matching degree between the simulated stress concentration area and the measured microseismic events reaches 87%. By using the Kalman filter to correct the φ value of the fault zone from 35° to 32°, the prediction error of the roadway deformation is reduced from 12% to 6.5%.

[0204] Figure 4 It shows the longitudinal section position of the engineering core area, which can clearly display the spatial layout of the strata and the in-situ stress as a whole and important structural features, providing an intuitive reference for subsequent analysis.

[0205] Example 5: In view of the engineering problem that a certain reservoir dam foundation crosses the F2 active fault, the method of the present invention is adopted to realize the three-dimensional in-situ stress field reconstruction with meter-level accuracy in the dam area. The implementation process is as follows:

[0206] S1 Model configuration; Geometric model: A 4km×3km×1.5km three-dimensional model is established, and the main dam axis intersects the fault strike at an angle of 35°. Define the F2 fault parameters: strike N45°E, dip 70°, slip rate 2.8mm / yr. Physical parameters: The fault zone adopts a strain-softening model, peak φ = 28°, residual φ = 18°; the elastic modulus of the dam concrete E = 30GPa, Poisson's ratio ν = 0.2.

[0207] S2 Mesh setting; Initial discretization: The overall mesh resolution is 100m, and a preset encryption area is set in the fault zone, and an unstructured hexahedron-tetrahedron hybrid mesh is used. Dynamic optimization: When the local formed curvature > 0.02 or strain gradient > 5% / m, trigger local encryption to 5m mesh, and a total of 4 levels of adaptive refinement are completed.

[0208] S3 Solver configuration; Adopt an explicit-implicit coupling algorithm: The dam body adopts an explicit central difference method (time step 0.1s), and the rock mass adopts an implicit Newmark-β method. Set contact surface elements to simulate the interaction between the fault and the dam foundation, and the friction coefficient μ is dynamically updated (0.6 → 0.4).

[0209] S4 Run simulation; Apply seismic load (PGA = 0.3g), and the simulation shows that the fault dislocation triggers the stress redistribution of the dam foundation. Figure 5 It shows that the refined mesh accurately captures the tensile stress area within the fault influence zone. Record the shear stress mutation process at the anti-seepage wall joint, and the peak value τ = 4.2MPa exceeds the design value by 23%.

[0210] S5 Data output; Generate a three-dimensional tensile stress isosurface and mark the dangerous area with < -5MPa. Export the time series curve of the slip amount at the dam-fault contact surface, and the maximum relative displacement reaches 12cm. Based on the Griffith criterion, it is recommended to add seismic hinges in the fault influence zone and increase the thickness of the anti-seepage wall from 2m to 2.5m.

[0211] S6 post-processing verification; by comparing distributed fiber optic monitoring data, the coincidence degree between the simulated dam foundation deformation mode and the measured value reaches 91%. The μ value of the fault is calibrated by the MCMC method, reducing the prediction error of the joint opening from 9 mm to 3 mm.

[0212] Figure 5 It is the refined simulation result of the in-situ stress in the dam foundation area across the fault in Embodiment 5. Figure 5 It can be seen that during the data assimilation process, after adaptive mesh refinement near the fault, the three-dimensional in-situ stress of the calculation area shows a more realistic in-situ stress distribution, enhancing the stress prediction ability of the local area and providing a more scientific basis for subsequent engineering design, structural optimization, etc.

[0213] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic bidirectional coupling, characterized in that, It includes the following steps: S1. Model configuration: Construct a three-dimensional geometric model based on the geological structure characteristics of the target area, define the elastoplastic rheological parameters of deep rocks, establish a surface erosion rate equation, and set the gravity field parameters, radioactive heat generation parameters, and layered temperature field configuration; S2. Mesh setting: Generate deep unstructured meshes and surface irregular meshes using a differential mesh strategy, and set boundary conditions, an adaptive time step control strategy, and a mesh adaptive optimization strategy; S3. Solver configuration: Configure a nonlinear solver to handle the thermo-mechanical coupling equation, establish a staged solution strategy, configure a deep and shallow dynamic two-way coupling mechanism, and achieve cross-scale data transfer through a bilinear interpolation algorithm; S4. Run simulation: Output physical field data including viscosity, density, lithostatic pressure, tectonic stress, and strain rate, and record terrain elevation evolution data; S5. Data output: Generate a superimposed visualization map of three-dimensional stress isosurfaces and terrain contour lines, mark the principal stress distribution and stress concentration areas, and export a key parameter table; S6. Post-processing: Calculate the error between the simulation results and geological observation data, and trigger a parameter optimization loop when the error exceeds the threshold; Use a filtering algorithm to correct the viscosity, friction coefficient, and calibrate the erosion equation coefficient; Generate a fault location and activation risk assessment report.

2. The three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic two-way coupling according to claim 1, wherein: In S1, the construction of the three-dimensional geometric model includes: based on the actual geological background and geographical environment of the target area, defining the length L, width W, and depth H of the model as the three-dimensional spatial range of the target research area from the deep crust to the shallow surface; among them, in the elastoplastic rheological parameters of deep rocks, the unit of viscosity η is Pa·s, which characterizes the flow resistance of rocks under high stress; the unit of the internal friction angle φ is degree, which is determined by rock mechanics experiments and is used to describe the non-linear relationship between the shear strength of rocks and the change of normal stress; the unit of cohesion C is MPa, which characterizes the bonding force between rock particles and is calibrated by triaxial compression experiments; the erosion rate equation of the shallow geodynamics model is , where Q is the surface water flow rate (m³ / s), S is the slope, k is an empirical coefficient ( ), x and y are exponential parameters fitted based on regional hydrological observation data; the gravitational acceleration g is set to 9.81 m / s², and the direction is vertically downward; the unit of the radioactive heat generation rate A is μW / m³, which is calculated and determined according to the abundance of radioactive elements and the decay heat generation rate in the target area; among them: The elastoplastic rheological parameters adopt the Drucker-Prager model, and its yield criterion is: ， Among them, is the second invariant of the deviatoric stress, is the first invariant of the stress tensor, , ; is the angle of internal friction, an angular parameter characterizing the shear strength of rocks, reflecting the internal friction characteristics of materials during shear failure; this model characterizes the creep behavior of deep rocks in a high-temperature environment by introducing temperature-dependent viscosity η(T), and the viscosity expression is: ， wherein, is the reference viscosity, is the activation energy, is the gas constant, is the absolute temperature (K).

3. The three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic two-way coupling according to claim 2, wherein: In S1, the temperature field configuration includes: The upper crust temperature is jointly determined by heat flow and radioactive heat generation, and the calculation formula is: where T is the absolute temperature, z is the depth, 0.065 is the heat flow density, 2.5 is the thermal conductivity of the rock, 1.5× is the radioactive heat generation rate; The lower crust temperature linearly decreases to 1573 K at the bottom of the lithosphere.

4. The three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic two-way coupling according to claim 1, wherein: In S2, the generation of the deep unstructured grid includes: initially discretizing the deep model using tetrahedral elements, with the element size dynamically adjusted according to the fracture zone density, and the grid resolution in the fracture zone area being encrypted to the order of 100 meters; the shallow irregular grid is generated by the Delaunay triangulation algorithm, with the horizontal resolution set to 50 - 100 meters, and according to the terrain curvature, the curvature threshold is set to , triggering local encryption to ensure the grid refinement in the terrain mutation area; the horizontal tectonic movement rate of the lateral boundary of the deep model is set according to the regional plate movement data, ranging from 0.1 - 10 mm / yr, to simulate the crustal extension or extrusion process.

5. The three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic two-way coupling according to claim 1, wherein: In S2, the adaptive time step control strategy includes: the deep dynamic step size is adaptively adjusted according to the stress change rate When MPa / yr, the step size is shortened to less than 100 years; when MPa / yr, the step size is extended to 1000 years; the shallow fixed step size follows the CFL condition, that is: ， wherein, is the minimum grid size (m), is the maximum surface water flow velocity (m / s) to ensure the numerical stability of terrain evolution; The grid adaptive optimization strategy includes: when the deep local strain rate exceeds grid encryption is triggered, and the grid resolution in the fracture zone area is increased to 200 meters; for the shallow model, local refinement of Delaunay triangulation is triggered according to the ground formation curvature threshold to ensure that the grid size is ≤ 50 meters in areas where the slope change rate ≥ 5° / 100 meters.

6. The three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic two-way coupling according to claim 1, wherein: In S3, the staged solution strategy includes: in the first stage, the algebraic multigrid (AMG) preprocessing technology is adopted to reduce the residual of the Stokes equation to 10% of the initial value; in the second stage, the Newton-Krylov iteration method is adopted, with a tolerance set to 1× , and the maximum number of iterations is 5000; for cross-scale data transfer, the bilinear interpolation algorithm is adopted to map the stress field of the deep kilometer-level grid to the shallow hundred-meter-level grid. The interpolation weight is calculated according to the element shape function, and the error is controlled within 5%; In S3, the process of configuring the shallow and deep coupling mechanism includes: the thickness d of shallow deposition sed Dynamically feedback the deep viscosity field through the ballast effect, and the viscosity correction formula is: ， where d sed is in meters, is the initial viscosity (Pa·s); the deep horizontal tectonic stress σ xy is transferred to the shallow model through bilinear interpolation to drive the topographic elevation h to calculate the settlement amount according to the formula Δh = σ xy / (ρg), where ρ is the rock density (kg / m³) and g is the acceleration due to gravity (m / s²).

7. The three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic two-way coupling according to claim 1, wherein: In S4, the lithostatic pressure is determined by calculating with the formula where: is the rock density, representing the mass distribution characteristics of the rock mass; is the acceleration due to gravity, with the direction vertically downward; is the current depth, representing the vertical distance from the target point to the surface of the earth; Tectonic stress , representing the dynamic stress field generated by crustal tectonic movement, is obtained by solving the Stokes equation and is superimposed with the lithostatic pressure to obtain the total pressure field ; Hidden fault displacement By calculating the time integral of the velocity field That is , the accuracy reaches the centimeter level; where Is the starting time of the simulation; Is the current simulation time.

8. The three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic two-way coupling according to claim 1, wherein: In S5, the generation of the three-dimensional stress isosurface includes: the maximum principal stress and the minimum principal stress are obtained through eigenvalue decomposition, and the regions with the principal stress ratio are marked as stress concentration areas; the rock burst risk areas are determined by the Griffith criterion, that is: ， Among them, is the tensile strength of the rock (MPa), and the threshold is set to 5 MPa; the displacement of the fault slip amount is recorded in the table , slip rate and activation probability , where is calculated based on the ratio of shear stress and normal stress .

9. The three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic two-way coupling according to claim 1, wherein: In S6, the parameter optimization loop includes: correcting the viscosity using the Kalman filter algorithm and the internal friction angle , with the observed data being the borehole in-situ stress measurement values; the erosion equation coefficients , , are calibrated by the Markov Chain Monte Carlo (MCMC) method, the objective function is the root mean square error (RMSE) between the simulated terrain elevation and the LiDAR measured data, and the threshold is set at ±2 meters; when the error exceeds the limit, 10 iterative optimizations are triggered until RMSE < 1 meter; The fault location report includes determining the spatial coordinates, strike angle , dip angle and slip type; the activation risk assessment is based on the probability distribution of the historical earthquake catalog and the simulated slip rate, and the risk levels are divided into low, medium, and high, and a risk heat map and emergency plan suggestions are output.

10. Use of the three-dimensional in-situ stress numerical simulation method based on deep and shallow dynamic two-way coupling according to any one of claims 1-9 in buried fault location and resource exploration.

Images