A method for constructing a kinetic model of carbon dioxide mineralization and storage reaction considering optimization algorithms

Abstract

This invention provides a method for constructing a kinetic model of carbon dioxide mineralization and storage reaction considering optimization algorithms, belonging to the field of mining technology. This invention obtains the mineral composition of reservoir cores and constructs a graph database to constrain parameter boundaries. It designs a multi-temperature gradient batch reaction experiment combined with an active learning mechanism to optimize the sampling scheme. A hierarchical optimization strategy is used to gradually invert kinetic parameters from single-mineral to multi-mineral systems. A Gaussian process proxy model is constructed to accelerate the process using coarse and fine grid proxying. The uncertainty of parameters is quantified through Markov chain Monte Carlo sampling. Finally, the laboratory parameters are corrected to the reservoir scale and imported into a reservoir simulator. This invention solves the technical problem of simultaneously ensuring the convergence and identification accuracy of kinetic parameter inversion calculations in multi-mineral coupled systems.

Description

Technical Field

[0001] This invention belongs to the field of mining technology, and more specifically, it relates to a method for constructing a kinetic model of carbon dioxide mineralization and storage reaction that considers optimization algorithms. Background Technology

[0002] Geological carbon dioxide sequestration is an important technology for mitigating greenhouse gas emissions. Mineralization sequestration relies on the dissolution and precipitation reaction between carbon dioxide and reservoir minerals. Accurately describing the reaction kinetics of this process is fundamental to constructing reliable reservoir simulation models. Currently, researchers typically use batch reaction experiments combined with geochemical software (such as PHREEQC and TOUGHREACT) for forward simulations. They then manually adjust the rate constants and activation energies of minerals such as quartz, plagioclase, and calcite using trial-and-error methods or simple gradient optimization algorithms to fit the experimentally observed time-series data of mineral mass fractions and solution ion concentrations. This process is acceptable in single-mineral systems, but in actual sandstone or carbonate reservoirs, multiple minerals coexist and are coupled with each other. Dissolution products can promote or inhibit the precipitation of other minerals, leading to a rigid equation set. Traditional explicit solvers frequently diverge or have limited step sizes in long-scale simulations, resulting in extremely long forward simulation times. Meanwhile, the dimension of the parameter space increases exponentially with the number of mineral species. Global search algorithms suffer from computational explosion due to excessive calls to precise simulations, while local optimization algorithms get trapped in local extrema due to a lack of reasonable initial values. Both fail to converge to an accurate solution within a reasonable timeframe. In other words, existing technologies face the technical challenge of simultaneously ensuring convergence and identification accuracy in the inversion calculation of kinetic parameters under multi-mineral coupled systems. Summary of the Invention

[0003] In view of this, the present invention provides a method for constructing a kinetic model of carbon dioxide mineralization and storage reaction considering optimization algorithms, which can solve the technical problem in the prior art that it is difficult to simultaneously ensure the computational convergence and parameter identification accuracy of multi-mineral coupled systems in the process of inverting kinetic parameters of carbon dioxide mineralization and storage.

[0004] This invention is implemented as follows: This invention provides a method for constructing a carbon dioxide mineralization and sequestration reaction kinetic model considering optimization algorithms, including the following steps: obtaining the porosity and permeability of the target reservoir core; determining the content of minerals such as quartz, plagioclase, calcite, mica, and kaolinite through whole-rock mineral quantitative analysis using X-ray diffraction; extracting the parameter boundaries of mineral dissolution and precipitation reactions from publicly available data and literature; constructing a graph database; designing a multi-temperature gradient experiment for carbon dioxide mineralization and sequestration; introducing an active learning mechanism to recommend the next set of experimental conditions through Bayesian optimization; establishing a multi-temperature, multi-time-point mineral-solution time series database; inputting the mineral-solution time series database into a reaction network recognition model for training and prediction; and through reverse... A simplified core reaction network was formed based on contribution analysis, and a stoichiometric correlation matrix of mineral-solution components was established. Dissolution-precipitation rate equations were written in geochemical software, and an adaptive implicit time-step algorithm was introduced in conjunction with the Gill method to handle rigid ordinary differential equation sets. The kinetic parameters of minerals such as quartz, plagioclase, calcite, mica, and kaolinite were inverted through a hierarchical optimization strategy. A physical information-constrained surrogate model was constructed to accelerate the process using coarse and fine grid surrogates. Uncertainty quantification was performed using a normalized objective function based on multi-source data fusion combined with Bayesian inference from Markov chain Monte Carlo sampling. Laboratory-scale parameters were converted into reservoir-scale kinetic parameters and imported into a simulation model to customize the long-term mineralization and storage process.

[0005] Specifically, the graph database extracts the pre-exponential factor range, activation energy range, reaction order range, and pH index range of mineral dissolution and precipitation reactions from public data and literature, and constructs a triplet containing mineral name, reaction type, and parameter boundary to store in the graph database.

[0006] Specifically, the active learning mechanism involves using a Bayesian optimization algorithm to establish a Gaussian process model between existing experimental data and the uncertainty of dynamic parameters after each set of temperature and time point experiments is completed. The recommended values ​​for the next set of experimental conditions are calculated by maximizing the expected improvement criterion or the upper confidence bound criterion, and the temperature-time combination that can minimize the variance of the posterior distribution of parameters is preferentially selected.

[0007] Specifically, the carbon dioxide mineralization and storage multi-temperature gradient experiment involves setting a first temperature point, a second temperature point, and a third temperature point based on the pressure and temperature of the target reservoir core, with the pressure being higher than the critical pressure of carbon dioxide. A logarithmic time interval sampling method is used to take samples at the first, second, third, fourth, fifth, and sixth moments. The experimental water is simulated formation water, and the target reservoir core is ground into powder.

[0008] Specifically, the structure of the reaction network identification model is as follows: the input layer receives mineral content change data and solution ion concentration time series data, which are then passed through a pulse coding layer, a temporal memory layer, a dynamic feedback synapse layer, and a mineral coupling attention layer. The output layer outputs a simplified core reaction network topology, where nodes represent minerals and ions participating in the reaction, edges represent stoichiometric relationships, and edge weights represent the reaction contribution.

[0009] Specifically, the steps for establishing the training dataset for the reaction network identification model involve collecting no fewer than 50 sets of mineral-solution time series under different temperature gradients and pressure conditions as original samples, labeling the true core reaction network topology, dividing it into a training set and a test set at an 8:2 ratio, and performing data augmentation on the time series data in the training set.

[0010] Specifically, the hierarchical optimization strategy consists of three stages: first, single-temperature single-mineral optimization using a genetic algorithm or particle swarm optimization to search for the approximate range of parameters globally and constraining parameter boundaries using a graph database; second, single-temperature multi-mineral optimization using a quasi-Newton method or Bayesian optimization for local fine-tuning and inversion of Arrhenius parameters; and third, multi-temperature multi-mineral optimization using multi-task learning or co-evolutionary algorithms to gradually release mineral parameters.

[0011] Specifically, the adaptive implicit time step algorithm calculates the rate of change of all species concentrations at the beginning of each time step. When the ratio of the maximum rate of change to the minimum rate of change exceeds a rigid threshold, an implicit solver is selected; otherwise, an explicit solver is selected. The time step is dynamically adjusted based on the local truncation error estimate.

[0012] Specifically, the rigid threshold is obtained by selecting a typical... Numerical experiments were conducted on the water-rock reaction system, recording the statistical distribution of species concentration change rates over different time periods. The ratio of change rates that would limit the explicit solver step size to sub-second levels or cause divergence was set as a rigid threshold. Empirical values ​​for this rigid threshold were obtained through statistical analysis of 100 sets of core reaction simulation experiments with different mineral compositions. .

[0013] Specifically, the implementation of the Gill method in the scheme involves using a fourth-order Gill method to divide the time interval into four sub-intervals, calculating the slope estimate of the species concentration in each sub-interval, and using the weighted average of the four slope estimates as the concentration update for the entire time step. The weighting coefficients are determined based on the Butch table of the Gill method, and the slope values ​​of each sub-interval are obtained by iteratively solving the implicit equation system.

[0014] Specifically, the surrogate model for physical information constraints uses Gaussian process regression to establish a mapping relationship between input parameters and simulation output. During the optimization iteration process, when the parameter space exploration enters a low uncertainty region, the surrogate model for physical information constraints is called to predict the results. When it enters a high uncertainty region or the prediction error of the surrogate model for physical information constraints exceeds the prediction uncertainty threshold, it switches back to accurate forward simulation and updates the training set of the surrogate model for physical information constraints.

[0015] Specifically, the uncertainty quantification involves dividing a multi-temperature, multi-time-point mineral-solution time series database into a training set and a validation set using cross-validation. Parameters are inverted using the training set and the prediction capability is evaluated using the validation set. The generalization capability of the model is evaluated by using data from different core samples or different experimental conditions through independent experiments. The posterior probability distribution of the parameters is generated, and the confidence interval of the parameters is given.

[0016] Specifically, the reaction contribution analysis involves constructing a chemical network containing all reactions using the solution species database and phase database of geochemical software. Minerals with a mass change of less than 1% within 60 days or whose saturation index is always far from 1 are transferred from kinetic control to equilibrium control or ignored, forming a simplified core reaction network.

[0017] Specifically, the reservoir-scale kinetic parameters are output as a table of reservoir-scale kinetic parameters after correcting mineral surface area, effective reaction volume, and reactant transport limitation factor. A grid is constructed based on porosity and permeability, initial conditions, and rock and fluid physical properties.

[0018] Specifically, the mineralization and storage process involves adding mineral components such as quartz, plagioclase, calcite, mica, and kaolinite to the carbon capture and storage function to create a custom equilibrium reaction kinetic formula that conforms to the format of geochemical software. The reservoir-scale kinetic parameter table is then imported into the simulation model to analyze the sensitivity of injection rate, injection pressure, and temperature parameters.

[0019] Specifically, the steps for training the reaction network recognition model are as follows: the network weights are initialized using the Xavier initialization method, the learning rate is set to 0.001 and the Adam optimizer is used, the loss function is the weighted sum of the network topology prediction error and the reaction contribution prediction error, the weight coefficients are determined to be 0.6 and 0.4 through grid search, and an early stopping strategy is set to terminate training when the validation set loss does not decrease for 20 consecutive rounds.

[0020] This invention addresses the challenge of simultaneously ensuring convergence and identification accuracy in kinetic parameter inversion calculations within multi-mineral coupled systems by constructing a framework combining a hierarchical optimization strategy with a physical information constraint surrogate model. First, the invention utilizes a reaction network identification model to automatically identify core reaction paths and rank mineral importance from experimental time-series data. This allows the hierarchical optimization strategy to progressively tighten the parameter range from single-mineral to multi-mineral systems, avoiding gradient vanishing or diverging issues in high-dimensional parameter spaces. Simultaneously, an adaptive implicit time-step algorithm and the Gill method are introduced to handle rigid equation sets, fundamentally eliminating numerical divergence caused by multi-mineral coupling and ensuring convergence in each forward simulation. Furthermore, a Gaussian process regression surrogate model replaces the time-consuming precise simulation, directly predicting the objective function value in low-uncertainty regions, significantly reducing the number of precise simulation calls. Thus, while ensuring convergence, parameter identification accuracy is quantified through Bayesian inference. In summary, this invention solves the technical problem mentioned in the background art of simultaneously ensuring convergence and identification accuracy in kinetic parameter inversion calculations within multi-mineral coupled systems. Attached Figure Description

[0021] Figure 1 This is a flowchart of the method of the present invention.

[0022] Figure 2 This is a posterior probability distribution diagram of the activation energy parameters of each mineral.

[0023] Figure 3 for -Schematic diagram of the water-rock core reaction network topology. Detailed Implementation

[0024] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0025] like Figure 1 The diagram shown is a flowchart of a method for constructing a carbon dioxide mineralization and storage reaction kinetics model considering optimization algorithms, provided by this invention. This method includes the following steps:

[0026] S1. Obtain the porosity and permeability of the target reservoir core, determine the mineral content of quartz, plagioclase, calcite, mica, kaolinite and other minerals by whole-rock mineral quantitative analysis by X-ray diffraction, and extract the pre-exponential factor range, activation energy range, reaction order range and pH index range of mineral dissolution and precipitation reaction from public data and literature, and construct a ternary set containing mineral name, reaction type and parameter boundary and store it in the graph database.

[0027] S2. Based on the pressure and temperature of the target reservoir core, a multi-temperature gradient experiment for carbon dioxide mineralization and storage is designed. A first temperature point, a second temperature point, and a third temperature point are set, and the pressure is higher than the critical pressure of carbon dioxide. A logarithmic time interval sampling method is used to sample at the first, second, third, fourth, fifth, and sixth time points. The first time point is 1 day, the second time point is 3 days, the third time point is 7 days, the fourth time point is 14 days, the fifth time point is 30 days, and the sixth time point is 60 days. The experimental water is simulated formation water, and the target reservoir core is ground to powder. An active learning mechanism is introduced to recommend the next set of experimental conditions through Bayesian optimization to minimize the uncertainty of kinetic parameters.

[0028] S3. Perform mineral analysis on the reaction powder at the first temperature point, the second temperature point, the third temperature point, the first time point, the second time point, the third time point, the fourth time point, the fifth time point, and the sixth time point, and collect the corresponding aqueous solutions to determine the pH value, anion content, and cation content. Establish a multi-temperature and multi-time point mineral-solution time series database, and input the multi-temperature and multi-time point mineral-solution time series database into the reaction network recognition model for training and prediction.

[0029] S4. Using the solution species database and phase database of geochemical software, construct a chemical network containing all reactions. Through reaction contribution analysis, minerals with a mass change of less than 1% within 60 days or whose saturation index is always far from 1 are transferred from kinetic control to equilibrium control or ignored, forming a simplified core reaction network and establishing a stoichiometric correlation matrix of mineral-solution components.

[0030] S5. Based on the transition state theory, write the dissolution and precipitation rate equations for the corresponding minerals such as quartz content, plagioclase content, calcite content, mica content, and kaolinite content in the reaction rate data block of the geochemical software. Introduce an adaptive implicit time step algorithm in conjunction with the Gill method to process the rigid ordinary differential equation system, and dynamically switch the solver according to the species concentration change rate.

[0031] S6. The dynamic parameters of the corresponding minerals, such as the quartz content, plagioclase content, calcite content, mica content, and kaolinite content, are determined through a hierarchical optimization strategy. In the first stage, single-temperature single-mineral optimization is performed to fix the parameters of other minerals except the target mineral. Genetic algorithms or particle swarm optimization are used to search the approximate range of parameters in the global scope and the parameter boundaries are constrained by the graph database. In the second stage, single-temperature multi-mineral optimization is performed. The parameter results obtained in the first stage are used as initial values. Quasi-Newton method or Bayesian optimization is used to perform local fine-tuning and inversion of Arrhenius parameters. In the third stage, multi-temperature multi-mineral optimization is performed. Considering the coupling effect between minerals, multi-task learning or co-evolutionary algorithms are used to fix the secondary mineral parameters, optimize the main mineral parameters, and then gradually release more mineral parameters.

[0032] S7. Construct a proxy model with physical information constraints to replace part of the time-consuming precise forward simulation to achieve coarse and fine mesh acceleration. The proxy model with physical information constraints uses Gaussian process regression to establish the mapping relationship between input parameters and simulation output. During the optimization iteration process, when the parameter space exploration enters the low uncertainty region, the proxy model with physical information constraints is called to predict the results. When entering the high uncertainty region or the prediction error of the proxy model with physical information constraints exceeds the prediction uncertainty threshold, the model is switched back to the precise forward simulation and the training set of the proxy model with physical information constraints is updated.

[0033] S8. Uncertainty is quantified by using a normalized objective function based on multi-source data fusion combined with Bayesian inference based on Markov chain Monte Carlo sampling. The multi-temperature, multi-time-point mineral-solution time series database is divided into a training set and a validation set using cross-validation. Parameters are inverted using the training set and the prediction ability is evaluated using the validation set. The generalization ability of the model is evaluated by using data from different core samples or different experimental conditions through independent experiments. The posterior probability distribution of the parameters is generated and the confidence interval of the parameters is given.

[0034] S9. Convert laboratory-scale parameters into parameters suitable for reservoir scale, correct mineral surface area, effective reaction volume, and reactant transport limitation factor, and output a reservoir-scale kinetic parameter table. Construct a grid, initial conditions, and rock and fluid physical properties based on the porosity and permeability. Add the mineral composition corresponding to the quartz content, plagioclase content, calcite content, mica content, and kaolinite content to the carbon capture and storage function to create a custom equilibrium reaction kinetic formula conforming to the geochemical software format. Import the reservoir-scale kinetic parameter table into the simulation model to customize the simulation of a long-term mineralization and storage process and analyze the sensitivity of injection rate, injection pressure, and temperature parameters.

[0035] The implementation steps of the active learning mechanism are as follows: after completing each set of temperature and time point experiments, a Gaussian process model between the existing experimental data and the uncertainty of the dynamic parameters is established using the Bayesian optimization algorithm. The recommended values ​​of the next set of experimental conditions are calculated by maximizing the expected improvement criterion or the upper confidence bound criterion. The temperature-time combination that can minimize the variance of the posterior distribution of the parameters is selected first, so as to obtain the highest parameter identification accuracy with the fewest number of experiments.

[0036] The structure of the reaction network recognition model is as follows: The input layer receives mineral content change data and solution ion concentration time series data from the multi-temperature, multi-time-point mineral-solution time series database, with the data dimension being temperature number × time point number × species number; the first hidden layer is a pulse coding layer, which converts non-uniform sampling data with logarithmic time intervals into pulse sequences, activating spiking neurons only when the ion concentration gradient or mineral content fraction change rate exceeds a dynamic threshold, which is adaptively adjusted based on the concentration standard deviation within the previous time window; the second hidden layer is a temporal memory layer, using a long short-term memory network unit to process the pulse sequence. The long short-term memory network unit contains three gating mechanisms: a forget gate, an input gate, and an output gate. The forget gate determines whether to retain or discard historical information based on the current concentration change rate; the input gate updates the cell state based on the pulse intensity; and the output gate combines the cell state with the current... The first hidden layer generates a hidden state; the third hidden layer is a dynamic feedback synaptic layer, which feeds back the hidden state of the temporal memory layer to the pulse coding layer to adjust the excitation threshold at the next moment. The feedback weights are obtained through training using an error backpropagation algorithm. When the prediction error is greater than a set threshold, the excitation threshold is increased to improve the sensitivity to sudden phase transitions. When the prediction error is less than the set threshold, the excitation threshold is decreased to capture slow precipitation processes. The fourth hidden layer is a mineral coupling attention layer, which calculates the temporal correlation weights between different mineral dissolution and precipitation reactions through a multi-head attention mechanism. The rows and columns of the weight matrix correspond to the primary and secondary minerals, respectively, and the element values ​​represent the degree of influence of the secondary minerals on the reaction rate of the primary minerals. The output layer is a fully connected layer, which outputs the simplified core reaction network topology. Nodes represent the minerals and ions participating in the reaction, edges represent stoichiometric relationships, and edge weights represent the contribution of the reaction.

[0037] The steps for establishing the training dataset for the reaction network identification model include: collecting core experimental data from different reservoir types as raw samples, wherein the raw samples contain no fewer than 50 sets of mineral-solution time series under different temperature gradients and pressure conditions; labeling the true core reaction network topology for each set of raw samples, wherein the labeling method is to manually select the main reactions according to the reaction contribution threshold by geochemical experts and record the stoichiometric correlation matrix; dividing the raw samples into the training set and the test set in an 8:2 ratio, wherein the training set is used for model parameter learning and the test set is used to evaluate the generalization performance of the model; and performing data augmentation on the time series data in the training set, wherein the augmentation methods include adding Gaussian noise that conforms to the experimental measurement error distribution to the concentration data, randomly scaling the time axis to simulate different sampling frequencies, and randomly perturbing the mineral content to simulate sample heterogeneity.

[0038] The training steps of the reaction network recognition model include: initializing network weights using the Xavier initialization method to ensure consistent variance of activation values ​​in each layer, setting the learning rate to 0.001 and using the Adam optimizer for gradient descent; defining the loss function as the weighted sum of the network topology prediction error and the reaction contribution prediction error, where the network topology prediction error uses cross-entropy loss to measure the difference between predicted edges and true edges, and the reaction contribution prediction error uses mean square error to measure the deviation between predicted weights and true weights, with weight coefficients determined to be 0.6 and 0.4 through grid search; randomly sampling 32 groups of samples in each iteration to form a mini-batch input network for forward propagation to calculate the loss function value, and calculating the gradient of each layer's weights and updating the weight parameters using the backpropagation algorithm; setting an early stopping strategy, terminating training when the validation set loss does not decrease for 20 consecutive rounds to prevent overfitting; evaluating model performance on the test set after training, with evaluation metrics including topology prediction accuracy, root mean square error of reaction contribution prediction, and core reaction recognition recall.

[0039] The technical effects of the reaction network identification model in the scheme are as follows: Through the processing of non-uniform logarithmic time sampling by the pulse coding layer, the reaction network identification model can automatically identify the concentration change characteristics of the rapid dissolution stage and the slow precipitation stage, avoiding the omission of key phase transition information by traditional uniform sampling methods; the gating mechanism of the temporal memory layer realizes the modeling of long-term span dependencies, accurately capturing the promoting or inhibiting effect of mineral dissolution products on subsequent precipitation reactions; the dynamic feedback synaptic layer adaptively adjusts the excitation threshold according to the prediction error, enabling the reaction network identification model to have self-calibration capability when facing different reservoir types or experimental conditions, significantly improving the robustness of the core reaction network identification; the mineral coupling attention layer quantifies the complex cooperative and competitive relationships in the multi-mineral system, providing an accurate mineral importance ranking basis for the subsequent hierarchical optimization strategy, thereby greatly reducing the search space for the inversion of kinetic parameters.

[0040] The implementation steps of the adaptive implicit time step algorithm are as follows: at the beginning of each time step, the rate of change of the concentration of all species is calculated. The rigidity of the equation system is judged based on the order of magnitude span of the rate of change of the species concentration. When the ratio of the maximum rate of change to the minimum rate of change exceeds the rigidity threshold, the implicit solver is selected; otherwise, the explicit solver is selected. The implicit solver uses a backward difference scheme to discretize the ordinary differential equations and solves the nonlinear algebraic equation system through the Newton-Raphson iterative method. The convergence criterion for the iteration is that the relative residual is less than the tolerance value. The time step is dynamically adjusted according to the local truncation error estimate. When the local truncation error is less than the target accuracy, the time step is increased to accelerate the calculation. When the local truncation error exceeds the target accuracy, the time step is decreased to ensure accuracy. The step adjustment factor is calculated by the error controller.

[0041] The rigid threshold is obtained by selecting a typical... Numerical experiments were conducted on the water-rock reaction system, recording the statistical distribution of the species concentration change rate over different time periods. The ratio of change rates that would cause the explicit solver step size to be limited to sub-second levels or diverge was set as the rigid threshold. Through 100 sets of core reaction simulation experiments with different mineral compositions, the empirical value of the rigid threshold was statistically obtained. .

[0042] The principle of the Gill method is as follows: The Gill method is a class of high-order implicit Runge-Kutta methods. It improves numerical stability by setting multiple intermediate points within the time step and constructing a linear combination scheme. The rigidity decay factor of the Gill method decreases with increasing order, thereby enabling stable solution of rigid equations with a larger step size.

[0043] The implementation of the Gill method in the scheme is as follows: the time interval is divided into four sub-intervals using the fourth-order Gill method, the slope estimate of species concentration is calculated in each sub-interval, the weighted average of the four slope estimates is used as the concentration update for the entire time step, the weight coefficients are determined according to the Butch table of the Gill method, and the slope value of each sub-interval is obtained by iteratively solving the implicit equation system.

[0044] The technical advantages of the Gill method are as follows: compared with the traditional explicit Euler method or Runge-Kutta method, the stability region of the Gill method is significantly expanded, allowing the use of time steps tens of times larger than those of explicit methods without causing numerical oscillations or divergences. This reduces the time required for a single forward simulation from several hours to several minutes while ensuring computational accuracy. The high-order accuracy characteristic of the Gill method causes the local truncation error to decrease with the fourth power of the time step, significantly reducing the total number of time steps required for the same accuracy requirements and further reducing computational costs.

[0045] The calculation formula for the reaction contribution analysis is expressed as follows:

[0046] ;

[0047] in minerals Ions in dissolution and precipitation reactions stoichiometric coefficients The unit of reaction rate is , The reference value for the stoichiometric coefficient is taken as 1. Take as a reference value for the reaction rate Molecules represent reactions Para ion The contribution of concentration change, with the denominator representing the total reaction response of ions. The total contribution of concentration change.

[0048] The steps for establishing the stoichiometric correlation matrix are as follows: arrange all participating minerals in rows, arrange all ions and gaseous components in the solution in columns, and the matrix elements are the stoichiometric coefficients of the corresponding components in the dissolution and precipitation reactions of the corresponding minerals. The components generated in the dissolution reaction are given positive values, the components consumed are given negative values, the precipitation reaction is the opposite, and the elements that do not participate in the reaction are given zero values. The number of independent reactions is determined by null space analysis of the stoichiometric correlation matrix, and the major and minor reactions are identified by singular value decomposition.

[0049] The equation for the dissolution-precipitation rate is expressed as follows:

[0050] ;

[0051] in The unit of reaction rate is , The rate constant is in units of , The unit for mineral specific surface area is , The saturation index is dimensionless. and The reaction order is an empirical parameter, dimensionless, and is usually taken as 1. Take the reaction rate reference value as , Take the reference value of the rate constant as , Take the reference value for the specific surface area of ​​minerals as .

[0052] The temperature dependence of the rate constant is expressed as follows:

[0053] ;

[0054] in The rate constant at 25℃ is in units of , The activation energy unit is , The gas constant is taken as 8.314. , The unit of absolute temperature is , The reference value for the rate constant at 25℃ is taken as follows. , The reference value for the temperature correction term is set to 1. .

[0055] The formula for calculating the saturation index is as follows:

[0056] ;

[0057] in It is the ion activity product. is the solubility product constant, and both are dimensionless.

[0058] The first stage of the hierarchical optimization strategy, the single-temperature single-mineral optimization steps, are as follows: The kinetic parameters of all minerals except the target mineral are fixed to the literature-recommended values. Only the 25℃ rate constant and activation energy of the target mineral are optimized. The objective function is the root mean square error between the experimentally measured mass fraction and the simulated predicted mass fraction of the target mineral at a single temperature. The genetic algorithm is used for global search, with the population size set to 50, the crossover probability set to 0.8, the mutation probability set to 0.1, and the number of generations set to 100. The parameter search range fluctuates by an order of magnitude above and below the literature values ​​provided by the graph database.

[0059] The graph database constrains the parameter boundaries by querying all triples of the same type as the target mineral in the graph database, extracting the upper and lower bounds of their dynamic parameter ranges, and restricting the search space of the genetic algorithm to between 1.5 times the upper bound and 0.5 times the lower bound. If the individual parameters generated during the search process exceed the parameter boundaries, they are reset to the parameter boundary values.

[0060] The second stage of the hierarchical optimization strategy, specifically the single-temperature multi-mineral optimization, involves the following steps: using the mineral parameters obtained in the first stage as initial values, simultaneously optimizing the 25℃ rate constant and activation energy of all major minerals. The objective function is the weighted root mean square error between the experimentally measured mass fraction and the simulated predicted mass fraction of all major minerals at a single temperature. The weights are determined based on the mineral content, with higher-content minerals receiving greater weights. The quasi-Newton method is used for local fine-tuning. This method calculates the search direction using an approximate update formula for the Hessian matrix. In each iteration, the gradient of the objective function is calculated, and the parameter values ​​are updated. The iteration terminates when the gradient norm is less than a certain value. Or the number of iterations exceeds 200.

[0061] The Arrhenius parameter is the pre-exponential factor and the activation energy.

[0062] The steps of the third stage of multi-temperature, multi-mineral optimization in the hierarchical optimization strategy are as follows: Experimental data at all temperatures are simultaneously incorporated into the optimization; the objective function is the global root mean square error (RMSE) between the experimental measurements and simulated predictions of all minerals at all temperatures; the multi-task learning framework treats different temperatures as different tasks, achieving cross-temperature information transfer by sharing underlying parameters; the co-evolutionary algorithm divides the minerals into two subpopulations: a primary mineral group and a secondary mineral group; first, the secondary mineral parameters are fixed while optimizing the primary mineral parameters until convergence; then, the primary mineral parameters are fixed while optimizing the secondary mineral parameters until convergence; the two subpopulations evolve alternately until the global root mean square error no longer decreases; during the optimization process, elemental mass conservation is checked using the stoichiometric correlation matrix; if the total elemental deviation after a certain iteration exceeds 0.1%, the parameter combination is rejected and resampling is performed.

[0063] The criteria for distinguishing between primary and secondary minerals are as follows: based on the weight matrix of the mineral coupling attention layer output by the reaction network recognition model, the weight of each mineral as a row vector is calculated, and the minerals with the top 30% of the weight sums are defined as primary minerals, and the rest are secondary minerals.

[0064] The steps for constructing the surrogate model constrained by physical information are as follows: Gaussian process regression is selected as the surrogate model framework, the input is the dynamic parameter vector, and the output is the mineral content fraction time series obtained from the geochemical software simulation; initial training samples are generated in the parameter space through Latin hypercube sampling, the number of initial training samples being 10 times the parameter dimension; the exact forward simulation is run on each initial training sample to obtain the corresponding output result; the Gaussian process regression is constructed using a radial basis function kernel, where the length scale parameter and signal variance parameter of the radial basis function kernel are determined by maximizing the marginal likelihood function; during the optimization iteration process, whenever the genetic algorithm or the particle swarm optimization algorithm generates a new parameter combination, the surrogate model constrained by physical information is first used to predict the corresponding output and calculate the prediction uncertainty. If the prediction uncertainty is less than the prediction uncertainty threshold, the prediction result is directly used to evaluate the fitness; otherwise, the exact forward simulation is called to obtain the true result, and the parameter combination is added to the training set of the surrogate model constrained by physical information for retraining.

[0065] The prediction uncertainty threshold is obtained as follows: leave-one-out cross-validation is performed on the initial training sample set, the root mean square error between the predicted value and the true value of the surrogate model constrained by the physical information is calculated, and the root mean square error is used as the initial value of the prediction uncertainty threshold; during the optimization process, the cross-validation error is recalculated and the prediction uncertainty threshold is updated after every 10 new samples are added, so that the prediction uncertainty threshold gradually decreases as the training set of the surrogate model constrained by the physical information is expanded.

[0066] The implementation of the coarse-fine mesh proxy acceleration is as follows: the parameter space is divided into a coarse mesh region and a fine mesh region. The coarse mesh region corresponds to the large-scale search stage of the parameters. At this time, the surrogate model constrained by physical information uses fewer training samples and a larger length scale parameter of the radial basis function to quickly provide approximate predictions. The fine mesh region corresponds to the local fine-tuning stage of the parameters. At this time, the surrogate model constrained by physical information uses more training samples and a smaller length scale parameter of the radial basis function to provide high-precision predictions. The boundary between the coarse mesh region and the fine mesh region is dynamically adjusted according to the convergence state of the optimization algorithm. When the improvement amount of the optimization objective function is less than a set value after multiple consecutive iterations, it is determined that the algorithm enters the fine mesh region. The training sample density of the fine mesh region is increased and the length scale parameter of the radial basis function is reduced.

[0067] The technical effects of the coarse and fine mesh proxy acceleration are as follows: In the global search phase, the proxy model constrained by the physical information of the coarse mesh region quickly filters out obviously unreasonable parameter combinations, avoiding time-consuming precise forward simulation of these parameter combinations and concentrating computing power on potential parameter regions; In the local fine-tuning phase, the proxy model constrained by the physical information of the fine mesh region provides high-precision prediction, enabling the optimization algorithm to accurately identify local optima and determine whether to exit; The dynamic mesh adjustment mechanism realizes adaptive allocation of computing resources, reducing the total computational load to less than one-tenth of the traditional method while ensuring global convergence.

[0068] The normalized objective function for the multi-source data fusion is expressed as follows: ;in The total objective function value, For the first Fitting error of each data source For the corresponding weights, The normalization coefficient is... For normalization items, , , The reference values ​​for the overall objective function, the errors of each data source, and the normalization term are all set to 1.

[0069] The data sources include mineral content fraction time series, solution ion concentration time series, and pH value time series. The corresponding weights are determined based on the measurement accuracy of each data source, with data sources having higher measurement accuracy having higher weights. The corresponding weights are obtained through experimental error analysis. The standard deviation is calculated by repeatedly measuring the same sample 10 times, and the reciprocal of the standard deviation is normalized and used as the corresponding weight value.

[0070] The normalization term is in the form of the dynamic parameter vector. The norm is used to penalize parameter values ​​that deviate too far from prior estimates, preventing overfitting. The normalization coefficients are determined through the cross-validation method, iterating through... to From multiple candidate values ​​within the range, select the value that minimizes the error of the validation set.

[0071] The steps of the Markov chain Monte Carlo sampling are as follows: Define the prior probability distribution of the parameters as a uniform distribution, with the upper and lower bounds of the prior probability distribution provided by the graph database; define the likelihood function as the probability density of the normally distributed error between the observed data and the simulated prediction; use the Metropolis-Hastings algorithm for sampling, randomly perturbing the current parameter value in each iteration to generate candidate parameters, calculating the ratio of the posterior probability of the candidate parameter to the posterior probability of the current parameter value; if the ratio of the posterior probabilities is greater than 1, the candidate parameter is accepted; otherwise, the candidate parameter is accepted with the probability of the posterior probability ratio. Rejected candidate parameters retain the current parameter value; set the sampling chain length to 10000, discard the first 2000 samples as the burning period, and the remaining samples constitute the posterior probability distribution of the parameters; extract the 2.5% quantile and 97.5% quantile from the posterior probability distribution of the parameters as the lower and upper bounds of the confidence interval of the parameters.

[0072] The technical benefits of the Markov chain Monte Carlo sampling are as follows: compared to point estimation methods which only provide a single optimal value for the parameters, the Bayesian inference fully characterizes the uncertainty range of the parameters through the posterior probability distribution of the parameters, enabling decision-makers to assess the reliability of prediction results based on the confidence interval of the parameters; the sampling process automatically explores multiple local optima in the parameter space, and reveals the non-uniqueness of the parameters through the multimodal characteristics of the posterior probability distribution of the parameters, providing a basis for model diagnosis and improvement.

[0073] The formula for correcting the surface area of ​​the mineral is as follows:

[0074] ;

[0075] in The specific surface area of ​​minerals at the reservoir scale is in units of , The unit for measuring specific surface area in the laboratory is: , The particle size of the laboratory core powder is measured in μm. The in-situ mineral grain size of the reservoir is expressed in μm. The surface roughness correction factor is dimensionless. and All .

[0076] The method for obtaining the in-situ mineral particle size of the reservoir is as follows: thin sections are prepared from the unground core and observed under an optical microscope. The particle size distribution of no less than 500 mineral particles is statistically analyzed, and the median particle size is taken as the representative value.

[0077] The surface roughness correction factor is obtained by measuring the root mean square value of the surface roughness of laboratory powder and reservoir core thin sections using an atomic force microscope. The ratio of the two values ​​is used as the surface roughness correction factor. Through measurements of 10 groups of cores of different reservoir types, the empirical value range of the surface roughness correction factor is statistically obtained to be 0.3 to 0.7.

[0078] The correction formula for the effective reaction volume is expressed as follows:

[0079] ;

[0080] in The effective reaction volume at the reservoir scale is in units of , The unit for effective reaction volume in the laboratory is , Since reservoir porosity is dimensionless, The porosity of the powder deposits in the laboratory core is dimensionless. The pore connectivity correction factor is dimensionless. Take the volume reference value .

[0081] The pore connectivity correction factor is obtained by measuring the pore throat radius distribution and pore connectivity of the reservoir core through mercury intrusion porosimetry, and using the ratio of the pore connectivity to the theoretical connectivity of laboratory powder as the pore connectivity correction factor.

[0082] The correction formula for the reactant transport limitation factor is expressed as follows: ;in The reactant transport restriction factor is dimensionless. The dimensionless Darmköhler number represents the ratio of reaction rate to transport rate.

[0083] The formula for calculating the Damcole number is as follows:

[0084] ;

[0085] in The unit of the reaction rate is , The characteristic length unit is , The diffusion coefficient is in units of , The unit of concentration is , Pick , Take 1 , Pick , Pick .

[0086] An excitation threshold adjustment function is defined to adjust the dynamic threshold of the pulse coding layer in the reaction network recognition model. The excitation threshold adjustment function calculates a threshold adjustment coefficient based on the concentration gradient variance within the current time window, the historical prediction error mean, and mineral phase transition detection indicators. When the threshold adjustment coefficient... At this time, a strategy of reducing the excitation threshold by 20% is adopted to enhance the ability of the reaction network recognition model to capture slow precipitation processes. When the threshold adjustment coefficient is... While maintaining the current sensitivity by keeping the excitation threshold unchanged, when the threshold adjustment coefficient The strategy of increasing the excitation threshold by 30% is adopted to suppress noise interference and highlight significant phase transition characteristics.

[0087] The calculation formula for the excitation threshold adjustment function is expressed as follows:

[0088] ;

[0089] in The unit of the current time window concentration gradient variance is... , The average prediction error over the last 5 time steps is expressed in units of . , The mineral phase transformation detection index is dimensionless and takes the value of 0 or 1. , , The weighting coefficients are set to 0.4, 0.4, and 0.2 respectively. Take the reference value for the concentration gradient variance , The historical prediction error mean reference value is taken as follows: , Take 1.

[0090] The mineral phase transition detection index is calculated as follows: if the saturation index of any mineral in the current time step crosses the threshold of 1, that is, changes from supersaturation to undersaturation, or vice versa, then the mineral phase transition detection index is set to 1; otherwise, it is set to 0.

[0091] The weight coefficients are obtained by performing a grid search on the validation set of the reaction network recognition model, traversing various combinations of the weight coefficients, and selecting the combination that gives the reaction network recognition model the highest comprehensive score in terms of sudden phase transition detection accuracy and slow precipitation prediction root mean square error. The comprehensive score is the weighted average of the two after normalization. The weight of sudden phase transition detection accuracy is 0.6, and the weight of slow precipitation prediction root mean square error is 0.4.

[0092] Alternatively, the present invention also provides a system implemented by a computer, forming a A system for constructing a kinetic model of mineralization and storage reaction is provided, wherein the computer is equipped with a readable storage medium, and the readable storage medium stores program instructions, which execute the above-described method when the program instructions are run in the computer.

[0093] The specific implementation of step S1 is as follows: First, physical property tests are performed on the target reservoir core. Porosity is measured using a helium porosimeter, and permeability is measured using a steady-state gas measurement method or a transient pressure decay method to obtain the basic physical property parameters of the reservoir. Then, the core is prepared into powder samples, and whole-rock mineral quantitative analysis is performed using X-ray diffraction. The mass percentage content of each mineral, such as quartz, plagioclase, calcite, mica, and kaolinite, is calculated using the Ritterwald refinement method. In the knowledge graph construction stage, a broad search is conducted of published... This study utilizes literature on water-rock reaction kinetics and international mineral kinetic databases (such as MINTEQ and EQ3 / 6) to extract the pre-exponential factor range, activation energy range, reaction order range, and pH index range for the dissolution and precipitation reactions of various minerals. Mineral names, reaction types, and parameter boundaries are organized into ternary structures and stored in a graph database. Relationships between minerals are also established to provide constraint boundaries for subsequent parameter optimization. The graph database supports quick querying of parameter statistical ranges for similar minerals by mineral type, providing a basis for limiting the parameter search space in hierarchical optimization strategies.

[0094] The specific implementation of step S2 is as follows: Based on the in-situ formation temperature and pressure conditions of the target reservoir, three representative temperature gradient points are designed, and the pressure is set higher than that of the target reservoir. Critical pressure (7.38) ),make sure In a supercritical state, it enhances its reactivity with formation water and minerals. Simulated formation water, prepared based on reservoir water chemical analysis data, is added to the high-pressure reactor according to the solid-liquid ratio with the core powder, and then charged... After reaching the target pressure, the system is sealed and kept at a constant temperature. Sampling times are set at logarithmic intervals of 1, 3, 7, 14, 30, and 60 days. This sampling strategy can simultaneously cover the concentration change characteristics of both the early rapid dissolution stage and the later slow precipitation stage. The active learning mechanism is activated after the first batch of fixed gradient experiments. It trains a Gaussian process model under a Bayesian optimization framework using existing experimental data. By maximizing the expected improvement criterion or the upper confidence bound criterion, it evaluates the uncertainty of kinetic parameters corresponding to different temperature-time combinations in the parameter space and recommends the next set of experimental conditions that can minimize the variance of the posterior distribution of parameters, thus achieving the goal of achieving the highest parameter identification accuracy with the fewest number of experiments.

[0095] The specific implementation of step S3 is as follows: At each sampling moment, a solid powder sample is taken out from the reaction vessel, filtered, dried, and then subjected to X-ray diffraction analysis to obtain the mass fraction of each mineral at that moment; simultaneously, liquid samples are collected, and the pH value is measured using a pH meter, and the anion concentration (including...) is measured using an ion chromatograph. , , (etc.), using inductively coupled plasma mass spectrometry to determine cation concentration (including , , , , , (etc.). All solid-liquid test data at three temperature points and six time points were compiled into a multi-temperature, multi-time-point mineral-solution time series database, with data dimensions of temperature number × time point number × species number. This time series database was input into a reaction network recognition model. The model automatically identified the active dissolution and precipitation stages of different minerals and the coupling relationships between minerals through layer-by-layer processing of pulse coding layer, temporal memory layer, dynamic feedback synaptic layer and mineral coupling attention layer, and output the core reaction network topology.

[0096] The specific implementation of step S4 is as follows: Using the solution species and phase databases built into geochemical software (such as PHREEQC), the saturation index of each mineral and the activity of each ion are calculated based on the measured solution chemical composition, constructing a complete chemical network encompassing all possible dissolution-precipitation reactions. Subsequently, the reaction contribution of each mineral over a 60-day experimental period is calculated using the following formula. For minerals whose mass change is less than 1% within 60 days or whose saturation index consistently deviates from 1, their control is shifted from kinetic to thermodynamic equilibrium or they are simply ignored, forming a simplified core reaction network. A stoichiometric correlation matrix is ​​established with minerals as rows and solution ions and gaseous components as columns. Components generated in dissolution reactions are assigned positive values, while components consumed are assigned negative values, and the opposite is true for precipitation reactions. Singular value decomposition is used to identify primary and secondary reactions, providing a basis for ranking mineral importance in subsequent stratified optimization strategies.

[0097] The specific implementation of step S5 is as follows: Based on the transition state theory, in the reaction rate data block (RATES module) of the geochemical software, dissolution and precipitation rate equations are written for the main minerals such as quartz, plagioclase, calcite, and kaolinite. The equations are in the form of... The temperature dependence of the rate constant is expressed by the Arrhenius equation. Description. To address the rigid ordinary differential equations in multi-mineral coupled reaction systems, an adaptive implicit time-step algorithm is introduced: at the beginning of each time step, the rate of change of concentration for all species is calculated, and the algorithm is applied when the ratio of the maximum to the minimum rate of change exceeds a rigid threshold. When the time interval is switched to the implicit solver, the fourth-order Gill method is used to divide the time interval into four sub-intervals. The weighted average of the slope estimates of each sub-interval (the weights are determined by the Butch table) is used as the concentration update. The implicit equation system is solved by Newton-Raphson iteration. The convergence criterion is that the relative residual is less than the tolerance value. When the stiffness ratio is lower than the threshold, the explicit solver is switched to save computational resources. The time step is dynamically adjusted according to the local truncation error.

[0098] The specific implementation of step S6 is as follows: The hierarchical optimization strategy is executed in three stages sequentially. In the first stage, a global search is performed on a single mineral at a single temperature. The parameters of other minerals are fixed to the recommended values ​​from the graph database. Only the 25℃ rate constant and activation energy of the target mineral are optimized. A genetic algorithm (population size 50, crossover probability 0.8, mutation probability 0.1, generation number 100) is used. The search range is constrained by the graph database to fluctuate by one order of magnitude above and below the literature values. Individuals exceeding the boundary are reset to boundary values. The optimization objective is the root mean square error between the experimental mass fraction and the simulated predicted mass fraction of the mineral. In the second stage, the parameters of each mineral from the first stage are used as initial values. Simultaneously, the Arrhenius parameters of all major minerals are optimized. A quasi-Newton method is used to calculate the search direction using the Hessian matrix approximation update formula. The weighted root mean square error (weights positively correlated with mineral content) is used as the objective function. The iteration termination condition is that the gradient norm is less than... Or the number of iterations exceeds 200. In the third stage, all temperature data are incorporated, and a multi-task learning framework is used to share the underlying parameters. The co-evolutionary algorithm divides the minerals into primary and secondary mineral groups according to the attention layer weights and rankings. The parameters of the two groups are alternately optimized until the global root mean square error converges. After each iteration, the element mass conservation is checked by the stoichiometric correlation matrix, and parameter combinations with a deviation of more than 0.1% are rejected.

[0099] The specific implementation of step S7 is as follows: Initial training samples with a parameter dimension 10 times larger are generated in the dynamic parameter space through Latin hypercube sampling. For each sample, an exact forward simulation is run to obtain the mineral content fraction time series as output. A Gaussian process regression surrogate model is constructed using a radial basis function kernel. The hyperparameters of the kernel function are determined by maximizing the marginal likelihood function. In subsequent optimization iterations, whenever the genetic algorithm or particle swarm optimization algorithm generates a new parameter combination, the surrogate model is first called to predict the output and estimate the prediction uncertainty. If the uncertainty is lower than the prediction uncertainty threshold, the fitness is directly evaluated using the predicted value; otherwise, the exact forward simulation is called, and the result is added to the training set to retrain the surrogate model. The prediction uncertainty threshold is determined through leave-one-out cross-validation and is dynamically updated gradually as the training set expands. The coarse-grid surrogate acceleration mechanism determines the current stage based on the continuous iterative improvement of the optimization objective function. In the coarse-grid stage, the surrogate model uses parameters with a larger length scale for rapid approximation; in the fine-grid stage, it uses parameters with a smaller length scale to provide high-precision predictions. The boundary between the two stages is dynamically adjusted according to the convergence state.

[0100] The specific implementation of step S8 is as follows: The normalization objective function for multi-source data fusion calculates the weighted sum of the fitting errors from three types of data sources: mineral content fraction time series, solution ion concentration time series, and pH value time series. The weight of each data source is obtained by normalizing the inverse of the standard deviation calculated after 10 repeated measurements. The normalization term is taken from the kinetic parameter vector. Norms, normalization coefficients are verified through cross-validation. to The range is traversed to select the value that minimizes the validation set error. Uncertainty quantification employs Bayesian inference using Markov chain Monte Carlo sampling. The prior distribution is taken as a uniform distribution within the constraints of the graph database. The likelihood function assumes that the prediction error follows a normal distribution. The Metropolis-Hastings algorithm executes a sampling chain of length 10,000, discarding the first 2,000 samples as the burn period. The remaining samples constitute the posterior probability distribution. The 2.5% and 97.5% quantiles are used as the lower and upper bounds of the parameter confidence interval, respectively, providing decision-makers with a quantitative basis for the reliability of the prediction.

[0101] The specific implementation of step S9 is as follows: The kinetic parameters obtained under laboratory core powder conditions are converted to in-situ reservoir conditions, and three corrections are performed sequentially. The mineral specific surface area correction uses the formula... The in-situ mineral grain size of the reservoir was obtained by statistically analyzing thin sections of no less than 500 mineral grains and taking the median value. The surface roughness correction factor was obtained by measuring the root mean square roughness ratio of the powder and the thin section using atomic force microscopy, with an empirical range of 0.3 to 0.7. The effective reaction volume correction was performed using the formula... The pore connectivity correction factor is determined by the ratio of the pore connectivity obtained from mercury intrusion porosimetry experiments to the theoretical pore connectivity obtained from laboratory powder experiments. The reactant transport limitation factor is determined by the Damköhler number. Substitute the calculation The process is as follows: After completing the three corrections, a reservoir-scale dynamic parameter table is output. This table is then imported into the reservoir simulator to build a grid model. Custom dynamic reaction equations for each mineral component are added, and a numerical simulation of the long-term mineralization and storage process is performed. Sensitivity analysis is conducted on key parameters such as injection rate, injection pressure, and temperature.

[0102] It should be noted that the key technologies of this invention include: firstly, a synergistic mechanism between the reaction network recognition model and the hierarchical optimization strategy. The reaction network recognition model transforms non-uniform time series data into a mineral importance ranking through a pulse coding layer and a mineral-coupled attention layer, providing an ordered solution path for the hierarchical optimization strategy and avoiding convergence failure caused by blind searching in high-dimensional parameter spaces in traditional methods; secondly, a rigid equation system solution mechanism composed of an adaptive implicit time step algorithm and the Gill method, with a rigidity ratio exceeding [missing information]. The system automatically switches to a higher-order implicit scheme, with a stability region covering the entire negative real half-axis, allowing step sizes to be increased by tens of times without causing numerical oscillations, fundamentally ensuring the convergence of multi-mineral coupled forward simulations. Thirdly, it employs a dynamic switching mechanism between the Gaussian process surrogate model and the precise forward simulation. In low-uncertainty regions, fitness is directly evaluated using the surrogate model's predicted values, significantly reducing the number of precise simulation calls and lowering the total computational cost of parameter inversion to a minimal proportion of traditional methods. The synergistic effect of these three key technologies is that the reaction network identification model provides effective initial values ​​for hierarchical optimization, the rigid solution mechanism ensures the reliable completion of each forward simulation, and the surrogate model acceleration makes high-precision Bayesian inference feasible under limited computing power. All three are indispensable, forming a complete technical chain that simultaneously ensures the convergence and identification accuracy of dynamic parameter inversion in multi-mineral coupled systems.

[0103] It should be noted that the present invention also solves the following technical problem: in multi-mineral... In water-rock reaction experiments, the choice of sampling time greatly affects the accuracy of parameter identification due to the significant differences in reaction rates among different minerals. Too sparse sampling can miss crucial phase transition information, while too dense sampling leads to a surge in experimental costs. This invention addresses the optimization problem of experimental sampling schemes through an active learning mechanism. After each batch of experiments, a Bayesian optimization framework is used to evaluate the uncertainty corresponding to each temperature-time combination in the parameter space, recommending the next set of experimental conditions that minimizes the posterior variance of the parameters, thus concentrating experimental resources on the sampling points with the highest information content. Simultaneously, the pulse coding layer in the reaction network identification model adaptively processes the non-uniform sampling data with logarithmic time intervals using dynamic thresholding. The excitation threshold adjustment function comprehensively adjusts the sensitivity based on the concentration gradient variance, historical prediction error, and mineral phase transition detection indicators, employing different excitation strategies in the rapid dissolution and slow precipitation stages to ensure that crucial phase transition information is not missed. This allows for the acquisition of highly accurate kinetic parameters even with limited experimental runs, solving the technical problem of insufficient parameter identification information caused by inappropriate sampling strategies at both the experimental design and data processing levels.

[0104] Specifically, the principle of this invention is as follows: The reason this invention can simultaneously guarantee convergence and identification accuracy is that it eliminates limiting factors from two independent dimensions. At the numerical solution level, the rigid equation system leads to divergence because the magnitude of the reaction rates of different minerals exceeds the numerical stability domain. The Gill method extends the stability domain to the entire negative real half-axis domain through a higher-order implicit scheme, and the adaptive step-size controller dynamically adjusts the step size based on local truncation errors, ensuring that each forward simulation can be completed stably, thus providing a reliable fitness assessment for optimization iterations. At the parameter search level, multi-mineral coupling results in a non-convex and high-dimensional parameter space. The reaction network identification model quantifies the coupling strength between minerals through attention weights, providing an ordered solution sequence for hierarchical optimization. The genetic algorithm performs a global coarse search within the boundaries constrained by the graph database, and the quasi-Newton method further refines it within the convergence neighborhood. These two approaches work together to avoid local extremum traps. The Gaussian process surrogate model constructs an accurate response surface approximation in the explored region, allowing the optimization algorithm to evaluate a large number of candidate parameters without frequently calling accurate simulations, significantly reducing the computational cost required to achieve the target accuracy. The technical logic of the two dimensions is independent yet synergistic, therefore the solution of this invention is logically self-consistent and sufficient.

[0105] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0106] The specific implementation of step S1 is as follows: the porosity and permeability of the target reservoir core are determined by mercury intrusion porosimetry or nitrogen adsorption method, and the mass fractions of quartz, plagioclase, calcite, mica and kaolinite are determined by whole-rock mineral quantitative analysis using X-ray diffraction. The pre-exponential factor range, activation energy range, reaction order range and pH index range of each mineral dissolution and precipitation reaction are extracted from published literature. The mineral name, reaction type and parameter boundary are used to form a ternary set and stored in the graph database to provide a priori knowledge basis for the boundary constraints of subsequent parameter search.

[0107] The specific implementation of step S2 is as follows: A multi-temperature gradient experiment is designed based on the pressure and temperature of the target reservoir, setting a first temperature point, a second temperature point, and a third temperature point, with the pressure higher than... Critical pressure (7.38) The experimental water used was simulated formation water, and the core samples were ground to powder. Logarithmic time intervals were used for sampling at 1, 3, 7, 14, 30, and 60 days. An active learning mechanism was introduced: after each set of temperature and time point experiments, a Gaussian process model between the existing experimental data and the uncertainty of the kinetic parameters was established using a Bayesian optimization algorithm. Recommended values ​​for the next set of experimental conditions were calculated by maximizing the expected improvement criterion or the upper confidence bound criterion, prioritizing the temperature and time combinations that minimized the posterior variance of the parameters to achieve the highest parameter identification accuracy with the fewest experiments.

[0108] The specific implementation of step S3 is as follows: Mineral analysis is performed on the reaction powder at each temperature point and time point. The corresponding aqueous solutions are collected to determine pH, anion content, and cation content, establishing a multi-temperature, multi-time-point mineral solution time-series database. This database is input into the reaction network recognition model for training and prediction. The data dimension of the model input layer is the number of temperatures multiplied by the number of time points multiplied by the number of species. The pulse coding layer converts the non-uniform sampling data of logarithmic time intervals into pulse sequences, and the dynamic threshold is adaptively adjusted according to the concentration standard deviation within the previous time window. The temporal memory layer uses long short-term memory network units to process the pulse sequences, including three gating mechanisms: forget gate, input gate, and output gate. The dynamic feedback synapse layer feeds back the hidden state of the temporal memory layer to the pulse coding layer to adjust the excitation threshold for the next time point. The feedback weights are obtained through training using the error backpropagation algorithm. The mineral coupling attention layer calculates the temporal correlation weights between different mineral dissolution and precipitation reactions through a multi-head attention mechanism. The rows and columns of the weight matrix correspond to the primary and secondary minerals, respectively. The training loss function is the weighted sum of the network topology prediction error and the reaction contribution prediction error, specifically described as follows:

[0109] ;

[0110] In the formula, The total loss function value is dimensionless. To measure the error in network topology prediction, cross-entropy loss is used to measure the difference between predicted and actual edges. This loss is dimensionless. To measure the error in the prediction of response contribution, the mean squared error is used to measure the deviation between the predicted weights and the actual weights; this is dimensionless. , The weighting coefficients, determined through grid search, are set to 0.6 and 0.4 respectively. The formula for calculating the excitation threshold adjustment function is as follows:

[0111] ;

[0112] In the formula, This is the threshold adjustment coefficient, dimensionless. The variance of the concentration gradient within the current time window, in units of Reference value Pick , dimensionless This represents the average prediction error over the last 5 time steps, in units of... Reference value Pick , dimensionless This is a mineral phase transition detection index, dimensionless, taking a value of 0 or 1. If the saturation index of any mineral in the current time step crosses 1 (i.e., changes from supersaturated to undersaturated or vice versa), it is set to 1; otherwise, it is set to 0. (Reference value) Take 1 , , To optimize the weighting coefficients of each term in the threshold adjustment function, a grid search was performed on the validation set. The coefficients were determined based on the highest combined score, calculated using a weight of 0.6 for the accuracy of sudden phase transition detection and 0.4 for the root mean square error of slow precipitation prediction. The coefficients were set to 0.4, 0.4, and 0.2 respectively. When the excitation threshold is reduced by 20%, Time remains unchanged, when Increase the excitation threshold by 30%.

[0113] The specific implementation method of step S4 is as follows: a chemical network containing all reactions is constructed using geochemical software, and the reaction contribution analysis formula is expressed as follows:

[0114] ;

[0115] In the formula, For the reaction Para ion Contribution of concentration change, dimensionless minerals Ions in dissolution and precipitation reactions Stoichiometric coefficients, dimensionless, reference values Take 1 For the reaction The volumetric reaction rate, in units of Reference value Pick , dimensionless For the reaction medium ion Stoichiometric coefficients, dimensionless, reference values Take 1 For the reaction The volumetric reaction rate, in units of Reference value Pick , dimensionless To achieve the summation index, all reactions are traversed. Minerals with a mass change of less than 1% over 60 days or whose saturation index consistently deviates from kinetic control are reclassified as equilibrium control or ignored, forming the core reaction network. The stoichiometric correlation matrix elements are assigned values ​​according to the following rules: components generated in dissolution reactions are assigned positive values, components consumed are assigned negative values, precipitation reactions are assigned the opposite, and elements not involved in the reaction are assigned zero values. The number of independent reactions is determined through null space analysis, and major and minor reactions are identified through singular value decomposition.

[0116] The specific implementation method of step S5 is as follows: The dissolution-precipitation rate equation is expressed as follows:

[0117] ;

[0118] In the formula, The mineral surface reaction rate, in units of Reference value Pick The rate constant is expressed in units of 1000 ppm. Reference value Pick , dimensionless Specific surface area of ​​minerals, in units of Reference value Pick , dimensionless The saturation index is calculated using the following formula: ,in It is the ion activity product. These are the solubility product constants, all dimensionless. and This is an empirical parameter for the reaction order, dimensionless, and typically takes a value of 1; the left side of the equation... Dimensionless, all three terms on the right side of the equation are dimensionless, and their dimensions correspond. The temperature dependence of the rate constant is expressed as follows:

[0119] ;

[0120] In the formula, The rate constant at 25℃, in units of Reference value Pick , dimensionless Activation energy, unit: Let be the gas constant, and take . Absolute temperature, unit: For reference temperature, take 298.15. For the temperature correction term reference value, take... Its function is to (Unit is) )and (Unit is) The product of ) multiplied by (Unit is) After this, the dimensions of the entire exponential term cancel each other out, resulting in a dimensionless equation; both sides of the equation are dimensionless, and their dimensions correspond. An adaptive implicit time step algorithm is introduced to determine the rigidity based on the order-of-magnitude range of the species concentration change rate. When the ratio of the maximum change rate to the minimum change rate exceeds the rigidity threshold... An implicit solver is selected, and the fourth-order Gear method is used to divide the time interval into four sub-intervals. The weighted average of the four slope estimates is used as the concentration update for the entire time step. The weight coefficients are determined by the Butch table. The slope values ​​of each sub-interval are obtained by iteratively solving the implicit equation system. Compared with the explicit method, the step size can be increased by tens of times without causing numerical oscillations.

[0121] The specific implementation of step S6 is as follows: In the first stage of single-temperature single-mineral optimization, a genetic algorithm is used for global search. The population size is 50, the crossover probability is 0.8, the mutation probability is 0.1, and the number of generations is 100. The parameter search range is limited to 1.5 times the upper bound and 0.5 times the lower bound of the graph database literature value. Individual parameters exceeding the boundary are reset to the boundary values. The optimization objective function is the root mean square error between the experimentally measured mass fraction and the simulated predicted mass fraction of the target mineral at a single temperature. In the second stage of single-temperature multi-mineral optimization, the optimization objective function is the weighted root mean square error between the experimentally measured mass fraction and the simulated predicted mass fraction of all major minerals at a single temperature, with weights... For the first The fitting weights of the major minerals are determined based on their content, taking the ratio of the mineral's mass fraction to the sum of the mass fractions of all major minerals. Minerals with higher content have greater weights. The search direction is calculated using a quasi-Newton method with an approximate update formula based on the Hessian matrix. The iteration terminates when the gradient norm is less than 1. Or the number of iterations exceeds 200. In the third stage of multi-temperature and multi-mineral optimization, a co-evolutionary algorithm is used to divide the minerals into two subpopulations: a primary mineral group and a secondary mineral group. The primary minerals are defined based on the weight of the mineral coupling attention layer and the top 30% of the minerals. The two subpopulations evolve alternately until the global root mean square error no longer decreases. During the optimization process, the element mass conservation is checked through the stoichiometric correlation matrix. If the total element deviation exceeds 0.1%, the parameter combination is rejected and resampling is performed.

[0122] The specific implementation of step S7 is as follows: Gaussian process regression is selected as the physical information constraint surrogate model framework. A radial basis function kernel function is used, and the length scale parameter and signal variance parameter of the kernel function are determined by maximizing the marginal likelihood function. The initial training samples are generated by Latin hypercube sampling, and the number is 10 times the parameter dimension. When the prediction uncertainty is less than the threshold, the fitness is directly evaluated by using the prediction results of the surrogate model; otherwise, exact forward simulation is called and new samples are added to the training set for retraining. The prediction uncertainty threshold is determined by the root mean square error of leave-one-out cross-validation and is updated every 10 new samples. In the coarse grid region, fewer training samples and larger length scale parameters are used to quickly filter out unreasonable parameter combinations, while in the fine grid region, more training samples and smaller length scale parameters are used to provide high-precision predictions. When the improvement of the objective function through multiple consecutive iterations is less than a set value, it is determined to enter the fine grid region and the boundary is dynamically adjusted. The total computational cost can be reduced to less than one-tenth of that of traditional methods.

[0123] The specific implementation of step S8 is as follows: The normalization objective function for multi-source data fusion is expressed as follows:

[0124] ;

[0125] In the formula, The value of the overall objective function is dimensionless. The reference value for the overall objective function is set to 1. For the first The fitting error of the data sources, which include time series of mineral content fraction, solution ion concentration, and pH value. For the first Each data source error reference value is set to 1. dimensionless For the first The fusion weights corresponding to each data source are determined by calculating the standard deviation of the same sample after repeated measurements 10 times, and then normalizing by taking the reciprocal of the standard deviation. These weights are dimensionless. The normalization coefficients are dimensionless and are obtained through cross-validation. to Determine candidate values ​​by traversing the range For the normalization term, take the dynamic parameter vector. Norms, specifically expressed as follows:

[0126] ;

[0127] In the formula, The vector of dynamic parameters to be optimized. For the prior estimated parameter vector, Representing vectors norm and The dimensions are consistent with the dynamic parameters. and Same dimensions Take 1 (and) (Reference value with the same dimensions) The entire sample is dimensionless. In Markov chain Monte Carlo sampling, the prior probability distribution is uniform, with upper and lower bounds provided by a graph database. The Metropolis-Hastington algorithm is used for sampling, with a sampling chain length of 10,000. The first 2,000 samples are discarded as the burning period. The 2.5% quantile and 97.5% quantile are extracted from the posterior probability distribution of the remaining samples as the lower and upper bounds of the parameter confidence interval.

[0128] The specific implementation method of step S9 is as follows: The mineral surface area correction formula is expressed as follows:

[0129] ;

[0130] In the formula, The specific surface area of ​​minerals at the reservoir scale, in units of Reference value Pick , dimensionless For laboratory measurement of specific surface area, the unit is... Reference value Pick , dimensionless The particle size of the laboratory core powder is given in units of 1. The in-situ mineral grain size in the reservoir is expressed in units of 1. The median particle size was obtained by preparing thin sections from unground rock cores, observing them under an optical microscope, and statistically analyzing the particle size distribution of at least 500 mineral grains. dimensionless This is a dimensionless surface roughness correction factor, obtained by measuring the root mean square value of the surface roughness of laboratory powder and reservoir core thin sections using atomic force microscopy. Its empirical range is 0.3 to 0.7; both sides of the equation are dimensionless, and their dimensions correspond. The effective reaction volume correction formula is expressed as follows:

[0131] ;

[0132] In the formula, The effective reaction volume is at the reservoir scale, in units of Reference value Pick , dimensionless The effective reaction volume in the laboratory is expressed in units of... , dimensionless Reservoir porosity, dimensionless The porosity of the powder deposits in the laboratory core is dimensionless. dimensionless The pore connectivity correction factor is dimensionless and obtained by measuring the pore throat radius distribution and pore connectivity of reservoir cores using mercury intrusion porosimetry, and then comparing it to the theoretical connectivity of laboratory powder. Both sides of the equation are dimensionless and dimensionally equivalent. The formula for the reactant transport limitation factor correction is as follows:

[0133] ;

[0134] ;

[0135] In the formula, The reactant transport limiting factor is dimensionless. The Darmquerel number is a dimensionless number representing the ratio of reaction rate to transport rate. The volumetric reaction rate is expressed in units of 1000 ppm. Reference value Pick , dimensionless The characteristic length is expressed in units of 10 ... Reference value Pick , dimensionless This is the diffusion coefficient, in units of... Reference value Pick , dimensionless The concentration of reactants in the solution, in units of Reference value Pick , dimensionless molecules With denominator All are dimensionless. The overall model is dimensionless. The reservoir-scale dynamic parameter table is imported into the simulation model to customize the long-term mineralization and storage process, and the sensitivity of injection rate, injection pressure, and temperature parameters is analyzed.

[0136] To better understand and implement this invention, the following is a specific application scenario of the invention, Example 2: To verify the effectiveness of the invention, technicians used a typical continental sandstone reservoir as a target and constructed it according to the method of the invention. A kinetic model for mineralization and storage reaction was developed, and complete parameter inversion and reservoir-scale simulation verification were carried out.

[0137] Following step S1, technicians conducted physical property tests on the obtained reservoir core, finding a porosity of 15% and a permeability of 50%. The results of whole-rock quantitative analysis by X-ray diffraction are shown in Table 1.

[0138] Table 1. Mineral composition of reservoir core

[0139]

[0140] As shown in Table 1, the reservoir is dominated by quartz, containing a certain amount of plagioclase and calcite, exhibiting typical clastic mineral assemblage characteristics, consistent with... Typical reservoir types for mineralization preservation research. Based on literature databases, the boundaries of kinetic parameters such as pre-exponential factor range and activation energy range for each mineral were organized into triplets and stored in a graph database. The activation energy range for quartz was set to 75,000 to 90,000. The activation energy range for plagioclase is set to be 50,000 to 70,000. The activation energy range for calcite is set to be 20,000 to 40,000. .

[0141] Following step S2, three temperature points (40℃, 60℃, 80℃) are designed, with the pressure uniformly set to 15. higher than Critical pressure 7.38 ,ensure It is in a supercritical state. The simulated formation water is prepared as follows: Type, mineralization 10 Core powder is ground to a particle size of less than 75 mm. The solid-liquid mixture was loaded into a high-pressure reactor at a solid-liquid ratio of 1:10. Logarithmic time interval sampling was performed at 1, 3, 7, 14, 30, and 60 days. After completing the first batch of 60-day data at three temperature points in the initial experiment, an active learning mechanism was used to establish a Gaussian process model using Bayesian optimization to assess parameter uncertainty. It was recommended to perform the model at 75℃ and 20℃. Under the given conditions, a supplementary set of experiments was conducted to further constrain the uncertainty of the plagioclase activation energy. After the completion of this set of experiments, the variance of the posterior distribution of the parameters was significantly reduced.

[0142] Following step S3, technicians performed X-ray diffraction analysis on the solid powder at 18 sampling times (3 temperature points × 6 time points), and measured the pH value of the liquid sample. , , , , , , , A mineral-solution time series database with equal component concentrations and multiple temperatures and time points was formed. This database was then input into a pre-trained reaction network recognition model. The model identified a significant acceleration in plagioclase dissolution between 7 and 14 days through the pulse coding layer. The mineral coupling attention layer calculated that the promoting weight of plagioclase on calcite precipitation was 0.67, and the inhibiting weight of kaolinite on plagioclase dissolution was 0.31.

[0143] Following step S4, the reaction contribution of minerals at each time point was calculated using PHREEQC software. The results showed that the mica mass change was only 0.3% over 60 days, and the saturation index remained below 0.05. Based on the reaction contribution analysis criteria, it was shifted from kinetic control to equilibrium control. The core reaction network includes four main pathways: plagioclase dissolution, calcite dissolution and precipitation, kaolinite precipitation, and quartz dissolution. Figure 3 As shown, plagioclase dissolution produces and Directly promotes the precipitation of kaolinite. Increased concentration promotes calcite precipitation, forming two key mineralization pathways. The stoichiometric correlation matrix has a dimension of 4×8, and singular value decomposition confirms that the number of independent reactions is 4.

[0144] Following step S5, rate equations were written for quartz, plagioclase, calcite, and kaolinite in the RATES module of PHREEQC. Initial parameter values ​​were derived from the recommended median values ​​in the graph database. The initial value for the quartz rate constant at 25°C was [value missing]. plagioclase calcite Kaolinite The adaptive implicit time-step algorithm detected that the stiffness ratio of the multi-mineral coupled system in the 7-14 day period was approximately... Exceeding the rigid threshold The system automatically switches to the fourth-order Gill method for solving the problem. The time step in this stage is extended from 0.01 days, which is limited by the explicit method, to 0.5 days, and the time taken for a single forward simulation is significantly reduced.

[0145] Following step S6, in the first stage of the hierarchical optimization strategy, the kaolinite parameters are fixed to the recommended values ​​from the graph database. A global genetic algorithm search (population size 50, generation 100) is then performed on plagioclase at a single temperature of 60℃, yielding a preliminary estimate of the plagioclase rate constant at 25℃. The activation energy is initially estimated to be 61,000. The second stage used the mineral results from the first stage as initial values, and simultaneously optimized the Arrhenius parameters of plagioclase and calcite based on the 60℃ data. A quasi-Newton method was used for fine-tuning, reducing the weighted root mean square error from 12.3% to 3.8% from the initial value, and achieving a gradient norm of [value missing]. The process terminates at a certain point. In the third stage, three temperature data sets of 40℃, 60℃, and 80℃ are simultaneously incorporated. The co-evolutionary algorithm classifies quartz and plagioclase as the primary mineral group (weighted by the attention layer and ranked in the top 30%), and kaolinite and calcite as the secondary mineral group. After five rounds of alternating optimization, the global root mean square error converges to 2.1%, and the element mass conservation deviation in each iteration step is less than 0.05%. The Gaussian process surrogate model initially trains with 80 samples under the condition of 8 parameter dimensions. During the optimization process, a total of 142 exact forward simulations are called (a reduction of about 78% compared to the traditional method without a surrogate model). The uncertainty threshold of the surrogate model prediction is dynamically updated from the initial 0.031 to 0.009 as the training set expands. The coarse and fine mesh boundaries automatically switch to fine mesh mode when the global root mean square error improvement is less than 0.1% for five consecutive iterations.

[0146] Following step S7, cross-validation results showed that the coefficient of determination for the predicted plagioclase mass fraction was 0.93, and for calcite it was 0.91, meeting the acceptance criteria. External validation experiments were conducted using an independent set of core samples from different batches (porosity 16%, similar mineral composition). The average error for the predicted mineral mass fraction was 6.2%, indicating good model generalization ability. Markov chain Monte Carlo sampling (chain length 10000, combustion period 2000) yielded the posterior probability distribution of the activation energies of each mineral. The 95% confidence interval for the plagioclase activation energy was (57200, 64800). The 95% confidence interval for the calcite activation energy is (26500, 34100). ,like Figure 2 As shown, the posterior distributions of each parameter are all unimodal, and the parameters are well identifiable.

[0147] Following step S8, during the mineral specific surface area correction, 547 mineral grains were counted in the core thin section, with the median grain size of plagioclase being 120 mm. The particle size of the laboratory powder is 38 mm. The surface roughness correction factor for atomic force microscopy was 0.52, and the effective specific surface area of ​​plagioclase after correction was [missing value]. Mercury intrusion porosimetry determined the pore connectivity to be 82%, while the theoretical connectivity of laboratory powder was 100%, with a pore connectivity correction factor of 0.82. The calculated Damcole number for plagioclase dissolution was 0.18, corresponding to a reactant transport limitation factor of 0.85, indicating that plagioclase dissolution is moderately restricted by transport under these reservoir conditions. The reservoir-scale kinetic parameter table includes corrected rate constants, activation energies, effective specific surface areas, and transport limitation factors for five minerals.

[0148] Following step S9, the reservoir-scale dynamic parameter table was imported into the CMG-GEM reservoir simulator, based on a porosity of 15% and a permeability of 50%. Build 100 ×100 ×20 A three-dimensional mesh model was created, incorporating the kinetic reactions of four minerals: quartz, plagioclase, calcite, and kaolinite, to simulate an injection that occurred over 1000 years. The mineralization and sequestration process under specific conditions. Simulation results show that mineral capture increases rapidly within 200 years after injection, then stabilizes after 500 years, with plagioclase dissolution leading to mineralization and sequestration. The mineralization was primarily captured as calcite, accounting for 34.2% of the total injected amount. This was observed at injection rates (0.5 to 3). Injection pressure (12 to 20) Sensitivity analysis was conducted on the reservoir development process and temperature (40 to 80°C). The results showed that temperature had the most significant impact on the mineralization rate, followed by injection pressure, while injection rate had the least impact. This provides a quantitative basis for optimizing reservoir development schemes.

[0149] Compared to traditional methods, this invention automatically ranks mineral importance and identifies core reaction paths using a reaction network identification model, avoiding the omission of key reactions caused by manual reliance on experience. The introduction of an adaptive implicit time step algorithm and the Gill method ensures the numerical stability of rigid equation sets at the theoretical level, rather than relying on a trial-and-error process of manually adjusting the time step. The hierarchical optimization strategy decomposes the high-dimensional parameter inversion problem into ordered low-dimensional subproblems, and with the acceleration of the surrogate model, the computational controllability of the parameter identification process is greatly improved. Bayesian inference provides a complete probabilistic description of parameter uncertainty, giving the prediction reliability assessment of reservoir simulation a theoretical basis, while traditional point estimation methods cannot quantify this uncertainty.

[0150] It should be noted that the variables involved in this invention are explained in detail in Tables 2 and 3.

[0151] Table 2. Variable Explanation Table (Part 1)

[0152]

[0153] Table 3. Variable Explanation Table (Part Two)

[0154]

[0155] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for constructing a kinetic model of carbon dioxide mineralization and storage reaction considering optimization algorithms, characterized in that, Includes the following steps: The porosity and permeability of the target reservoir core were obtained. The contents of quartz, plagioclase, calcite, mica, and kaolinite were determined by whole-rock mineral quantitative analysis using X-ray diffraction. Parameter boundaries of mineral dissolution and precipitation reactions were extracted from publicly available data and literature, and a graph database was constructed. A multi-temperature gradient experiment for carbon dioxide mineralization and storage was designed, and an active learning mechanism was introduced to recommend the next set of experimental conditions through Bayesian optimization. A multi-temperature, multi-time-point mineral-solution time series database was established, and this database was input into a reaction network identification model for training and prediction. A simplified core reaction network was formed through reaction contribution analysis, and a stoichiometric correlation matrix of mineral-solution components was established. Dissolution and precipitation rate equations were written in geochemical software, and an adaptive implicit time step algorithm combined with the Gill method was introduced to handle rigid ordinary differential equation systems. The kinetic parameters of minerals corresponding to the contents of quartz, plagioclase, calcite, mica, and kaolinite were inverted through a hierarchical optimization strategy. A surrogate model with physical information constraints was constructed to accelerate the process using coarse and fine grid proxying. Uncertainty quantification is achieved by combining a normalized objective function derived from multi-source data fusion with Bayesian inference based on Markov chain Monte Carlo sampling. Laboratory-scale parameters are converted into reservoir-scale kinetic parameters and imported into a simulation model to customize the long-term mineralization and storage process. The structure of the reaction network identification model is as follows: the input layer receives mineral content variation data and solution ion concentration time series data, which are then processed through a pulse coding layer, a temporal memory layer, a dynamic feedback synapse layer, and a mineral coupling attention layer. The output layer outputs a simplified core reaction network topology, where nodes represent the minerals and ions participating in the reaction, edges represent stoichiometric relationships, and edge weights represent the reaction contribution.

2. The method according to claim 1, characterized in that, The graph database specifically extracts the pre-exponential factor range, activation energy range, reaction order range, and pH index range of mineral dissolution and precipitation reactions from public data and literature, and constructs a triplet containing mineral name, reaction type, and parameter boundaries to store in the graph database.

3. The method according to claim 2, characterized in that, The active learning mechanism specifically involves using a Bayesian optimization algorithm to establish a Gaussian process model between existing experimental data and the uncertainty of dynamic parameters after each set of temperature and time point experiments is completed. The recommended values ​​for the next set of experimental conditions are calculated by maximizing the expected improvement criterion or the upper confidence bound criterion, and the temperature-time combination that can minimize the variance of the posterior distribution of parameters is preferentially selected.

4. The method according to claim 3, characterized in that, The carbon dioxide mineralization and storage multi-temperature gradient experiment specifically involves setting a first temperature point, a second temperature point, and a third temperature point based on the pressure and temperature of the target reservoir core, with the pressure being higher than the critical pressure of carbon dioxide. A logarithmic time interval sampling method is used to take samples at the first, second, third, fourth, fifth, and sixth moments. The experimental water is simulated formation water, and the target reservoir core is ground into powder.

5. The method according to claim 4, characterized in that, The steps for establishing the training dataset for the reaction network identification model are as follows: collect no less than 50 sets of mineral-solution time series under different temperature gradients and pressure conditions as original samples, label the real core reaction network topology, divide it into training set and test set in an 8:2 ratio, and perform data augmentation on the time series data in the training set.

6. The method according to claim 5, characterized in that, The hierarchical optimization strategy consists of three stages: First, single-temperature single-mineral optimization is performed using a genetic algorithm or particle swarm optimization to search for the approximate range of parameters globally and constrain the parameter boundaries using a graph database. Second, single-temperature multi-mineral optimization is performed using a quasi-Newton method or Bayesian optimization to perform local fine-tuning and inversion of Arrhenius parameters. Third, multi-temperature multi-mineral optimization is performed using multi-task learning or co-evolutionary algorithms to gradually release mineral parameters.

7. The method according to claim 6, characterized in that, The adaptive implicit time step algorithm specifically calculates the rate of change of all species concentrations at the beginning of each time step. When the ratio of the maximum rate of change to the minimum rate of change exceeds a rigid threshold, the implicit solver is selected; otherwise, the explicit solver is selected. The time step is dynamically adjusted based on the local truncation error estimate.

8. The method according to claim 7, characterized in that, The rigid threshold is obtained by selecting a typical water-rock reaction system for numerical experiments, recording the statistical distribution of species concentration change rates over different time periods, setting the ratio of change rates that cause the explicit solver step size to be limited to sub-second levels or divergent as the rigid threshold, and obtaining the empirical value of the rigid threshold through statistical analysis of 100 sets of core reaction simulation experiments with different mineral compositions.

9. The method according to claim 8, characterized in that, The implementation of the Gill method in the scheme is as follows: the time interval is divided into four sub-intervals using the fourth-order Gill method. The slope estimate of the species concentration is calculated in each sub-interval. The weighted average of the four slope estimates is used as the concentration update for the entire time step. The weight coefficients are determined according to the Butch table of the Gill method. The slope values ​​of each sub-interval are obtained by iteratively solving the implicit equation system.

Images