A Three-Dimensional Simulation Method for Ore-Forming Fluids Based on Multi-Field Coupling Effect

Abstract

This invention relates to the field of computer-aided engineering applications in geological science, and particularly to a three-dimensional simulation method for ore-forming fluids based on multi-field coupling effects. The method includes the following steps: acquiring and calculating critical temperature thresholds based on viscosity-temperature response for grid cells in a three-dimensional geological grid model data, then classifying the grid cell states to obtain grid state labeling data; expanding geological grid attributes based on the grid state labeling data to obtain a viscosity-temperature partitioned grid; constructing a simulation environment using the viscosity-temperature partitioned grid and extracting temperature thresholds to form comprehensive temperature threshold data; detecting state changes in the comprehensive temperature threshold data to obtain state-change cells; and analyzing and recording pulse events in the viscosity-temperature partitioned grid to obtain flow pulse event sequences. This invention, through computer-aided engineering recording of fluid migration phenomena, significantly improves the scientific rigor and accuracy of ore-forming target area prediction.

Description

Technical Field

[0001] This invention relates to the field of computer-aided engineering applications in geological science, and in particular to a three-dimensional simulation method for ore-forming fluids based on multi-field coupling effects. Background Technology

[0002] Existing technologies, when simulating the dynamic behavior of geological fluids, simplify the viscosity of ore-forming fluids to a constant value or a simple function of temperature (such as a linear or exponential function). This simplification fails to reflect the nonlinear abrupt changes in viscosity of real geological fluids, which can drop by several orders of magnitude within a specific temperature range. This leads to fundamental deviations between simulation results and actual geological processes. Traditional fluid simulation methods, based on the assumption of a continuous medium and smoothly changing parameters, cannot capture abrupt changes in fluid states, especially the instantaneous transition from a "high-viscosity stagnant state" to a "low-viscosity, high-flowability state." This makes it impossible for models to explain the "pulse-like" or "episode-like" fluid migration phenomena commonly found in geological records, which are precisely the key processes controlling the formation of major mineralization events. Existing mineralization prediction models often treat fluid events as independent entities, ignoring the spatiotemporal coupling relationships and cumulative enhancement effects between multiple fluid events. However, the accumulation of ore-forming materials in actual geological processes is often formed by the combined action of multiple interconnected fluid events through "chain reactions" or "cascade triggering." This simplification significantly reduces the accuracy and reliability of prediction models.

[0003] In summary, existing technologies lack comprehensive ore-forming fluid simulation methods that can accurately simulate abrupt changes in fluid state, capture pulsed transport characteristics, and consider the coupling effects of multiple events. This severely limits the accuracy and practical value of mineral exploration prediction. Summary of the Invention

[0004] Therefore, it is necessary to provide a three-dimensional simulation method for ore-forming fluids based on multi-field coupling effects to solve at least one of the above-mentioned technical problems.

[0005] To achieve the above objectives, a three-dimensional simulation method for ore-forming fluids based on multi-field coupling effects is proposed, comprising the following steps:

[0006] Step S1: Obtain and calculate the critical temperature threshold based on viscosity-temperature response for the grid cells of the 3D geological grid model data, then classify the grid cell states to obtain grid state label data; expand the geological grid attributes based on the grid state label data to obtain a viscosity-temperature partitioned grid.

[0007] Step S2: Construct a simulation environment using a viscosity-temperature partitioned grid and extract temperature thresholds to form comprehensive temperature threshold data; perform state change detection on the comprehensive temperature threshold data to obtain state change units; analyze and record pulse events in the cells of the viscosity-temperature partitioned grid that have undergone state changes based on the state change units to obtain a flow pulse event sequence;

[0008] Step S3: Construct an initial flux evaluation grid based on the pulse event data, identify the event cascade effect, and obtain the event superposition enhancement factor; superimpose the event superposition enhancement factor onto the initial flux evaluation grid to simulate the cumulative quantization of the fluid transport channel efficiency within the simulation time period, forming a cumulative transport flux field;

[0009] Step S4: Perform fluid flux and ore-forming condition matching analysis on the cumulative transport flux field to obtain the ore-forming condition matching degree field; delineate and visualize the ore-forming potential target area based on the ore-forming condition matching degree field to obtain a three-dimensional ore-forming probability cloud map.

[0010] This invention, by introducing a piecewise nonlinear viscosity response function and a "three-state partitioning" mechanism, accurately quantifies the nonlinear abrupt change in the viscosity of geological fluids within the critical temperature range, overcoming the simplification limitations of traditional simulations and laying the physical foundation for simulating instantaneous fluid state transitions. By real-time monitoring of cell state changes, instantaneous amplification of permeability, and reconstruction of the local pressure field, this method accurately captures and quantitatively simulates the instantaneous transition of fluids from stagnation to high fluidity. The generated sequence of flow pulse events records "pulsating" flow information in the form of "spatiotemporal fingerprints," providing a new perspective for understanding the correlation between major mineralization events and fluid dynamics.

[0011] This invention innovatively constructs an event spatiotemporal correlation matrix and trigger chain analysis to identify coupling and cascading enhancement effects among multiple events. By quantifying the "resonance enhancement region" and generating superimposed enhancement factors, it more realistically reflects the formation of dominant fluid channels under the synergistic effect of multiple pulse events. The resulting cumulative transport flux field transforms discrete events into a macroscopic evaluation of the fluid transport "highway," providing direct physical evidence for accurately locating ore-forming accumulation paths. Furthermore, this invention, by calculating the mineralization condition matching degree field, closely integrates fluid dynamics simulation with the physicochemical conditions of mineralization (such as cooling rate), overcoming the one-sidedness of relying solely on fluid flux prediction. The weighted fusion of the flux field and the matching degree field ensures that the prediction results take into account both "quantity" and "quality," significantly improving the accuracy of target area delineation. The final generated three-dimensional mineralization probability cloud map intuitively displays the mineralization potential and prediction reliability, greatly improving the accuracy and efficiency of mineral exploration prediction.

[0012] In summary, this invention, through systematic and innovative simulation of fluid property abrupt changes, pulsed migration events, multi-event cascade effects, and mineralization conditions, constitutes a complete solution that can accurately reproduce the "pulsed" fluid migration phenomenon in the geological record, significantly improving the scientific rigor and accuracy of mineralization target area prediction. Attached Figure Description

[0013] Figure 1 This is a schematic diagram of the steps in a three-dimensional simulation method for ore-forming fluids based on multi-field coupling effects. Detailed Implementation

[0014] The objectives, features, and advantages of this invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings.

[0015] The technical method of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0016] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0017] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0018] To achieve the above objectives, please refer to Figure 1 This invention provides a three-dimensional simulation method for ore-forming fluids based on multi-field coupling effects, comprising the following steps:

[0019] Step S1: Obtain and calculate the critical temperature threshold based on viscosity-temperature response for the grid cells of the 3D geological grid model data, then classify the grid cell states to obtain grid state label data; expand the geological grid attributes based on the grid state label data to obtain a viscosity-temperature partitioned grid.

[0020] In this embodiment of the invention, three-dimensional geological grid model data is analyzed. For different lithologies, piecewise nonlinear viscosity-temperature response functions are defined, and the flow threshold temperature T_p and the complete flow temperature T_m are accurately calculated through derivative analysis. Subsequently, based on the relationship between the current temperature of each grid cell and these critical thresholds, it is divided into one of three states: "high viscosity stagnant zone," "critical activation zone," or "low viscosity free flow zone," and its critical proximity and state stability index are calculated. Finally, these state labels and dynamic viscosity parameters are integrated into the original grid as extended attributes to construct an innovative viscosity-temperature partitioned grid.

[0021] Step S2: Construct a simulation environment using a viscosity-temperature partitioned grid and extract temperature thresholds to form comprehensive temperature threshold data; perform state change detection on the comprehensive temperature threshold data to obtain state change units; analyze and record pulse events in the cells of the viscosity-temperature partitioned grid that have undergone state changes based on the state change units to obtain a flow pulse event sequence;

[0022] In this embodiment of the invention, a viscosity-temperature partitioned grid is used as the initial state, and the three-dimensional temperature field is iteratively updated within discrete time steps. At each time step, by comparing the temperature with the critical threshold at previous and subsequent times, cells where state label transitions occur are detected and identified, and spatially adjacent strongly changing cells are clustered into macroscopic "state change cells". For these cells, a pulse event analysis process is triggered: the permeability is instantaneously and nonlinearly amplified, the local pressure field is reconstructed, and the resulting pulse flux vector is calculated. Finally, the spatiotemporal characteristics of all pulse events and their chain reaction propagation are recorded to form a flow pulse event sequence.

[0023] Step S3: Construct an initial flux evaluation grid based on the pulse event data, identify the event cascade effect, and obtain the event superposition enhancement factor; superimpose the event superposition enhancement factor onto the initial flux evaluation grid to simulate the cumulative quantization of the fluid transport channel efficiency within the simulation time period, forming a cumulative transport flux field;

[0024] In this embodiment of the invention, a sequence of flow pulse events is analyzed, and a cumulative flux evaluation grid is initialized. By analyzing the spatiotemporal correlation matrix and event triggering chains between events, a "resonance enhancement region" is identified where multiple pulse events synergistically enhance the nonlinear fluid flux, and the event superposition enhancement factor is calculated accordingly. Then, the flux contribution of each pulse event (considering distance and structural anisotropy) is multiplied by the corresponding enhancement factor and accumulated over the entire simulation time. Finally, by thresholding and connectivity analysis of the accumulation results, a final cumulative transport flux field is generated, which clearly outlines the dominant fluid transport driven by pulse events.

[0025] Step S4: Perform fluid flux and mineralization condition matching analysis on the cumulative transport flux field to obtain the mineralization condition matching degree field; delineate and visualize the mineralization potential target area based on the mineralization condition matching degree field to obtain a three-dimensional mineralization probability cloud map.

[0026] In this embodiment of the invention, the temperature and pressure history, cooling efficiency, and fluid residence time during the simulation process are analyzed and compared with a preset ore-forming condition template to calculate a three-dimensional ore-forming condition matching degree field. Next, this matching degree field is weighted and fused with the cumulative transport flux field to generate a comprehensive ore-forming potential field. Based on multi-level threshold segmentation and spatial connectivity analysis of this potential field, ore-forming potential target areas are automatically delineated, optimized, and ranked. After uncertainty assessment, all results are finally integrated, and a visually intuitive three-dimensional ore-forming probability cloud map is generated using volumetric rendering technology, where color represents the level of potential and transparency represents the reliability of the prediction.

[0027] Preferably, step S1 includes the following steps:

[0028] Step S11: Analyze the three-dimensional geological grid model data to generate basic geological grid data;

[0029] Step S12: Based on different lithological types and pressure conditions in the geological grid base data, define viscosity-temperature response curves to form viscosity response functions;

[0030] Step S13: Determine the critical temperature threshold for the viscosity response function to obtain the critical temperature threshold;

[0031] Step S14: Traverse the cells of the geological grid basic data, classify the grid cell status according to the critical temperature threshold, and obtain grid status label data;

[0032] Step S15: Extend the geological grid base data with viscosity-temperature properties using grid state label data and viscosity response function to obtain a viscosity-temperature partitioned grid.

[0033] In an embodiment of the present invention, three-dimensional geological grid model data including lithology, structure, and initial temperature and pressure fields is received. This data contains 100×100×100 hexahedral grid cells. A spatial index structure is used to parse the input data and convert it into a three-dimensional array structure to ensure efficient access. For each grid cell, its spatial coordinates (x, y, z), lithology type identifier lithID, distance to the nearest fault zone distFault, initial temperature T_init, and initial pressure P_init are extracted. lithID is converted into a specific set of lithology physical parameters through a predefined mapping table, including thermal conductivity λ, specific heat capacity C_p, porosity φ, and base permeability k_0. The temperature field and pressure field data are respectively organized into three-dimensional arrays T(i, j, k) and P(i, j, k). The finally generated geological grid basic data contains the complete physical properties of each grid cell and maintains the topology of the original grid.

[0034] According to each lithology type and pressure condition extracted from the geological grid basic data, a unique non-linear viscosity-temperature response curve is defined. This curve uses a three-segment function to describe the change of fluid viscosity μ with temperature T: when T < T_p, μ = μ_high, representing the "high-viscosity stagnant zone"; when T_p ≤ T < T_m, μ = μ_high×exp(-α×(T - T_p) / (T_m - T_p)), representing the rapid decrease of viscosity in the "critical activation zone"; when T ≥ T_m, μ = μ_low, representing the "low-viscosity free-flow zone". Among them, μ_high takes the value of 10³ Pa·s, μ_low takes the value of 1 Pa·s, and α is the viscosity mutation coefficient taking the value of 10. For different pressure conditions, the T_p and T_m thresholds are linearly corrected, T_p = T_p° + β_p×(P - P_0), T_m = T_m° + β_m×(P - P_0), where T_p° and T_m° are the reference temperatures under the standard pressure P_0, and β_p and β_m are the pressure correction coefficients. Corresponding reference temperatures, pressure correction coefficients, and viscosity mutation parameters are configured for different lithologies (such as granite, sandstone) to form a complete set of viscosity response functions.

[0035] Numerically analyze each viscosity response function formed in step S12 to determine the critical temperature points at which the viscosity changes drastically. Within the preset temperature range [50°C, 500°C], calculate the viscosity value sequence μ(T) and its first derivative sequence dμ / dT with a step size of 1°C. Identify the temperature intervals where the absolute value of dμ / dT exceeds a specific threshold (e.g., dμ / dT < -0.05 × μ_high). Based on this analysis, accurately determine the flow threshold temperature T_p and the full flow temperature T_m for each lithology type under a given pressure. T_p is defined as the starting temperature point at which the viscosity begins to decrease significantly, and T_m is defined as the temperature point at which the viscosity reaches a stable low value (e.g., μ(T) < 1.1 × μ_low). For fluids with multiple viscosity mutation stages, record all eligible (T_p, T_m) temperature pairs. Finally, the critical temperature pairs of all lithology types under different pressure conditions are organized into a critical temperature threshold table.

[0036] Traverse each cell in the geological grid basic data to obtain its current temperature T, lithology type lithID, and pressure P. Query the corresponding T_p and T_m from the critical temperature threshold table according to lithID and P. Based on the comparison relationship between T and T_p and T_m, mark the cell with the status code "1" (T < T_p, high viscosity stagnation area), "2" (T_p ≤ T < T_m, critical activation area), or "3" (T ≥ T_m, low viscosity free flow area). At the same time, calculate the critical proximity parameter δ = (T - T_p) / (T_m - T_p) to quantify the relative position of its temperature in the critical interval. Analyze the status distribution of 26 neighboring cells around each cell to form a three-dimensional status frequency vector F = [f_1, f_2, f_3]. Based on F and δ, calculate the status stability index S = 1 - H(F) - |δ - 0.5|, where H(F) is the normalized entropy value. The finally generated grid status marking data includes the status code, δ value, and S value of each cell.

[0037] Integrate the grid status marking data (including status code, critical proximity δ, status stability index S) with the relevant parameters (μ_high, μ_low, α, T_p, T_m) in the viscosity response function set into the geological grid basic data. For each cell in the geological grid basic data, based on the original geological attributes, add the following extended attribute fields: the status label indicating the viscosity interval it belongs to, the T_p and T_m values corresponding to this cell, the viscosity response function parameters μ_high, μ_low, and α, as well as the calculated δ and S. In addition, calculate and store the current viscosity value μ_current of this cell according to the current temperature and the viscosity response function. For large three-dimensional models, adopt a sparse storage structure and a spatial octree index to optimize the storage and access efficiency of the viscosity-temperature partition grid. The finally generated viscosity-temperature partition grid accurately characterizes the non-linear characteristics of the fluid viscosity of each grid cell and its potential dynamic change information.

[0038] Preferably, step S14 includes:

[0039] Read the temperature field data and lithology distribution data in the geological grid basic data, perform basic state division according to the critical temperature threshold, and obtain the preliminary state classification result;

[0040] Calculate the degree of the cell's distance from the state transition in the preliminary state classification result to form the critical area sensitivity data;

[0041] Analyze the neighborhood state structure of the cells in the preliminary state classification result to form the neighborhood state structure;

[0042] Perform state stability evaluation according to the neighborhood state structure and the critical area sensitivity data to obtain the state stability index;

[0043] Integrate the critical area sensitivity data and the state stability index to generate the grid status marking data.

[0044] In the embodiment of the present invention, read the temperature T and lithology lithID of the geological grid basic data. Query the corresponding flow threshold T_p and complete flow temperature T_m from the critical temperature threshold table according to lithID. Subsequently, compare the T of each cell with T_p and T_m, and divide it into a "high-viscosity stagnant area" (T < T_p), a "critical activation area" (T_p ≤ T < T_m), or a "low-viscosity free-flow area" (T ≥ T_m) to obtain the preliminary state classification result.

[0045] For the cells in the "critical activation area", calculate the critical proximity parameter δ = (T - T_p) / (T_m - T_p) to quantify the relative position of its temperature within the critical activation interval. For the cells in other areas, the δ value is set to a specific value indicating being far from the critical state. All the calculated δ values constitute the critical area sensitivity data.

[0046] For each grid cell, a neighborhood analysis is performed to examine the preliminary state classification results of its 26 neighboring cells. The number of cells in the neighborhood at states "1", "2", and "3" is counted and normalized to a frequency, forming a three-dimensional state frequency vector F=[f_1,f_2,f_3], which represents the neighborhood state structure of the central cell. The F vectors of all cells are then aggregated to form the neighborhood state structure data.

[0047] After obtaining the sensitivity data of the critical region and the neighborhood state structure, the state stability of each cell is evaluated. The state stability index S = 1 - H(F) - |δ - 0.5|. Here, H(F) is the normalized entropy value of the neighborhood state frequency vector F, representing the degree of disorder in the neighborhood state; |δ - 0.5| measures the distance of the temperature from the midpoint of the critical interval. The larger the S value, the more stable the state. The S values ​​of all cells together constitute the state stability index.

[0048] The preliminary state classification results, critical region sensitivity data (δ), and state stability index (S) are integrated. The data structure for each grid cell includes: a state label (state code 1, 2, or 3), a critical proximity parameter δ, and a state stability index S. This integrated data provides potential information about the cell's current state and its dynamic behavior, generating the final grid state labeling data.

[0049] Preferably, step S2 involves constructing a simulation environment using a viscosity-temperature partitioned grid and extracting temperature thresholds, including:

[0050] Using the viscosity-temperature partitioned grid as the initial state, the spatiotemporal iterative framework is obtained and initialized according to the preset simulation parameters to obtain the initial simulation environment;

[0051] The temperature field is dynamically updated based on the initial simulation environment and viscosity-temperature partitioned grid to obtain the updated temperature field;

[0052] Temperature time-series comparison data are extracted from the updated temperature field and combined with lithological information and viscosity critical thresholds extracted from the viscosity-temperature zoning grid to form comprehensive temperature threshold data.

[0053] In this embodiment of the invention, a viscosity-temperature partitioned grid is loaded as the initial condition at time t=0. This grid includes the initial temperature and pressure, lithology, state label, and associated viscosity response function and critical temperature threshold (T_p, T_m) for each cell. The total simulation duration is set to 10,000 years, the time step Δt is 10 years, and the boundary conditions are 20°C at the surface and 60 mW / m² heat flow at the bottom. Simultaneously, an empty flow pulse event sequence is initialized. A three-dimensional implicit finite difference solver is established for iterative temperature field calculations. This process constructs a structured execution environment containing all parameters and data structures.

[0054] Within each time step Δt, the temperature field change is calculated based on the current simulation environment and the thermophysical parameters stored in the viscosity-temperature partitioned grid. The heat conduction equation is solved using a three-dimensional implicit finite difference method. The calculations consider the differences in thermal conductivity λ and specific heat capacity C_p for different lithologies (e.g., granite λ=2.5W / (m·K), C_p=800J / (kg·K); sedimentary rock λ=1.8W / (m·K), C_p=950J / (kg·K)). Simultaneously, the pressure P of each cell is updated. At the end of each time step, an updated temperature field containing the new temperature values ​​T(t+Δt) of all grid cells is obtained.

[0055] The temperature values ​​of the current time step T(t+Δt) and the previous time step T(t) are extracted from the updated temperature field to construct temperature time series comparison data ΔT_time = T(t+Δt) - T(t). For each grid cell, its lithological information (lithID) and the viscosity critical thresholds T_p and T_m determined by that lithology are extracted from the viscosity-temperature partitioned grid. These thresholds indicate the key temperature points of fluid state transitions. ΔT_time, the current temperature T(t+Δt), T_p, T_m, and lithID are integrated to form comprehensive temperature threshold data. This data structure provides necessary information for subsequent state change detection, such as temperature change trends, the current temperature range, and the benchmark point for determining state transitions, ensuring accurate capture of abrupt changes in fluid viscosity.

[0056] Preferably, step S2, which involves detecting state changes in the comprehensive temperature threshold data, includes:

[0057] Determine whether temperature changes in the comprehensive temperature threshold data have led to changes in the status label, and obtain candidate change units;

[0058] Candidate change units are classified by transformation type to obtain classified change units;

[0059] A quantitative assessment of the change intensity of the classification change units is performed to obtain an intensity assessment change table;

[0060] The intensity assessment change table is divided into change clusters and state change units are generated.

[0061] In this embodiment of the invention, each grid cell in the temperature threshold integrated data is traversed to obtain its previous time step state label (State_old), current temperature T(t+Δt), and corresponding T_p and T_m. Based on T(t+Δt), its new state label (State_new) is re-determined. If State_new is inconsistent with State_old, the cell is identified as a candidate change cell, and its coordinates, State_old, State_new, current temperature, and temperature change rate ΔT_time are recorded. For example, a cell that changes from a stagnant region (State_old=1) to a critical activation region (State_new=2) due to a temperature increase above T_p is listed as a candidate change cell.

[0062] All candidate change units are meticulously categorized based on their State_old and State_new combinations. The main transition types include: Type A (stagnant to critical), Type B (critical to free), Type C (stagnant to free, the most drastic transition), and Type D (free to critical, etc., reverse transitions). Each candidate change unit is labeled with its corresponding transition type, forming a categorized change unit list. For example, Type C cells are additionally labeled "high priority" to indicate their potentially greatest impact on fluid transport.

[0063] Based on the transformation type, temperature change rate ΔT_time, and critical proximity parameter δ of the classification change unit, the intensity of the change is quantitatively assessed. For "positive" transformations of types A, B, and C, the intensity assessment function is F_intensity = (ΔT_time / Δt) × (State_new - State_old) × f_nonlinearity(δ). Here, f_nonlinearity(δ) is a nonlinear amplification function with a high value when δ approaches the critical point. Transformations of type C are additionally multiplied by a weighting coefficient (e.g., 5.0) to reflect their stronger instantaneous impact. The F_intensity for reverse transformations is set to a smaller negative value. The F_intensity value of each cell, along with its coordinates and transformation type, is recorded to form an intensity assessment change table.

[0064] Spatial clustering algorithms (such as DBSCAN or connected component labeling algorithms) are used to process cells in the intensity assessment change table. A spatial distance threshold (e.g., the side length of two grid cells) and a minimum number of clustering units (e.g., five cells) are set to group spatially adjacent change cells with intensity above the average into change clusters. For example, 100 spatially continuous cells that simultaneously transition from stagnant zones are identified as a large cluster. The average intensity, center coordinates, volume, and principal direction of each cluster are calculated. These clusters represent local or regional "hotspot events" that may occur during geological processes. Each cluster is abstracted as a "state change unit," and its characteristic information is recorded, ultimately generating state change unit data, transforming discrete change events into macroscopic events with spatially organized characteristics.

[0065] Preferably, step S2, which involves analyzing and recording pulse events in cells of the viscosity-temperature grid where state changes occur based on the state change units, includes:

[0066] Based on the state change unit, the permeability mutation of the viscosity-temperature partition grid is calculated to obtain the permeability mutation data;

[0067] The updated temperature field and status labels are obtained, and the fluid property parameters of each changed unit are recalculated in combination with the permeability mutation data to form a fluid property update table.

[0068] The local pressure field is reconstructed based on the fluid property update table to obtain local pressure field data;

[0069] The pulse flux vector of the modified unit is calculated based on the local pressure field data to obtain the pulse flux vector table;

[0070] Chain reaction propagation simulation was performed based on the pulse flux vector table to obtain the chain reaction activation table;

[0071] The spatiotemporal features of pulse events are extracted based on the pulse flux vector table and the chain reaction activation table, and pulse event data is generated.

[0072] Based on the pulse event data and viscosity-temperature partitioned grid, the degree of influence of each pulse event on its surrounding cells is calculated to form extended pulse influence data;

[0073] The flow pulse event sequence is obtained by recording the extended pulse influence data and pulse event data.

[0074] In this embodiment of the invention, each state change unit is traversed, and its instantaneous permeability mutation factor F_k is calculated based on the transformation type of its internal grid cells (e.g., from stagnant to critical, or from stagnant to free). The permeability mutation factor F_k is calculated using a nonlinear function; for example, for a "stagnant to critical" transformation, F_k = exp(β_1 × ΔT_actual / ΔT_threshold), where β_1 is the mutation coefficient, ΔT_actual is the actual temperature rise, and ΔT_threshold is the width of the critical activation interval. For a "stagnant to free" transformation, F_k includes an additional strengthening factor K_boost. The old permeability k_old is multiplied by F_k to obtain the new permeability k_new. The calculation process considers lithological differences and tectonic effects; for example, F_k is further amplified near fault zones. All changed units and their new instantaneous permeability values ​​constitute permeability mutation data.

[0075] Obtain the updated temperature field and the latest state labels of the grid cells at the current time step. For each cell within a state-changing cell, recalculate the cell's fluid properties, including instantaneous fluid viscosity μ_fluid (calculated based on the new temperature and viscosity response function), density ρ_fluid, and specific heat capacity C_p_fluid, by combining the new permeability value k_new from the permeability mutation data with its current temperature T and pressure P. For example, when the cell state changes to "low viscosity free flow region," the μ_fluid value drops sharply, and the k_new value increases significantly. These updated fluid properties (k_new, μ_fluid, ρ_fluid, C_p_fluid), along with the cell's spatial coordinates and current temperature, form the fluid property update table.

[0076] Based on the fluid property update table, a local pressure field reconstruction is performed in and around the state-changing cell (e.g., within a radius of 3 grid cells). Because instantaneous and abrupt changes in permeability and viscosity cause significant changes in local fluid resistance, the pressure field adjusts rapidly. Using the new permeability and viscosity values ​​from the fluid property update table, a fast iterative solver (e.g., a multigrid method) is employed to solve the seepage control equations within the local region. The boundary conditions are provided by the ambient pressure of the surrounding stable region. The reconstructed result is a refined local pressure field data that reflects the instantaneous pressure gradient established within the region of abrupt permeability change, which is the direct driving force behind the "pulsating" flow.

[0077] After obtaining the local pressure field data, for each grid cell within a state-change unit, the instantaneous fluid flux vector q = -(k_new / μ_fluid)∇P is calculated based on Darcy's law. Here, k_new and μ_fluid are obtained from the fluid property update table, and ∇P is the pressure gradient calculated from the local pressure field data. The flux vector q has both direction and magnitude, precisely describing the direction and intensity of fluid flowing out of or into the cell at the instant of a sudden change. These instantaneous fluid flux vector values ​​(q_x, q_y, q_z) and their corresponding grid cell coordinates are recorded in the pulse flux vector table, quantifying the instantaneous fluid flow rate and direction caused by the viscosity change.

[0078] High-flux regions in the pulse flux vector table are used as "activation sources" to simulate the propagation of a "chain reaction" to surrounding regions. A chain reaction threshold is set (e.g., fluid flux exceeding 10). -7 (m / s) identifies all grid cells exceeding this threshold. For these cells, the temperature perturbation ∆T_perturbation on neighboring cells in the "critical activation region" but without state change is calculated. If the ∆T_perturbation is sufficient to cause the temperature of a neighboring cell to cross T_p or T_m, the neighboring cell is marked as "chain reaction activation". This process iterates within a finite propagation distance (e.g., 5 grid cells) until no new cells are activated. The chain reaction activation table records all activated cells and their activation time and intensity, depicting the spatial spread and cascading effects of the impulse event.

[0079] By integrating the pulse flux vector table and the chain reaction activation table, the complete spatiotemporal characteristics of each independent pulse event are extracted. A pulse event is defined as a spatially continuous high-flux region triggered by one or more state-changing units and propagated through a chain reaction. For each identified pulse event, the timestamp, event center coordinates (e.g., coordinates of the maximum flux unit), pulse intensity (flux integral or average), event type, influence range (chain reaction activation region), and duration (e.g., 2×Δt) are recorded. This structured information constitutes pulse event data, accurately depicting the spatiotemporal fingerprint of "pulsed" mineralization events.

[0080] We receive pulse event data and, combined with lithological and structural information stored in the viscosity-temperature zoning grid, calculate the long-term impact of each pulse event on its surrounding cells. This impact is quantified using the influence function I(d,θ,lith), where d is the distance, θ is the angle with the main pulse direction, and lith is the lithology. The influence function considers distance attenuation, structural orientation (slow attenuation along fault zones), and lithological permeability. The I(d,θ,lith) value is multiplied by the pulse event intensity to obtain the cumulative flux contribution potential of that pulse event to the surrounding cells. The cumulative flux contribution potential values ​​of all cells under each pulse event constitute the extended pulse impact data.

[0081] The core spatiotemporal features of the pulse event data, along with the cumulative flux contribution potential of surrounding cells calculated from the extended pulse impact data, are integrated into a structured record. Each record represents a complete flow pulse event, including a timestamp, center coordinates, pulse intensity, impact range, duration, and a cumulative impact weight matrix on the surrounding environment. These records are appended to the flow pulse event sequence in real time and indexed for rapid retrieval. At the end of the simulation, the flow pulse event sequence will contain all detected pulse events with spatiotemporal fingerprints, forming a complete dataset describing the "pulsating" migration history of geological fluids.

[0082] Preferably, step S3 includes the following steps:

[0083] Step S31: Parse the pulse event data of the flow pulse event sequence to obtain the index pulse events;

[0084] Step S32: Construct the initial flux evaluation grid using index pulse events;

[0085] Step S33: Traverse the index pulse events, calculate the impact value of each pulse time on each point in space based on the initial flux evaluation grid and the preset impact function, and form an event impact distribution set;

[0086] Step S34: Analyze the spatiotemporal relationship between different pulse events in the index pulse event based on the event impact distribution set, identify the existing event cascade effect, and obtain the event superposition enhancement factor;

[0087] Step S35: Superimpose the event superposition enhancement factor onto the initial flux evaluation grid to calculate the cumulative flux and obtain the original cumulative flux field;

[0088] Step S36: Perform flux threshold filtering on the original cumulative flux field to obtain the filtered flux field;

[0089] Step S37: Perform channel connectivity analysis on the screened flux field to obtain channel performance characteristics;

[0090] Step S38: Generate a cumulative transport flux field based on channel performance characteristics and the selected flux field.

[0091] In this embodiment of the invention, a simulated sequence of flowing pulse events is loaded and parsed. Each record contains an event ID, an occurrence timestamp (t_event), three-dimensional center coordinates (x_event, y_event, z_event), an instantaneous pulse intensity (I_pulse), and ellipsoidal parameters of the affected area. After parsing, a list of event objects residing in memory is constructed. To improve query efficiency, a B+ tree temporal index (with t_event as the key) and an R-tree spatial index (with the event's affected area bounding box as the key) are simultaneously established. These indexes support fast retrieval of pulse events in specific spatiotemporal regions. Finally, a set of pulse events containing the index structure is generated, facilitating subsequent spatiotemporal analysis.

[0092] Based on the geometry and resolution of the original geological model, a matching three-dimensional discrete evaluation grid is constructed. Each cell of this grid corresponds to a physical location in the geological model. The cumulative flux value of each cell in the initial flux evaluation grid is initialized to zero, while retaining its original coordinate information. This grid is specifically designed to accumulate and record the fluid flux contributed by all impulsive events throughout the simulation time history. For example, for a 100×100×100 original geological grid, the initial flux evaluation grid will also be a three-dimensional array of the same size, with all elements initialized to 0.0.

[0093] Iterate through each pulse event in the indexed pulse event set. Based on its center coordinates, pulse intensity, and influence range parameters, calculate the cumulative flux contribution of that event to all cells in the initial flux assessment grid. The calculation uses a pre-defined three-dimensional Gaussian decay influence function. The function has the highest intensity at the event center, decays with distance, and takes into account the anisotropy of geological structures (e.g., slower decay along fault zones). Influence value .in This represents the distance of the cell in the event's local coordinate system. The attenuation coefficient is along the axial direction and is affected by the event's influence range and structural orientation. All cells represent the influence distribution of the event, with the flux contribution values ​​generated by a single pulse event forming the event influence distribution set.

[0094] Spatiotemporal proximity analysis is performed on impulsive events with concentrated event impact distribution. First, an event spatiotemporal correlation matrix A is constructed, where element A_ij = exp(-Δt_ij / T_decay - Δs_ij / S_decay), representing the spatiotemporal correlation between events i and j (Δt_ij is the time interval, Δs_ij is the spatial distance, and T_decay and S_decay are decay constants). If A_ij exceeds a preset threshold (e.g., 0.5), events i and j are considered correlated. Next, connectivity analysis is performed on matrix A to identify event triggering chains composed of mutually correlated events (e.g., event A triggers event B). For each cell in the identified event triggering chain, the superposition enhancement factor F_enhancement = 1 + N_overlap × f_overlap(Δt_avg) is calculated. Here, N_overlap is the number of overlapping correlated events in that cell, and f_overlap(Δt_avg) is a function that indicates a stronger enhancement effect when the average time interval Δt_avg is smaller. This method captures the local flux accumulation enhancement caused by the spatiotemporal interaction of multiple impulsive events. The F_enhancement value of all cells forms the event superposition enhancement factor.

[0095] The process iterates through each cell in the initial flux assessment grid. For each cell, all impulse events affecting that cell are retrieved (determined by the event influence distribution set), and the flux contribution of these events at that cell location is summed. For cells with event cascade effects, an event superposition enhancement factor is applied to nonlinearly amplify the summed result. The cumulative flux C_total(i,j,k) = Σ(I_contrib_m(i,j,k) × F_enhancing_m(i,j,k)), where m represents the m-th impulse event affecting that cell, I_contrib_m is the flux contribution of that event, and F_enhancing_m is the corresponding superposition enhancement factor. This process accumulates and integrates the instantaneous flux and cascade effects generated by all impulse events at various points in space throughout the entire simulation history. The final result is a three-dimensional scalar field—the original cumulative flux field—where the value of each cell represents the total fluid flux experienced by that point in the simulation history, without filtering.

[0096] Statistical analysis was performed on the original cumulative flux field, calculating its mean (μ_total) and standard deviation (σ_total). A flux threshold T_cutoff = μ_total + 1.5 × σ_total was set. For each cell in the original cumulative flux field, if its cumulative flux value C_total(i,j,k) was lower than T_cutoff, its flux value was set to 0; if it was higher than or equal to T_cutoff, its original value was retained. This screening aimed to remove background noise and low-intensity random fluctuations, highlighting high-flux channels formed by significant impulse events and cascading effects. The stringency of the screening was controlled by adjusting a coefficient of 1.5 times the standard deviation. The screened results formed the screened flux field, mainly containing information on dominant fluid transport channels contributed by strong impulse events.

[0097] A three-dimensional connected component labeling algorithm was applied to the flux field to identify and extract continuous high-fluidity regions, defined as potential fluid transport channels. The algorithm started with seed points above a threshold and recursively expanded to all neighboring cells also above the threshold until an independent connected region was identified. For each identified channel, a series of performance characteristic parameters were calculated: channel volume, length (through skeletonization), average flux value, maximum flux value, flux variation coefficient, and principal direction (through principal component analysis). Furthermore, the start, end, and branch points of the channels were identified, as these locations are often favorable areas for mineral deposition. The performance characteristic parameters of all channels were recorded in a channel performance characteristic table.

[0098] The final cumulative transport flux field is generated by integrating the selected flux field and channel efficiency characteristic tables. In this step, the selected flux field is optimized and enhanced based on the channel efficiency characteristics. For example, the flux values ​​of identified main channels can be nonlinearly enhanced to highlight their importance in visualization. For complex channel networks with multiple branches and confluences, integration is performed based on their overall connectivity and intensity. Finally, all processed flux values ​​are reassembled into a three-dimensional scalar field. This cumulative transport flux field not only quantifies the total intensity of fluid flow at various points in space throughout the simulation timescale but also clearly delineates the dominant fluid transport channel network driven by impulse events, providing the most direct physical basis for subsequent mineralization potential assessment.

[0099] Preferably, step S34 includes:

[0100] Construct a spatiotemporal correlation matrix of events based on the event impact distribution set, and then identify related event groups;

[0101] Perform event trigger chain analysis on the associated event groups to obtain event trigger chain data;

[0102] Identify the resonance enhancement region in the event trigger chain data;

[0103] Calculate the event superposition enhancement factor according to the resonance enhancement region.

[0104] In the embodiments of the present invention, all pulse events in the concentrated event influence distribution are traversed. For any two events i and j, calculate the spatio-temporal correlation intensity A_ij = exp(-(Δt_ij / T_decay)-(Δs_ij / S_decay)). Δt_ij is the time interval between event occurrences, Δs_ij is the spatial distance, and T_decay and S_decay are the time and spatial decay constants. If A_ij is greater than the preset correlation threshold (such as 0.6), it is considered that there is a significant spatio-temporal correlation between events i and j. The A_ij values of all event pairs are constructed into a symmetric event spatio-temporal correlation matrix A. Subsequently, perform a connected component analysis on matrix A, and divide the mutually correlated events into multiple non-overlapping correlated event groups, where each group represents a cluster of pulse events that are closely connected in space and time.

[0105] After identifying the correlated event groups, analyze the events within each group to determine the causal triggering relationship. For events i and j within a group, if the occurrence time of i is earlier than that of j (t_event_i < t_event_j), and the influence range of i overlaps with the triggering region of j, and the intensity of i is sufficient to cause temperature and pressure fluctuations in the j region, then i is a potential trigger for j. Use a directed graph to represent the trigger chain, where the nodes are pulse events and the edges represent the triggering relationship, and the weights can be set according to the triggering intensity and time delay. By traversing the directed graph, identify all event trigger chains, describe the complete propagation path from the initial pulse to the cascade response, and form the event trigger chain data.

[0106] Analyze the event trigger chain data to identify the "resonance enhancement regions" where multiple trigger chains or multiple events of a single trigger chain act in the same region in space and time. These regions refer to the regions where the influences of multiple pulse events are superimposed at the same spatial position in a short period of time, resulting in non-linear enhancement of fluid flux and ore-forming material migration ability. The identification process includes: for each grid cell, count the number of different pulse events covering the cell within a given time window (such as 5 time steps), and evaluate the average trigger intensity and time interval of these events. If the number of covering events exceeds the threshold (such as 3), and the average time interval is less than the threshold (such as 2 time steps), then mark the cell as a resonance enhancement region. All identified resonance enhancement regions constitute the resonance enhancement region data.

[0107] For each grid cell, the event superposition enhancement factor F_superposition is calculated based on whether it is identified as a resonance enhancement region and the number of affected pulse events. If the cell is within a resonance enhancement region and is affected by N_c pulse events, then F_superposition = 1 + α × N_c + β × f_resonance(T_resonance, S_resonance). Here, α and β are weighting coefficients, N_c is the number of covered events, and f_resonance is the temporal and spatial intensity factor based on the resonance region. For cells not within a resonance enhancement region, F_superposition = 1. This enhancement factor is a value greater than or equal to 1, used to nonlinearly amplify the cumulative effect of fluid flux, thereby quantifying the spatiotemporal coupling relationship and cumulative enhancement effect between multiple events, providing key correction parameters for subsequent cumulative flux calculations.

[0108] Preferably, step S4, which involves performing a fluid flux-ore-forming condition matching analysis on the cumulative transport flux field, includes:

[0109] Temperature and pressure data are extracted from the viscosity-temperature partitioned grid, and then the ore-forming condition template parameters are defined.

[0110] Calculate the rate of change of temperature and pressure data;

[0111] The cooling efficiency index is calculated based on the temperature and pressure change rate data and the cumulative transport flux field.

[0112] The effective residence time index is obtained by evaluating the fluid residence time of the cumulative transport flux field.

[0113] The matching degree of ore-forming conditions is calculated based on the cooling efficiency index, effective residence time index, and ore-forming condition template parameters to obtain the ore-forming condition matching degree field.

[0114] In this embodiment of the invention, the temperature (T) and pressure (P) change history of each grid cell throughout the entire simulation time is extracted from the viscosity-temperature partitioned grid, including the final T_final and P_final values. Simultaneously, based on the target mineral type (e.g., gold ore), a set of mineralization condition template parameters is defined, including ideal mineralization temperature and pressure ranges [T_min_ore, T_max_ore] and [P_min_ore, P_max_ore], as well as ideal ranges for processes such as fluid cooling rate and residence time, serving as a benchmark for evaluating mineralization potential.

[0115] The temperature and pressure change trajectories of each grid cell are analyzed, and the average and instantaneous temperature-time change rate (dT / dt), temperature gradient (∇T), pressure-time change rate (dP / dt), and pressure gradient (∇P) are calculated. Special attention is paid to regions of abrupt temperature and pressure drops, as these are conducive to the supersaturation and precipitation of minerals. These calculated temperature and pressure change rate data are stored in a three-dimensional array, recording the average and maximum rates of change.

[0116] Combining temperature and pressure change rate data with the cumulative transport flux field, the cooling efficiency index (CE_cool) for each grid cell is calculated. This index measures the product of the rate of temperature decrease and the cumulative flux as the fluid passes through a region. CE_cool = C_total × (dT / dt_avg) × (T_final_initial_diff / T_max_ore_min_ore_diff). Where C_total is the cumulative transport flux, dT / dt_avg is the average cooling rate, T_final_initial_diff is the absolute value of the initial and final temperature difference, and T_max_ore_min_ore_diff is the ideal temperature difference range for the mineralization template. A high CE_cool value indicates that the region is conducive to mineral precipitation caused by rapid cooling. The CE_cool values ​​of all cells constitute the cooling efficiency index field.

[0117] The fluid residence time is evaluated for each grid cell based on the cumulative transport flux field, yielding the effective residence time index (TE_res). TE_res = V_cell / (C_total / T_total), where V_cell is the cell volume and T_total is the total simulation time. This index reflects the average time a unit volume of fluid spends within the cell. The TE_res values ​​are then matched against the ideal residence time interval [TE_min_ore, TE_max_ore] defined in the mineralization condition template. Cells with a high degree of matching receive higher TE_res values. The effective residence time indices of all cells constitute the effective residence time index field.

[0118] The cooling efficiency index, effective residence time index, and ore-forming condition template parameters are weighted and fused to calculate the comprehensive ore-forming condition matching degree for each grid cell. M_ore = w_1 × normalized(CE_cool) + w_2 × normalized(TE_res) + w_3 × normalized(T_match) + w_4 × normalized(P_match). Here, w_1, w_2, w_3, and w_4 are weighting coefficients (e.g., 0.4, 0.3, 0.2, 0.1, which can be adjusted according to the mineral type), and each term is a normalized index or a function of the degree of matching with the ideal temperature and pressure range of the ore-forming template. A higher M_ore value indicates that the physicochemical conditions of that cell are closer to the ideal ore-forming environment of the target mineral. The M_ore values ​​of all cells constitute a three-dimensional ore-forming condition matching degree field.

[0119] Preferably, step S4, which involves delineating and visualizing the mineralization potential target area based on the mineralization condition matching field, includes:

[0120] A flux distribution characteristic analysis was performed on the cumulative transport flux field to obtain a flux distribution characteristic report;

[0121] Potential classification thresholds for different mineralization potential levels were determined based on flux distribution characteristics analysis.

[0122] The cumulative transport flux field and the mineralization condition matching degree field are weighted and fused to calculate the comprehensive mineralization potential field;

[0123] Spatial connectivity analysis of the comprehensive mineralization potential field is performed based on the potential grading threshold to obtain the spatial entity of the target area;

[0124] The target area spatial entity is optimized by performing target area boundary optimization to obtain the optimized target area entity;

[0125] The optimized target area entities are classified and sorted to obtain graded target area data.

[0126] Uncertainty assessment is performed on the graded target area data to obtain a target area reliability assessment report;

[0127] By integrating comprehensive mineralization potential field, graded target area data, and target area reliability assessment reports, a three-dimensional mineralization probability cloud map is generated.

[0128] In this embodiment of the invention, statistical analysis is performed on the cumulative transport flux field to construct histograms and cumulative frequency curves. Kernel density estimation is applied to identify distribution inflection points as natural boundaries of flux levels. Simultaneously, spatial autocorrelation of the flux field (such as the Moran index) is calculated to determine the clustering effect in high-flux regions. These analytical results, including mean, median, standard deviation, skewness, kurtosis, natural boundaries, and spatial autocorrelation reports, constitute a flux distribution characteristic report.

[0129] Based on flux distribution characteristic reports, a set of multi-level thresholds is automatically determined to classify mineralized areas into four potential levels: extremely high, high, medium, and low. Clustering-based methods (such as K-means or Gaussian mixture models) are used to classify flux value distributions; for example, the highest level threshold is set as a flux field value P95, the second highest as P80, and the medium level as P60. An interactive interface is provided for geological experts to fine-tune and calibrate. Finally, four distinct potential classification thresholds, T_rank1, T_rank2, T_rank3, and T_rank4, are determined, forming a potential classification threshold table.

[0130] The cumulative transport flux field and the mineralization condition matching degree field are weighted and fused to generate a comprehensive mineralization potential field. The comprehensive mineralization potential value of each cell is P_total = w_flow × P_flow_norm + w_match × P_match_norm. Here, P_flow_norm and P_match_norm are the normalized flux field and matching degree field values ​​(range [0,1]), respectively, and w_flow and w_match are weighting coefficients, which sum to 1. For example, w_flow can be set to 0.6 and w_match to 0.4, and the weighting coefficients are adjusted according to the mineralization mechanism of the target mineral. The fusion result is a new three-dimensional scalar field, and the cell values ​​represent the normalized comprehensive mineralization potential.

[0131] The comprehensive mineralization potential field is segmented using a potential grading threshold. Grid cells with potential values ​​higher than the highest grading threshold T_rank1 are identified, and a three-dimensional connected component labeling algorithm (such as a flood-fill algorithm based on BFS or DFS) is used to aggregate spatially adjacent high-potential cells into independent, continuous three-dimensional entities, i.e., target area spatial entities. For each target area, its volume, shape factor, surface area, average potential value, and maximum potential value are calculated. This step ensures the target area has geological continuity. All identified target area spatial entities and their geometric and potential characteristic data constitute the target area spatial entity set.

[0132] The initially identified target area spatial entities undergo boundary smoothing and morphological optimization. For jagged boundaries caused by mesh discretization, adaptive surface smoothing algorithms (such as Laplacian smoothing or Taubin smoothing) are used to process the 3D boundaries, making their shapes more consistent with the characteristics of natural geological bodies (e.g., smoother or extending along tectonic trends). Geological structural information is considered during the smoothing process; for example, boundaries along fault zones have lower smoothness. Simultaneously, the possibility of merging adjacent but low-potential cells with high total potential values ​​and close proximity is evaluated. Finally, optimized target area entities with smooth and reasonable geological morphology are generated.

[0133] All optimized target areas are comprehensively evaluated, classified, and ranked. Evaluation indicators include target area volume, average comprehensive potential value, maximum comprehensive potential value, shape regularity, and spatial relationship with known mineral deposits.

[0134] Calculate a comprehensive score Score_target for each optimized target region:

[0135] Score_target = w_vol × V_norm + w_avgP × P_avg_norm + w_maxP × P_max_norm + w_shape × Shape_norm, where w is the weight and each term is a normalization parameter. Based on the comprehensive score, the target area is divided into three levels: Level 1 target area (score higher than 0.8), Level 2 target area (score between 0.6 and 0.8), and Level 3 target area (score between 0.4 and 0.6). The data for each level of target area includes detailed information such as level, score, geometric attributes, and internal potential distribution statistics, forming structured leveled target area data.

[0136] Uncertainty assessment was performed using Monte Carlo simulations of graded target area data. Multiple (e.g., 100) complete simulations were conducted using key input parameters of a stochastic perturbation model (such as initial temperature field, lithological permeability, and viscosity-temperature response function parameters). The frequency with which each grid cell was identified as a high-potential target area (Level 1 or Level 2) in the 100 simulations was statistically analyzed, and a reliability score was assigned accordingly (High: >80%, Medium: 50%-80%, Low: <50%). For each optimized target area entity, the average reliability score of its internal cells was calculated. The assessment results were integrated into a target area reliability assessment report, quantifying the stability and reliability of the prediction results.

[0137] By integrating comprehensive mineralization potential field data, graded target area data, and target area reliability assessment reports, a final 3D mineralization probability cloud map is generated. The cloud map employs volumetric rendering visualization technology, intuitively representing mineralization potential through pixel transparency and color mapping: potential values ​​are mapped to colors (e.g., a blue-to-red gradient), and reliability scores are mapped to transparency (high-reliability areas have high opacity). Primary target areas are highlighted with solid bounding boxes, while secondary and tertiary target areas are displayed with semi-transparent outlines. Users can interactively manipulate the map (rotate, zoom, section) and query information for any point. The cloud map supports overlaying of geological structures, lithology, and other background information, providing comprehensive 3D mineralization prediction results to guide exploration decisions.

[0138] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0139] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A three-dimensional simulation method for ore-forming fluids based on multi-field coupling effects, characterized in that, Includes the following steps: Step S1: Obtain and calculate the critical temperature threshold based on viscosity-temperature response for the grid cells of the three-dimensional geological grid model data, and then classify the grid cell states to obtain grid state labeling data; Geological grid properties are extended based on grid state labeling data to obtain viscosity-temperature partitioned grids; Step S2: Construct a simulation environment using a viscosity-temperature partitioned grid and extract temperature thresholds to form comprehensive temperature threshold data; State change detection is performed on the comprehensive temperature threshold data to obtain state change units; pulse events are analyzed and recorded in the cells of the viscosity-temperature partition grid where state changes occur based on the state change units to obtain the flow pulse event sequence; Step S3: Construct an initial flux assessment grid based on the pulse event data, identify the event cascade effect, and obtain the event superposition enhancement factor; The event superposition enhancement factor is superimposed on the cumulative quantization of the fluid transport channel efficiency within the simulation time history of the initial flux evaluation grid to form a cumulative transport flux field; Step S4: Perform fluid flux and ore-forming condition matching analysis on the cumulative transport flux field to obtain the ore-forming condition matching degree field; Based on the matching degree field of mineralization conditions, the target area of ​​mineralization potential is delineated and visualized to obtain a three-dimensional mineralization probability cloud map.

2. The three-dimensional simulation method for ore-forming fluids based on multi-field coupling effect according to claim 1, characterized in that, Step S1 includes the following steps: Step S11: Analyze the three-dimensional geological grid model data to generate basic geological grid data; Step S12: Based on different lithological types and pressure conditions in the geological grid base data, define viscosity-temperature response curves to form viscosity response functions; Step S13: Determine the critical temperature threshold for the viscosity response function to obtain the critical temperature threshold; Step S14: Traverse the cells of the geological grid basic data, classify the grid cell status according to the critical temperature threshold, and obtain grid status label data; Step S15: Extend the geological grid base data with viscosity-temperature properties using grid state label data and viscosity response function to obtain a viscosity-temperature partitioned grid.

3. The three-dimensional simulation method for ore-forming fluids based on multi-field coupling effect according to claim 2, characterized in that, Step S14 includes: Read the temperature field data and lithological distribution data from the geological grid base data, divide the basic state according to the critical temperature threshold, and obtain the preliminary state classification results; Calculate the degree of state transition of the unit distance in the preliminary state classification results to form critical region sensitivity data; Analyze the neighborhood state structure of the units in the preliminary state classification results to form a neighborhood state structure; State stability is assessed based on neighborhood state structure and critical region sensitivity data to obtain a state stability index. By integrating critical region sensitivity data and state stability index, grid state labeling data is generated.

4. The three-dimensional simulation method for ore-forming fluids based on multi-field coupling effect according to claim 1, characterized in that, In step S2, a simulation environment is constructed using a viscosity-temperature partitioned grid, and temperature thresholds are extracted, including: Using the viscosity-temperature partitioned grid as the initial state, the spatiotemporal iterative framework is obtained and initialized according to the preset simulation parameters to obtain the initial simulation environment; The temperature field is dynamically updated based on the initial simulation environment and viscosity-temperature partitioned grid to obtain the updated temperature field; Temperature time-series comparison data are extracted from the updated temperature field and combined with lithological information and viscosity critical thresholds extracted from the viscosity-temperature zoning grid to form comprehensive temperature threshold data.

5. The three-dimensional simulation method for ore-forming fluids based on multi-field coupling effect according to claim 1, characterized in that, Step S2 involves detecting state changes in the comprehensive temperature threshold data, including: Determine whether temperature changes in the comprehensive temperature threshold data have led to changes in the status label, and obtain candidate change units; Candidate change units are classified by transformation type to obtain classified change units; A quantitative assessment of the change intensity of the classification change units is performed to obtain an intensity assessment change table; The intensity assessment change table is divided into change clusters and state change units are generated.

6. The three-dimensional simulation method for ore-forming fluids based on multi-field coupling effect according to claim 1, characterized in that, Step S2 involves analyzing and recording impulse events in cells of the viscosity-temperature grid where state changes occur, based on the state change units. Based on the state change unit, the permeability mutation of the viscosity-temperature partition grid is calculated to obtain the permeability mutation data; The updated temperature field and status labels are obtained, and the fluid property parameters of each changed unit are recalculated in combination with the permeability mutation data to form a fluid property update table. The local pressure field is reconstructed based on the fluid property update table to obtain local pressure field data; The pulse flux vector of the modified unit is calculated based on the local pressure field data to obtain the pulse flux vector table; Chain reaction propagation simulation was performed based on the pulse flux vector table to obtain the chain reaction activation table; The spatiotemporal features of pulse events are extracted based on the pulse flux vector table and the chain reaction activation table, and pulse event data is generated. Based on the pulse event data and viscosity-temperature partitioned grid, the degree of influence of each pulse event on its surrounding cells is calculated to form extended pulse influence data; The flow pulse event sequence is obtained by recording the extended pulse influence data and pulse event data.

7. The three-dimensional simulation method for ore-forming fluids based on multi-field coupling effect according to claim 1, characterized in that, Step S3 includes the following steps: Step S31: Parse the pulse event data of the flow pulse event sequence to obtain the index pulse events; Step S32: Construct the initial flux evaluation grid using index pulse events; Step S33: Traverse the index pulse events, calculate the impact value of each pulse time on each point in space based on the initial flux evaluation grid and the preset impact function, and form an event impact distribution set; Step S34: Analyze the spatiotemporal relationship between different pulse events in the index pulse event based on the event impact distribution set, identify the existing event cascade effect, and obtain the event superposition enhancement factor; Step S35: Superimpose the event superposition enhancement factor onto the initial flux evaluation grid to calculate the cumulative flux and obtain the original cumulative flux field; Step S36: Perform flux threshold filtering on the original cumulative flux field to obtain the filtered flux field; Step S37: Perform channel connectivity analysis on the screened flux field to obtain channel performance characteristics; Step S38: Generate a cumulative transport flux field based on channel performance characteristics and the selected flux field.

8. The three-dimensional simulation method for ore-forming fluids based on multi-field coupling effect according to claim 7, characterized in that, Step S34 includes: Construct a spatiotemporal correlation matrix of events based on the event impact distribution set, and then identify related event groups; Perform event trigger chain analysis on the associated event groups to obtain event trigger chain data; Identify the resonant enhancement regions in the event trigger chain data; The event superposition enhancement factor is calculated based on the resonance enhancement region.

9. The three-dimensional simulation method for ore-forming fluids based on multi-field coupling effect according to claim 1, characterized in that, Step S4 involves performing a fluid flux-ore-forming condition matching analysis on the cumulative transport flux field, including: Temperature and pressure data are extracted from the viscosity-temperature partitioned grid, and then the ore-forming condition template parameters are defined. Calculate the rate of change of temperature and pressure data; The cooling efficiency index is calculated based on the temperature and pressure change rate data and the cumulative transport flux field. The effective residence time index is obtained by evaluating the fluid residence time of the cumulative transport flux field. The matching degree of ore-forming conditions is calculated based on the cooling efficiency index, effective residence time index, and ore-forming condition template parameters to obtain the ore-forming condition matching degree field.

10. The three-dimensional simulation method for ore-forming fluids based on multi-field coupling effect according to claim 1, characterized in that, Step S4, which involves delineating and visualizing the mineralization potential target area based on the mineralization condition matching field, includes: A flux distribution characteristic analysis was performed on the cumulative transport flux field to obtain a flux distribution characteristic report; Potential classification thresholds for different mineralization potential levels were determined based on flux distribution characteristics analysis. The cumulative transport flux field and the mineralization condition matching degree field are weighted and fused to calculate the comprehensive mineralization potential field; Spatial connectivity analysis of the comprehensive mineralization potential field is performed based on the potential grading threshold to obtain the spatial entity of the target area; The target area spatial entity is optimized by performing target area boundary optimization to obtain the optimized target area entity; The optimized target area entities are classified and sorted to obtain graded target area data. Uncertainty assessment is performed on the graded target area data to obtain a target area reliability assessment report; By integrating comprehensive mineralization potential field, graded target area data, and target area reliability assessment reports, a three-dimensional mineralization probability cloud map is generated.

Images