Simulation and prediction method of saturation field and pressure field in underground carbon dioxide injection process based on deep learning

Abstract

This invention relates to a deep learning-based method for simulating and predicting the saturation and pressure fields during underground carbon dioxide injection, belonging to the field of multiphase flow numerical simulation technology in underground carbon sequestration. The method includes: Step 1: Acquisition and construction of high-fidelity numerical simulation data for CO2 multiphase flow; Step 2: Preprocessing and normalization of the original spatiotemporal data; Step 3: Construction and training of an MRU-FNO model based on a fusion structure of FNO and residual U-Net; Step 4: Simulation and prediction of the saturation and pressure fields during underground carbon dioxide injection using the trained MRU-FNO model based on a fusion structure of FNO and residual U-Net. This invention maintains high prediction accuracy and exhibits excellent temporal generalization performance.

Description

Technical Field

[0001] This invention relates to a method for simulating and predicting the saturation field and pressure field during underground carbon dioxide injection based on deep learning, belonging to the field of multiphase flow numerical simulation technology in underground carbon sequestration. Background Technology

[0002] Geological carbon dioxide sequestration (GCS) captures, compresses, and injects CO2 emitted from industrial sources into deep underground geological structures for long-term storage, making it one of the most feasible carbon reduction methods currently available. In GCS, CO2 is typically injected in a supercritical state into aquifers, oil and gas reservoirs, or saline aquifers, forming a complex multiphase flow system with formation fluids. The distribution characteristics and storage efficiency of the injected CO2, influenced by various physical processes such as migration, diffusion, dissolution, and capillary action within the formation, directly affect the safety and long-term stability of the storage.

[0003] To study the transport behavior of CO2 in underground reservoirs, multiphase flow numerical simulation methods are commonly used in scientific research and engineering to model and predict underground fluid flow and mass exchange processes. However, this traditional numerical simulation method usually suffers from significant drawbacks such as high computational cost and weak time extrapolation capability.

[0004] Traditional numerical simulation methods for multiphase flow are essentially mathematical models built upon the first principles of physics. This method combines Darcy's law, the mass conservation equation, and component transport equations—partial differential equations—to characterize the CO2-water two-phase flow and mass exchange within porous media. Since these equations cannot be solved analytically under complex geological conditions, numerical techniques such as the finite volume method or finite element method are commonly used in practice. These methods discretize the continuous reservoir space and time into grids and time steps, transforming the physical equations into a large-scale algebraic equation system for iterative solution. This has led to the development of industry-standard software such as TOUGH2 and ECLIPSE. However, this traditional "physical modeling-numerical discretization" approach suffers from several fundamental technical drawbacks. First, it is computationally expensive: achieving high-fidelity prediction requires fine mesh generation and tiny time steps to analyze highly heterogeneous reservoirs and rapidly changing fluid fronts. This results in extremely high computational costs for a single simulation, severely limiting its application in critical engineering scenarios requiring massive simulations, such as uncertainty quantification, parameter inversion, and real-time optimization. Secondly, there is an inherent contradiction between physical simplification and simulation accuracy: while simplified models introduced to improve computational efficiency (such as the vertical equilibrium assumption) can achieve dimensionality reduction and acceleration, they systematically ignore or distort key physical mechanisms. For example, ignoring capillary pressure hysteresis significantly underestimates the residual CO2 capture and fails to accurately characterize complex migration behaviors such as fingering and flow around in heterogeneous reservoirs. Finally, there is a lack of long-term prediction capability and reliability: CO2 sequestration requires assessment of safety over long time scales, but traditional methods suffer from cumulative and amplified effects of numerical discretization errors, model parameter uncertainties, and insufficiently characterized physicochemical processes (such as long-term geochemical reactions) when extrapolating over long periods, casting doubt on the reliability of ultra-long-term prediction results. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention proposes a deep learning-based method for simulating and predicting the saturation field and pressure field during underground carbon dioxide injection.

[0006] This invention proposes a hybrid neural operator framework, MRU-FNO, which integrates Fourier neural operators and residual convolutional structures to efficiently predict the evolution characteristics of pressure and saturation fields during CO2 sequestration. While retaining the global modeling capabilities of Fourier neural operators, this method introduces a U-Net-like multi-scale decoding structure integrating residual convolutional blocks to enhance the model's sensitivity to boundary information and local front evolution. Numerical experiments under different geological conditions (such as permeability and porosity distributions) demonstrate that this method outperforms the traditional finite volume method (FVM) and various deep learning benchmark models proposed in recent years (such as CNN and U-Net) in terms of boundary error control, front propagation prediction accuracy, and overall generalization ability, showcasing its application potential in GCS numerical simulations under complex geological conditions.

[0007] The technical solution of this invention is as follows:

[0008] A deep learning-based method for simulating and predicting the saturation and pressure fields during underground carbon dioxide injection includes:

[0009] Step 1: Acquisition and construction of high-fidelity numerical simulation data for CO2 multiphase flow;

[0010] Step 2: Preprocessing and normalization of raw spatiotemporal data;

[0011] Step 3: Construct and train the MRU-FNO model based on the fusion structure of FNO and residual U-Net;

[0012] Step 4: Simulate and predict the saturation field and pressure field during underground carbon dioxide injection using the trained MRU-FNO model based on the fusion structure of FNO and residual U-Net.

[0013] According to a preferred embodiment of the present invention, the specific implementation process of step 1 includes:

[0014] Step 1-1: Establish a multi-scale heterogeneous reservoir physical model; including:

[0015] In a porous multiphase flow system, for any component α, the mass conservation equation is established as follows:

[0016] (1);

[0017] In equation (1), M α Components α Accumulated amount over time; t For time, s; This refers to the convective term of mass transport caused by fluid flow; This is the diffusion term representing the molecular dispersion effect; qa Indicates the mass rate of external injection or extraction;

[0018] Component cumulative term M α Further expressed as:

[0019] (2);

[0020] in, Porosity is the ratio of pore volume to total volume in a porous medium. S p For the sake of the prime minister p Saturation, dimensionless, represents phase. p The proportion of pore volume occupied; ρ p For the sake of the prime minister p Density, kg / m³ 3 ; Mass fraction, dimensionless, represents the component α In phase p The proportion of quality in;

[0021] Assuming in equation (1) diff If the diffusion term is 0, then after simplification... α The mass flux of the component is:

[0022] (3);

[0023] in, For the sake of the prime minister p The flow velocity is calculated using Darcy's law:

[0024] (4);

[0025] in, k Let mD be the absolute permeability tensor, which describes the ability of a porous medium to conduct fluids. For the sake of the prime minister p The fluid pressure gradient, Pa / m; g The acceleration vector due to gravity, in m / s² 2 ; k rp Relative permeability indicates the phase p Effective penetration rate; μ p Indicates phase p The dynamic viscosity, mPa·s;

[0026] In a multiphase flow of CO2 and brine, the capillary pressure is determined by the interfacial effect between the gaseous CO2 and liquid brine phases, and the relationship is as follows:

[0027] (5);

[0028] in, P c Capillary pressure, Pa, depends on interfacial tension and pore structure; P g , P w These represent the pressures of the CO2 phase and the brine phase, respectively, in Pa;

[0029] Steps 1-2: Design a multi-condition injection scheme and simulate to generate a dataset;

[0030] Input parameters are divided into two categories: spatially dependent field variables and global scalar parameters;

[0031] Among them, the space-dependent field variables include: horizontal permeability field k x Material anisotropy ratio k x / k y Porosity field Dimensionless; global scalar parameter scalars include: initial reservoir pressure at the top of the subsurface reservoir. P init reservoir temperature T Injection rate Q Capillary pressure scaling factor 𝜆 Irreducible water saturation S wi and the top position of the perforation. Perf top and bottom position Perf bottom ;

[0032] The output variables are the simulated gas phase saturation field and the cumulative pressure change field.

[0033] According to a preferred embodiment of the present invention, step 2 includes the following specific implementation process:

[0034] Step 2-1: Data normalization processing;

[0035] The max-min normalization method is used to map the characteristic values ​​of horizontal permeability, vertical permeability, porosity, injection location, and gas saturation to the [0,1] interval; the calculation formula is as follows:

[0036] (6);

[0037] in, For the first i The result after standardization of each parameter For the first i The values ​​of the parameters, For the first i The maximum value of each parameter. For the first i The minimum value of each parameter;

[0038] Step 2-2: Divide the training set and the test set;

[0039] The samples are divided into training and test sets; each sample includes a set of randomly sampled partial differential equation control parameters and their corresponding numerical solutions.

[0040] According to a preferred embodiment of the present invention, the MRU-FNO model includes the multi-input neural operator MIONet, i.e., a multi-branch structure, and a ResNet-enhanced U-FNO module;

[0041] The ResNet enhanced U-FNO module includes multiple coupled layers containing Fourier modules and residual convolutional blocks;

[0042] The input data is pre-encoded using MIONet, which includes: two branch networks extracting field variables and scalar variables features respectively, and a backbone network extracting time step features. The features are then fused through feature addition and multiplication operations to form an initial feature representation under time-space conditions.

[0043] The initial feature representation is input into the ResNet-enhanced U-FNO module, which retains FNO's ability to perform global modeling in the frequency domain, while enhancing the local spatial representation by introducing U-Net-style skip connections and multi-layer residual blocks; the final output is the temporal prediction result of the target physical field.

[0044] Further preferably, the specific implementation process of step 3 includes:

[0045] Step 3-1: Construct a multi-input neural operator fusion FNO module;

[0046] The multi-input neural operator MIONet learns nonlinear operator mappings between multiple input functions and a target output function through a branch-trunk architecture; including:

[0047] For a given set of input functions Branch network { B i Feature extraction is performed for each input function separately;

[0048] The backbone network consists of a multi-layer fully connected network FNN. The input is the target space-time coordinates, which are transformed by linear transformation and non-linear activation mapping in the hidden layers into a high-dimensional basis function representation. The backbone network encodes time information into the output basis function vector, thereby extracting time step features.

[0049] Finally, the basis function vectors and the input function features extracted by each branch network are fused through inner product to form a joint feature representation under time-space conditions; a basis function representation related to the target position x is generated, which is a vector generated by the backbone network for each target coordinate, including time step features and spatial location features; the final output is obtained through inner product fusion, as shown below:

[0050] (7);

[0051] Where ⟨⋅ and ⋅> represent the inner product. b 0 For bias terms; u ( x The joint feature representation under time-space conditions represents the mapping output at the target location x, i.e., the predicted output function value;

[0052] Building upon MIONet, an FNO layer is introduced. The FNO layer operates on the initial feature representation under time-space conditions, first performing a frequency-domain projection FFT, then executing pointwise multiplication convolution in the frequency domain to capture long-range dependencies, and finally mapping back to the spatial domain via inverse Fourier transform (IFFT) to obtain the enhanced feature representation. The final output is the prediction function value. A Fourier-MIONet module, i.e., a multi-input neural operator fusion FNO module, is constructed. The operations of the FNO layer include:

[0053] 1) Perform a Fast Fourier Transform (FFT) on the input features to project them into the frequency domain representation;

[0054] 2) Perform convolution operations, i.e., point-by-point multiplication, in the frequency domain to capture long-range dependencies across regions;

[0055] 3) By using Inverse Fourier Transform (IFFT), long-range dependencies across regions are mapped back to the original space;

[0056] Step 3-2: Construct the ResNet enhanced U-FNO module;

[0057] A residual learning structure, namely the ResNet module, is introduced into the U-FNO decoder;

[0058] The ResNet module is integrated into multiple scales of the U-FNO decoding path, with each scale containing multiple residual convolutional blocks. The core structure is as follows:

[0059] (8);

[0060] in, x The input feature map is obtained from the output of the previous network layer in the U-FNO decoding path; This represents a nonlinear transformation consisting of two-dimensional convolution, batch normalization, nonlinear activation, and random deactivation operations. It is an activation function. y This represents the output feature map of the current residual convolutional module;

[0061] In each U-Net downsampling layer, multiple consecutive residual blocks are added. Specifically, each residual block contains convolution operations with the same input and output channels; the input path is preserved and added to the convolution output; the output is activated by ReLU and then enters the next residual block or upsampling layer.

[0062] In implementation, the output of each layer of the decoder is:

[0063] (9);

[0064] In equation (9), This represents the hierarchical index in the network decoding path, where, Indicates the current layer. Indicates the previous layer; Indicates the first The intermediate feature map output by the layer decoding module is used as the first... Input to the layer decoding module; Indicates the first Fourier spectral convolution operator for layers; Indicates the first Pointwise convolution operator for each layer; Indicates embedded in the first The residual U-Net module in the layer decoding path includes a downsampling encoding path, a residual convolution bottleneck, and an upsampling decoding path, and achieves multi-scale feature fusion through skip connections; This represents a linear activation function, used to introduce non-linear mapping capabilities into the feature fusion results; Indicates the first Output feature map of the layer decoding module.

[0065] According to a preferred embodiment of the present invention, the MRU-FNO model employs a weighted average. l p -loss function As shown below:

[0066] ;

[0067] Where y represents the true value. This represents the predicted value from the MRU-FNO model. p It is the order of the norm. β Hyperparameters for adjusting the weights of different error terms.

[0068] A computer device includes a memory and a processor, the memory storing a computer program, and the processor executing the computer program to implement the steps of the above-described deep learning-based simulation and prediction method for saturation field and pressure field during underground carbon dioxide injection.

[0069] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the above-described deep learning-based method for simulating and predicting saturation and pressure fields during underground carbon dioxide injection.

[0070] The beneficial effects of this invention are as follows:

[0071] This invention effectively enhances the model's ability to extract and map features from complex coupled physical fields. In CO2 saturation and pressure field prediction tasks under high heterogeneity and multi-scale conditions, it can effectively reduce prediction errors and significantly improve accuracy. Compared with traditional deep learning models, this invention has a significant advantage in boundary fitting accuracy, effectively reducing the accumulation of boundary numerical errors and morphological shifts. Furthermore, this invention maintains high prediction accuracy even at time steps not covered by the training data, demonstrating excellent temporal generalization performance. Attached Figure Description

[0072] Figure 1 This is a graph showing the distribution of some parameters in the dataset;

[0073] Figure 2 Example dataset (SG represents gas saturation, dimensionless; dP represents pressure accumulation, MPa);

[0074] Figure 3 Here is a structural diagram of the MRU-FNO neural operator model;

[0075] Figure 4 Example 1 shows a schematic diagram comparing the gas saturation prediction results with the reference solution at different times.

[0076] Figure 5 Example 2 shows a schematic diagram comparing the gas saturation prediction results with the reference solution at different times.

[0077] Figure 6 Example 3 shows a schematic diagram comparing the gas saturation prediction results with the reference solution at different times.

[0078] Figure 7 Example 1 shows a schematic diagram comparing the predicted pressure accumulation results with the reference solution at different times.

[0079] Figure 8 Example 2 shows a schematic diagram comparing the predicted pressure accumulation results with the reference solution at different times.

[0080] Figure 9Example 3 shows a schematic diagram comparing the predicted pressure accumulation results with the reference solution at different times.

[0081] Figure 10 A comparison chart of saturation and pressure of various models after 30 years of injection;

[0082] Figure 11 A comparison chart showing the prediction of the gas saturation change front;

[0083] Figure 12 R is a time step interval of 1 2 Comparison chart with MAE;

[0084] Figure 13 R is a time step interval of 2 2 Comparison chart with MAE;

[0085] Figure 14 R is a time step interval of 3 2 Comparison chart with MAE. Detailed Implementation

[0086] The present invention will be further defined below with reference to the accompanying drawings and embodiments, but is not limited thereto.

[0087] Example 1

[0088] A deep learning-based method for simulating and predicting the saturation and pressure fields during underground carbon dioxide injection includes:

[0089] Step 1: Acquisition and construction of high-fidelity numerical simulation data for CO2 multiphase flow;

[0090] Step 2: Preprocessing and normalization of raw spatiotemporal data;

[0091] Step 3: Construct and train the MRU-FNO model based on the fusion structure of FNO and residual U-Net;

[0092] Step 4: Simulate and predict the saturation and pressure fields during underground carbon dioxide injection using the trained MRU-FNO model based on the fusion structure of FNO and residual U-Net. This includes:

[0093] Input Preparation and Model Loading: For the target geological carbon sequestration simulation scenario, geological parameter field information of the reservoir is acquired, including spatially dependent field variables such as permeability and porosity fields; simultaneously, injection engineering parameters are acquired, including global scalar parameters such as injection location, injection rate, and initial reservoir pressure, and the time series to be predicted is determined. The above input data are processed according to the same preprocessing and normalization rules as in the model training phase to construct a standardized input tensor. Subsequently, the pre-trained and stored MRU-FNO neural operator model parameters are loaded into the computing system to complete model initialization.

[0094] Model Inference and Output: The constructed input tensor is fed into the MRU-FNO model for forward inference computation. The model first uses a multi-input encoding network to extract and fuse features from the heterogeneous geological parameter field (field variables and scalar variables) and the time variable introduced as an explicit condition, forming an initial feature representation containing time-space constraint information. This feature representation is then used as input to the operator learning backbone module, where Fourier neural operator layers perform global spectral convolution operations in the frequency domain to capture the overall patterns and long-range spatial dependencies in the multiphase flow evolution process. Simultaneously, an embedded multi-scale residual U-Net network performs hierarchical extraction and residual correction of local spatial features. Through the synergistic effect of the global operator and local feature extraction modules, the modeling of the nonlinear evolution characteristics of multiphase flow is completed layer by layer. Finally, the decoder outputs the spatial distribution results of the CO2 gas phase saturation field or pressure field within the entire computational region at a specified time point.

[0095] Example 2

[0096] The difference between the deep learning-based simulation and prediction method for saturation and pressure fields during underground carbon dioxide injection described in Example 1 and the method described in Example 1 is as follows:

[0097] The specific implementation process of step 1 includes:

[0098] Step 1-1: Establish a multi-scale heterogeneous reservoir physical model; including:

[0099] In a porous multiphase flow system, for any component α, the mass conservation equation is established as follows:

[0100] (1);

[0101] In equation (1), M α Components α (CO2 or water) cumulative amount over time (dimensionless); t For time, s; This refers to the convective term of mass transport caused by fluid flow; This is the diffusion term representing the molecular dispersion effect; q a Indicates the mass rate of external injection or extraction;

[0102] Component cumulative term M α Further expressed as:

[0103] (2);

[0104] in, Porosity is the ratio of pore volume to total volume in a porous medium (dimensionless). ); S p For the sake of the prime minister p Saturation, dimensionless, represents phase. p The proportion of pore volume occupied (in this invention) p ∈{ CO 2 , salt}, and the sum of their saturation values ​​is 1); ρ p For the sake of the prime minister p Density, kg / m³ 3 ; Mass fraction, dimensionless, represents the component α In phase p The mass percentage of each component (the sum of the mass fractions of all components is 1);

[0105] Neglecting molecular diffusion and hydrodynamic dispersion, i.e., assuming equation (1) diff If the diffusion term is 0, then after simplification... α The mass flux of the component is:

[0106] (3);

[0107] in, For the sake of the prime minister p The flow velocity is calculated using Darcy's law:

[0108] (4);

[0109] in, k For absolute permeability tensor, describing the ability of porous media to conduct fluids (tensor in anisotropic cases), mD; For the sake of the prime minister p The fluid pressure gradient, Pa / m; g The acceleration vector due to gravity, in m / s² 2 ; k rp Relative permeability indicates the phase p The effective penetration rate (dimensionless, with a value range of [0,1]); μ p Indicates phase p The dynamic viscosity, mPa·s;

[0110] In a multiphase flow of CO2 and brine, the capillary pressure is determined by the interfacial effect between the gaseous CO2 and liquid brine phases, and the relationship is as follows:

[0111] (5);

[0112] in, P c Capillary pressure, Pa, depends on interfacial tension and pore structure; P g , P w These represent the pressures of the CO2 phase and the brine phase, respectively, in Pa; generally speaking, P g Greater than P w .

[0113] The multi-scale heterogeneous reservoir physical model is a comprehensive physical model system composed of multiphase mass conservation equations, interphase coupling relationships, and flow control equations. Specifically, the spatiotemporal evolution of CO2 and water in porous media is described by establishing component mass conservation equations. Furthermore, parameters such as porosity, phase saturation, and phase density in the component accumulation terms are used to characterize the influence of reservoir property differences at different spatial locations on flow behavior. Simultaneously, the absolute permeability tensor and relative permeability function with spatially heterogeneous distribution are introduced into Darcy's law to achieve a unified characterization of fluid migration characteristics under reservoir structural conditions at different scales. In addition, by setting the capillary pressure relationship between CO2 and water, the influence of interface effects on fluid distribution and pressure evolution during multiphase flow is further reflected.

[0114] Equation (1) constitutes the basic governing equation of the multi-scale heterogeneous reservoir physical model, which is used to describe the mass conservation relationship between CO2 and the components of brine in the heterogeneous porous medium.

[0115] In equation (2), porosity φ Spatial location related parameters φ(x) It is used to characterize the structural differences of reservoirs at different spatial scales.

[0116] In equations (3) and (4), by introducing the spatially heterogeneous permeability tensor into Darcy's law, the model can simultaneously describe the migration behavior of fluids in complex reservoir structures at different scales.

[0117] In equation (5), the capillary pressure relationship is used to characterize the interfacial interaction between CO2 and water at the pore scale. It is an important component of the multiphase flow coupling mechanism and is sensitive to pore structures of different scales.

[0118] Multi-scale characteristics are reflected in: (1) on a spatial scale, reservoir properties (porosity, permeability, etc.) vary significantly at different spatial locations; (2) on a temporal scale, rapid short-term evolution and slow long-term migration coexist after CO2 injection. Heterogeneity is manifested through the spatial distribution of porosity. φ(x) Penetration rate k(x) These parameters are described by other parameters, which are discretized at a grid scale to form a non-uniform parameter field.

[0119] Steps 1-2: Design a multi-condition injection scheme and simulate to generate a dataset;

[0120] The simulation was conducted using the commercial reservoir simulator ECLIPSE (e300) developed by Schlumberger. The simulation system settings were as follows: reservoir thickness ranged from 12.5 to 200 m, extending radially to 100,000 m, with an approximately infinite boundary condition assumption; closed boundary conditions were used at the top and bottom of the reservoir; and the simulation period was set to 30 years. Each sample in the dataset contains a set of randomly sampled partial differential equation control parameters and their corresponding numerical solutions. The input parameters were divided into two categories: spatially dependent field variables and global scalar parameters.

[0121] Figure 1 This is a graph showing the distribution of some parameters in the dataset; Figure 1 In the figure, (a), (b), (c), (d), (e), and (f) show the frequency distribution of key reservoir parameters, namely horizontal permeability field, porosity, injection rate, initial pressure, reservoir temperature, and irreducible water saturation, respectively.

[0122] Among them, the space-dependent field variables include: horizontal permeability field k r (10) -3 ≤ k r ≤ 10 5 mD Material anisotropy ratio k x / k z (1 ≤ k x / k z ≤ 150 (dimensionless); porosity field ϕ Dimensionless; global scalar parameter scalars include the initial reservoir pressure at the top of the subsurface reservoir. P init (10 ≤ P init ≤ 30 MPa), reservoir temperature T (35 ≤ T ≤ 170℃), injection rate Q (0.2 ≤ Q ≤ 2, MT / year), capillary pressure scaling factor 𝜆 (0.3 ≤ 𝜆 ≤ 0.7 (dimensionless), irreducible water saturation S wi (0.1 ≤ Swi ≤ 0.3 (dimensionless) and the top position of the perforation. Perf top (0 ≤ Perf top ≤ 200,m) and bottom position Perf bottom (0 ≤ Perf bottom ≤ 200,m);

[0123] The output variables are the simulated gas-phase saturation field and the cumulative pressure change field. To satisfy the Gaussian distribution assumption, all permeability fields underwent logarithmic transformation before use. Furthermore, the entire dataset underwent rigorous data quality control, including grid convergence verification and physical plausibility screening, to ensure the stability and reliability of the simulation results. The specific distribution of each parameter is as follows: Figure 1 As shown.

[0124] It is important to note that, in each case, the reservoir thickness... b It is a randomly sampled variable that determines the size of the reservoir input, ranging from 12.5 to 200 meters. When b For depths less than 200m, zero-fill is used for the area outside the actual reservoir to ensure all inputs have the same shape. Then, in each case, from 0 to... b Randomly select from within the range Perf top and Perf bottom All input and output fields in the dataset are stored on a spatial grid of (96, 200). The output is saved as 24 time snapshots, from 1 day to 30 years: {1 day, 2 days, 4 days, ..., 14.8 years, 21.1 years, 30 years}. The training set contains 4500 cases, and the test set contains 500 cases.

[0125] The specific implementation process of step 2 includes:

[0126] Step 2-1: Data normalization processing;

[0127] All input parameters are divided into two categories: spatially dependent field variables, such as the horizontal permeability field ( k x ), porosity field ( φ ) etc.; global scalar variables, such as initial reservoir pressure ( P init ), reservoir temperature ( T ), injection rate ( QAll permeability fields undergo logarithmic transformation before use to satisfy the Gaussian distribution assumption. This invention employs the Min-Max Normalization method to map the characteristic values ​​of horizontal permeability, vertical permeability, porosity, injection location, and gas saturation to the [0,1] interval. Characteristic values ​​refer to the numerical values ​​of various physical parameters constituting the model input tensor in the sample dimension, specifically including the values ​​of spatially dependent field variables in each discrete grid cell, and the values ​​of global scalar variables under corresponding sample conditions. For spatially dependent field variables, the characteristic values ​​correspond to the physical property parameters such as the horizontal permeability field and porosity field at each spatial location; for global scalar variables, the characteristic values ​​correspond to the scalar input values ​​such as initial reservoir pressure, reservoir temperature, and injection rate under a single simulation condition. These characteristic values ​​collectively constitute the input feature set of the model, and in the normalization process, they are linearly mapped based on their global minimum and maximum values ​​in the training dataset, thereby achieving a unified representation of parameters with different dimensions and numerical scales. The calculation formula is:

[0128] (6);

[0129] in, For the first i The result after standardization of each parameter For the first i The values ​​of the parameters, For the first i The maximum value of each parameter. For the first i The minimum value of each parameter; the normalization parameters are saved during the training phase for subsequent inverse normalization processing and result visualization.

[0130] Step 2-2: Divide the training set and the test set;

[0131] The samples were divided into training and test sets in a 9:1 ratio. Each sample included a set of randomly sampled partial differential equation control parameters and their corresponding numerical solutions. The partial differential equation control parameters refer to the set of input conditions used to determine the solution of the multiphase flow control equations, including spatially dependent field variables (such as permeability field and porosity field) and injection condition scalar variables (such as injection rate and injection location). The numerical solution refers to the spatial distribution of the pressure field or CO2 saturation field at a given time step obtained by solving the multi-scale heterogeneous reservoir physical model established in step 1-1 using a high-fidelity numerical simulation method under the above control parameter conditions.

[0132] Before being input into the neural operator model, the aforementioned control parameters are preprocessed and normalized according to the method in step 2-1. The normalized values ​​constitute the input feature values ​​of the model. Therefore, each sample can be formally represented as a set of normalized input feature values. The mapping relationship between it and its corresponding high-fidelity numerical solution is used for the neural operator model to learn the partial differential equation solution operator.

[0133] It covers a variety of reservoir thickness, permeability, and injection conditions to ensure the diversity and representativeness of data distribution. The final dataset is as follows: Figure 2 As shown.

[0134] The MRU-FNO model proposed in this invention consists of two parts, and the detailed architecture of the model is as follows: Figure 3 As shown;

[0135] The MRU-FNO model includes the multi-input neural operator MIONet, i.e., a multi-branch structure, and a ResNet-enhanced U-FNO module;

[0136] The ResNet enhanced U-FNO module includes multiple coupled layers containing Fourier modules and residual convolutional blocks;

[0137] The input data is pre-encoded using MIONet, which includes: two branch networks extracting features of field variables (such as permeability, porosity, etc.) and scalar variables (injection rate, reservoir temperature, etc.) respectively, and a backbone network extracting time step features. The features are then fused through feature addition and multiplication operations to form an initial feature representation under time-space conditions.

[0138] The initial feature representations are input into the ResNet-enhanced U-FNO module. On the one hand, it retains FNO's ability in global modeling in the frequency domain; on the other hand, it enhances the local spatial representation by introducing U-Net-style skip connections and multi-layer residual blocks, especially in terms of response capability in the boundary region. The final output is the temporal prediction result of the target physical field.

[0139] The overall forward propagation of the model can be expressed as follows: the input multi-physical quantity control parameters are fed into a multi-branch network to extract features, and after fusion, they are fed into a hierarchical network containing a Fourier module and ResConv. Spatiotemporal coupling modeling is performed in a multi-scale space, and finally the temporal prediction results of the target physical field are output.

[0140] The specific implementation process of step 3 includes: This invention proposes a fusion model structure (MRU-FNO) based on multi-input neural operators and residual-enhanced U-FNO to solve the problem that traditional Fourier neural operators are insufficient in modeling local high-frequency changes and boundary effects in complex multiphase flow systems such as geological CO2 sequestration.

[0141] Step 3-1: Construct a multi-input neural operator fusion FNO module;

[0142] The multi-input neural operator MIONet learns nonlinear operator mappings between multiple input functions and a target output function through a branch-trunk architecture; including:

[0143] For a given set of input functions Branch network { B i Feature extraction is performed for each input function separately; the branch network is divided into two categories:

[0144] branch1 extracts the local physical property features of each spatial unit for spatial field variables through linear mapping, boundary filling, and tensor rearrangement.

[0145] branch2, for scalar variables, maps the original scalar to a high-dimensional feature vector through a fully connected layer, capturing its global impact on the system output.

[0146] The backbone network consists of a multi-layer fully connected network FNN. The input is the target space-time coordinates, which are transformed by linear transformation and non-linear activation mapping in the hidden layers into a high-dimensional basis function representation. The backbone network encodes time information into the output basis function vector, thereby extracting time step features.

[0147] Finally, the basis function vectors and the input function features extracted by each branch network are fused through inner product to form a joint feature representation under time-space conditions; a basis function representation related to the target position x is generated, which is a vector generated by the backbone network for each target coordinate, including time step features and spatial location features; the final output is obtained through inner product fusion, as shown below:

[0148] (7);

[0149] Where ⟨⋅ and ⋅> represent the inner product. b 0 This serves as a bias term; the advantage of this framework lies in its ability to simultaneously handle multiple input types, such as spatially dependent field variables (e.g., permeability field, porosity field) and global control parameters (e.g., injection rate, initial pressure) in geological systems, providing flexibility for modeling flow processes under complex geological conditions. In the MIONet framework, the branch network extracts input function features, and the backbone network extracts space-time conditional basis functions. The two are fused through inner product to generate an initial feature representation under time-space conditions. Ultimately, u ( x The joint feature representation under time-space conditions represents the mapping output at the target location x, i.e., the predicted output function value;

[0150] Building upon MIONet, an FNO layer is introduced to enhance global modeling capabilities. The FNO layer operates on the initial feature representation under time-space conditions, first performing a frequency-domain projection FFT, then executing pointwise multiplication convolutions in the frequency domain to capture long-range dependencies, and finally mapping back to the spatial domain via inverse Fourier transform (IFFT) to obtain the enhanced feature representation, ultimately outputting the prediction function value. A Fourier-MIONet module, i.e., a multi-input neural operator fusion FNO module, is constructed, significantly enhancing the model's global modeling capabilities. The operations of the FNO layer include:

[0151] 1) Perform Fast Fourier Transform (FFT) on the input features to project them into the frequency domain representation; the input features are the spatial field variables (permeability field, porosity field), scalar variables (injection rate, perforation location, etc.) and spatiotemporal coordinates processed by the branches and backbone network mentioned above;

[0152] 2) Perform convolution operations, i.e., point-by-point multiplication, in the frequency domain to capture long-range dependencies across regions;

[0153] 3) By using Inverse Fourier Transform (IFFT), long-range dependencies across regions are mapped back to the original space;

[0154] In this invention, the Fourier-MIONet module, through joint modeling of spatially dependent field variables and global control parameters, can achieve efficient prediction of complex flow processes while maintaining physical consistency, laying the foundation for further construction of a more refined hybrid neural operator model (MRU-FNO).

[0155] Step 3-2: Construct the ResNet enhanced U-FNO module;

[0156] The ResNet-enhanced U-FNO module is a hybrid neural operator decoding structure centered on Fourier neural operators (FNOs), designed to simultaneously characterize the global operator properties and local multi-scale details of partial differential equation solutions. This module uses spectral convolution as its backbone to model long-range spatial dependencies and low-frequency dominant behaviors among physical field variables, and employs a residual learning mechanism to ensure the training stability and mapping consistency of deep networks. Building upon this, a multi-scale U-Net is introduced as a local correction branch to compensate for the limitations of traditional FNOs in characterizing high-frequency structures and non-stationary regions.

[0157] Structurally, the ResNet-enhanced U-FNO module includes multi-layer Fourier spectral convolutions, pointwise residual convolution paths, and an embedded residual U-Net. Each layer's feature update is obtained by adding three parts: the spectral convolution branch handles global information propagation, the pointwise convolution branch provides local linear residual mapping, and the U-Net branch introduced in the latter half of the network performs multi-scale spatial reconstruction of the features. These three paths are fused at the same spatial resolution and output through a non-linear activation function, thus forming a unified residual enhancement operator representation.

[0158] The U-Net module employed is an encoder-decoder structure with a residual bottleneck, achieving multi-scale feature interaction through downsampling, upsampling, and skip connections. Its core function is not to replace the operator learning capability of FNO, but rather to act as a local correction operator, finely correcting high-gradient, high-frequency physical features such as those near injection wells, phase interfaces, and strongly heterogeneous regions. The multi-layer residual convolution introduced at the bottleneck further enhances the nonlinear representation capability and effectively alleviates the gradient degradation problem in deep networks.

[0159] By organically combining the global Fourier neural operator with the local residual U-Net, this ResNet-enhanced U-FNO module significantly improves the modeling accuracy for complex, multi-scale physical processes while maintaining the operator's generalization ability and computational efficiency. This structure is particularly suitable for applications involving multiphase flow problems such as geological carbon sequestration, where long-range pressure propagation and strong local nonlinear responses coexist. It provides an effective network architecture design approach for efficient and high-precision PDE proxy modeling.

[0160] While traditional FNO and its extended models (such as U-FNO) possess powerful global modeling capabilities in the frequency domain, their ability to capture local high-frequency variations and boundary effects is limited. To enhance the model's ability to effectively extract local spatial details and mitigate prediction errors in boundary regions, this invention introduces a residual learning structure, namely the ResNet module, into the U-FNO decoder. This improvement enables the network to better maintain and propagate spatial gradient information, thereby enhancing the accuracy and continuity of local predictions. It is particularly suitable for complex environments with abrupt changes in the frontier and boundaries during saturation and pressure accumulation scenarios in geological carbon sequestration processes.

[0161] The ResNet module is integrated into multiple scales of the U-FNO decoding path, with each scale containing multiple residual convolutional blocks. The core structure is as follows:

[0162] (8);

[0163] in, x The input feature map is obtained from the output of the previous network layer in the U-FNO decoding path; This represents a nonlinear transformation consisting of two-dimensional convolution, batch normalization, nonlinear activation, and random deactivation operations. It is an activation function. y This represents the output feature map of the current residual convolutional module;

[0164] Input feature map in the residual module here x Instead of the original input data, it is obtained from the output of the previous network layer in the U-FNO decoding path. Specifically, the input data first undergoes frequency domain mapping and low-frequency pattern learning through multiple layers of Fourier spectral convolution operators, while linear mapping compensation is performed on each spatial location through pointwise convolutional residual paths. Subsequently, the outputs of each branch are superimposed and nonlinear activation is applied to obtain an intermediate feature map containing global operator information and local feature information. x The intermediate feature maps maintain the same spatial resolution and number of channels as the current decoding scale and serve as input to subsequent residual U-Net modules.

[0165] Nonlinear transformation function of residual convolution module It consists of the following operations in sequence:

[0166] (1) Two-dimensional convolution operation, used to extract local spatial features from the input feature map;

[0167] (2) Batch normalization operation is used to normalize the convolution output to improve the stability of network training;

[0168] (3) Nonlinear activation operation, using the LeakyReLU activation function to enhance feature representation and avoid gradient vanishing;

[0169] (4) Dropout is used to suppress overfitting.

[0170] y This represents the output feature map of the current residual convolutional module. Output Feature Map y In the input feature map x After nonlinear transformation Then, compared with the input feature map x The result is obtained by performing element-wise residual stacking and applying the activation function σ(·) to the stacked result. This output feature map maintains the same spatial resolution and number of channels as the input feature map and is used as the input to the next network layer in subsequent residual convolution modules, upsampling modules, or decoding paths.

[0171] In each U-Net downsampling layer, multiple consecutive residual blocks are added. Specifically, each residual block contains convolution operations with the same input and output channels; the input path is preserved and added to the convolution output; the output is activated by ReLU and then enters the next residual block or upsampling layer.

[0172] By fusing spectral domain modeling, linear channel mixing, and local residual convolution, the decoder can take into account both global structure and local details when reconstructing complex saturation and pressure fields.

[0173] In implementation, the output of each layer of the decoder is:

[0174] (9);

[0175] In equation (9), This represents the hierarchical index in the network decoding path, where, Indicates the current layer. Indicates the previous layer; Indicates the first The intermediate feature map output by the layer decoding module has its features fused with operator information from the preceding network layers, serving as the first... Input to the layer decoding module; Indicates the first The Fourier spectral convolution operator of the layer is used to perform Fourier transform on the input feature map, perform learnable linear mapping on a preset number of low-frequency patterns in the frequency domain, and return to the physical space through inverse Fourier transform. Indicates the first The pointwise convolution operator of the layer performs a linear mapping on the channel features at each spatial location while maintaining the spatial resolution. This is equivalent to a one-dimensional convolution operation and is used to provide local residual compensation and enhance the stability of the network. Indicates embedded in the first The residual U-Net module in the layer decoding path includes a downsampling encoding path, a residual convolution bottleneck, and an upsampling decoding path. It also achieves multi-scale feature fusion through skip connections to correct local high-frequency features in the input feature map. This represents a linear activation function, used to introduce non-linear mapping capabilities into the feature fusion results, thereby improving the training stability of deep networks; Indicates the first The output feature map of the layer decoding module. This output feature map maintains the same spatial resolution and number of channels as the input feature map and serves as the input to subsequent decoding layers or output mapping layers for further feature processing.

[0176] This structure explicitly integrates frequency domain features, channel-level information fusion, and multi-scale spatial information reconstruction, giving the model stronger representation capabilities, especially in handling multi-scale front advancements and non-stationary boundary changes.

[0177] The final output of the MRU-FNO model consists of two physical fields: a gas saturation field. SG ( r, z, t ) and pressure field dP ( r, z, t ).in, r The horizontal radial coordinate represents the distance from the injection center; z The vertical coordinates represent the position in the reservoir thickness direction; t The time coordinate represents different points in time during the simulation.

[0178] The MRU-FNO model uses weighted averages. l p -loss function As shown below:

[0179] ;

[0180] Where y represents the true value. This represents the predicted value from the MRU-FNO model. p It is the order of the norm. β To adjust the hyperparameters for different error term weights, the initial learning rate was set to 0.001 during training, and gradually decreased according to a set decay rate (0.05%) to avoid the model falling into oscillations or overfitting in the later stages. Simultaneously, the training process combined a masking mechanism with a learning rate scheduler, effectively enhancing the model's robustness under multi-scale and complex geological conditions.

[0181] The detailed structural parameters of the MRU-FNO model proposed in this invention are shown in Table 1.

[0182] Table 1. Detailed structural parameters of MRU-FNO;

[0183]

[0184] In terms of data structure design, considering the significant differences in reservoir thickness among different samples, and to ensure consistency in input and output dimensions and stability of the training process, this invention employs zero-padding for non-reservoir regions during the data preprocessing stage. Simultaneously, to prevent invalid regions from interfering with model optimization, a separate spatial mask matrix is ​​generated for each sample. During training, the loss function is calculated only within the mask region, effectively constraining the model to focus solely on the physical properties of the effective reservoir region, significantly improving training efficiency and generalization performance.

[0185] In terms of input feature design, in addition to including basic geological physical parameters (such as permeability and porosity), the MRU-FNO model also explicitly introduces spatial grid coordinate information as an additional input channel to enhance the model's ability to perceive the spatial correlation of heterogeneous geological bodies and improve its analytical level of complex spatial distribution characteristics.

[0186] To verify the performance advantages of the proposed MRU-FNO model in predicting multiphase flow in CO2 geological storage, this step designed a systematic verification scheme, including performance comparison experiments with the Fourier-MIONet model, time extrapolation capability tests, and structural ablation experiments. All experiments were based on the high-fidelity numerical simulation data generated in step 1, and the training and test sample distributions were kept consistent to ensure the objectivity and consistency of the comparison results.

[0187] Gas saturation field prediction verification:

[0188] To verify the effectiveness of the proposed MRU-FNO model in the gas saturation field prediction task, the Fourier-MIONet model, a representative high-performance model, was selected as the baseline model for comparison, and its overall prediction accuracy was evaluated. The experiment used three metrics—mean absolute error (MAE), mean squared error (MSE), and coefficient of determination (R²)—to systematically compare the performance of different models on the same test set.

[0189] Table 2 shows the evaluation metrics of the two models in the saturation prediction task. It can be observed that the MRU-FNO model proposed in this invention exhibits superior performance across all metrics. Compared to the Fourier-MIONet model, the mean absolute error of MRU-FNO decreases by 67.5%, demonstrating its advantage in capturing the characteristics of complex spatial fields under multivariate coupling.

[0190] Table 2 MAE, MSE, and R of MRU-FNO and Fourier-MIONet 2 Comparison table;

[0191]

[0192] To further demonstrate the model's fitting effect on the gas saturation field at different time scales, Figure 4 , Figure 5 , Figure 6 Visual comparison charts are provided showing the prediction results of the MRU-FNO and Fourier-MIONet models with the reference solution at three key moments (323 days, 7.3 years, and 30 years). Figure 4 Example 1 shows a schematic diagram comparing the gas saturation prediction results with the reference solution at different times. Figure 5 Example 2 shows a schematic diagram comparing the gas saturation prediction results with the reference solution at different times. Figure 6 Example 3 shows a schematic diagram comparing the gas saturation prediction results with the reference solution at different times. Figures 4 to 6 In the graph, the horizontal axis represents the radial coordinate r, with units of meters; the vertical axis represents the vertical coordinate z, with units of meters; the color scale on the right side of the graph represents the gas saturation, which is dimensionless.

[0193] Each column represents a fixed time point. The first row shows the actual simulation results, followed by the prediction outputs of the two models and their corresponding error assessment metrics (MAE and MSE).

[0194] It can be clearly observed from the figure that: (1) In the early stage (323 days), both MRU-FNO and Fourier-MIONet models can capture the saturation front well, but MRU-FNO has a smaller error and a smoother boundary transition; (2) In the middle stage (7.3 years), Fourier-MIONet shows obvious front lag, while MRU-FNO can still fit the overall saturation band shape more accurately; (3) In the long-term evolution stage (30 years), MRU-FNO still maintains good restoration of the saturation front and low value region, while the comparative model Fourier-MIONet shows edge distortion and numerical shift.

[0195] Validation of pressure field accumulation field prediction:

[0196] Besides gas saturation, the pressure field is an important physical quantity for evaluating the safety and sequestration efficiency of CO2 injection processes. Its accumulation mode and diffusion trend directly affect fault stability and spillover risk. Therefore, high-precision pressure prediction capability is also a key performance indicator for multiphase flow simulation models.

[0197] Figure 7 , Figure 8 and Figure 9 This presentation compares the performance of the MRU-FNO and Fourier-MIONet models in the pressure field prediction task. Figure 7 Example 1 shows a schematic diagram comparing the predicted pressure accumulation results with the reference solution at different times. Figure 8 Example 2 shows a schematic diagram comparing the predicted pressure accumulation results with the reference solution at different times. Figure 9 Example 3 shows a schematic diagram comparing the predicted pressure accumulation results with the reference solution at different times. Figures 7 to 9 In the graph, the horizontal axis represents the radial coordinate r, with units of meters (m); the vertical axis represents the vertical coordinate z, with units of meters (m); the color scale on the right side of the graph represents the pressure P, with units of MPa. From... Figures 7 to 9 As can be observed, MRU-FNO can accurately capture pressure accumulation centers and their diffusion boundaries, and maintains good numerical stability during long-term evolution. Table 3 shows the MAE, MSE, and R² indices for each model. MRU-FNO leads in all indices, with an R² of 99.1%, higher than Fourier-MIONet (98.6%).

[0198] Table 3 MAE, MSE, and R of MRU-FNO and Fourier-MIONet 2 Comparison table;

[0199]

[0200] It is worth noting that in regions with significant changes in boundary pressure, MRU-FNO, relying on its deep residual mechanism, can more effectively characterize the nonlinear response, effectively mitigating the numerical drift problem that occurs in traditional neural operators in regions with drastic gradient changes. In summary, the MRU-FNO model also demonstrates excellent accuracy and stability in modeling pressure accumulation processes.

[0201] Comparison results of different models:

[0202] To further evaluate the effectiveness and generalization ability of the proposed MRU-FNO, this invention introduces several representative baseline models for supplementary comparison, including traditional convolutional neural networks (CNN), Fourier neural operators (FNO), and their multi-scale improved version U-FNO as controls.

[0203] Table 4 presents the quantitative error comparison results of each model in the gas saturation field and pressure field prediction tasks.

[0204] Table 4. Quantitative comparison results of MAE and MSE of different models;

[0205]

[0206] The results show that CNN exhibits relatively large errors in the prediction of both types of physical quantities, especially in the pressure field prediction task, where its MAE and MSE are significantly higher than those of operator-based learning methods. This is mainly because CNN relies on local convolutional kernels for feature extraction, making it difficult to effectively capture the long-range spatial dependencies and global physical constraints present in the CO2-water two-phase seepage process.

[0207] In contrast, FNO, by directly learning the global operator structure of the input-output mapping in the Fourier domain, significantly improves prediction accuracy compared to CNN, validating the applicability of operator learning methods in complex multiphase seepage problems. However, standard FNO still has certain limitations in characterizing high-frequency local structures, resulting in relatively high errors in regions with large saturation boundaries and pressure gradient changes. U-FNO, by introducing a U-Net encoder-decoder structure to extract multi-scale information features, further improves the prediction accuracy of saturation and pressure fields, indicating that multi-scale information plays an important role in describing the heterogeneous structure and multi-scale seepage behavior within reservoirs.

[0208] like Figure 10As shown, MRU-FNO achieved the best results across all evaluation metrics. In the gas saturation prediction task, its MAE and MSE decreased to 0.0027 and 0.000305, respectively; in the pressure field prediction task, its MAE and MSE also significantly outperformed the baseline model. These results demonstrate that MRU-FNO can more accurately characterize the nonlinear coupling relationships between multiple physics fields during CO2 injection.

[0209] Frontier prediction and boundary accuracy testing:

[0210] The gas saturation front (i.e., the evolutionary boundary where saturation rapidly transitions from low to high) in multiphase flow simulations is a concentrated manifestation of the interaction between phase interface transport, strongly nonlinear transport, and complex geological conditions. This region not only involves rapid changes in relative permeability and capillary pressure but is also often strongly modulated by abrupt changes in stratigraphic heterogeneity, leading to significant prediction biases in traditional numerical models due to discretization errors and linearization assumptions. Therefore, accurately capturing the spatiotemporal evolution trajectory of the front, especially its location, morphology, and dynamic stability, becomes a key indicator for evaluating the physical generalization ability and stability of simulators.

[0211] To evaluate the model's performance in the frontier region, such as Figure 11 As shown, a color map of the saturation field evolution after 30 years of continuous injection was plotted, comparing the prediction performance of MRU-FNO and Fourier-MIONet. The results show that MRU-FNO can stably and continuously capture the leading edge propagation process, with its predicted leading edge position highly consistent with the actual solution and a smooth and consistent interface morphology. In contrast, Fourier-MIONet exhibits local ambiguity and positional shifts at the leading edge, reflecting the limitations of the baseline model in handling strongly nonlinear, multi-scale coupled boundaries.

[0212] Further observation reveals that the prediction errors of MRU-FNO are mainly concentrated in regions far from the leading edge and where saturation changes gradually, while the error at the leading edge itself is relatively small. This indicates that the model can concentrate its limited representational capabilities on the most complex and critical interface regions of the physical process, effectively preventing the spread of errors in sensitive areas. These visualization results intuitively demonstrate the high accuracy and robustness of MRU-FNO in leading edge tracking.

[0213] The superior performance of MRU-FNO can be attributed to the synergistic effect of its multi-scale residual structure and physics-guided design. At the architectural level, the residual-enhanced U-Net encoder-decoder path achieves cross-level aggregation of information from the pore scale to the field scale, enabling the model to keenly identify abrupt local saturation changes induced by abrupt formation changes. Meanwhile, the core FNO module ensures efficient modeling of long-range dependencies through global convolution in Fourier space, thus maintaining the overall coherence of the saturation front over long-term evolution. At the physical level, the essence of the saturation front is the dynamic equilibrium boundary of momentum and mass exchange between phases. MRU-FNO, through its deep nonlinear mapping, implicitly and accurately learns the abrupt changes in the relative permeability curve and the hysteresis relationship between capillary pressure and saturation within this region, among other complex coupling processes. This fusion of data-driven and physical characteristics results in predictions that are not only numerically accurate but also possess a high degree of physical consistency.

[0214] In summary, MRU-FNO demonstrated outstanding modeling capabilities in the critical region of evolutionary boundaries, validating its practical potential in scenarios such as CO2 geological sequestration and multiphase flow modeling.

[0215] Model performance validation at unknown time steps:

[0216] In long-term carbon sequestration simulations, the predicted time often does not perfectly match the time range of the training data. Therefore, the model needs to have the ability to extrapolate to unseen time steps. To verify the generalization performance of the proposed MRU-FNO model in the time dimension, this study designed three sets of training experiments with different time step intervals: using intervals of {1,2,3}, we selected {12,8,6} time steps as the training set.

[0217] Figure 12 R is a time step interval of 1 2 Comparison chart with MAE; Figure 13 R is a time step interval of 2 2 Comparison chart with MAE; Figure 14 R is a time step interval of 3 2 Comparison chart with MAE. (See attached image.) Figure 12-14As shown, even with large time step intervals in the training data, the MRU-FNO model maintains high prediction accuracy. This excellent generalization ability is attributed to the MRU-FNO architecture's use of temporal information as a core input feature, ensuring the continuity and physical consistency of saturation / pressure changes over time. Compared to traditional FNO models that use independent channels to predict different time steps, MRU-FNO naturally captures continuous temporal evolution patterns, thus providing reliable predictions even with sparse training time or testing time exceeding the training range. This characteristic enables the model to effectively capture the time-varying patterns of CO2 sequestration, providing reliable predictive capabilities for long-term carbon sequestration simulations.

[0218] Example 3

[0219] A computer device includes a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement the steps of the deep learning-based simulation and prediction method for saturation field and pressure field during underground carbon dioxide injection as described in Embodiment 1 or 2.

[0220] Example 4

[0221] A computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of the deep learning-based simulation and prediction method for saturation field and pressure field during underground carbon dioxide injection as described in Embodiment 1 or 2.

Claims

1. A method for simulating and predicting the saturation field and pressure field during underground carbon dioxide injection based on deep learning, characterized in that, include: Step 1: Acquisition and construction of high-fidelity numerical simulation data for CO2 multiphase flow; Step 2: Preprocessing and normalization of raw spatiotemporal data; Step 3: Construct and train the MRU-FNO model based on the fusion structure of FNO and residual U-Net; Step 4: Simulate and predict the saturation field and pressure field during underground carbon dioxide injection using the trained MRU-FNO model based on the fusion structure of FNO and residual U-Net; The MRU-FNO model includes the multi-input neural operator MIONet, i.e., a multi-branch structure, and a ResNet-enhanced U-FNO module; The ResNet enhanced U-FNO module includes multiple coupled layers containing Fourier modules and residual convolutional blocks; The input data is pre-encoded using MIONet, which includes: two branch networks extracting field variables and scalar variables features respectively, and a backbone network extracting time step features. The features are then fused through feature addition and multiplication operations to form an initial feature representation under time-space conditions. The initial feature representation is input into the ResNet-enhanced U-FNO module, which retains FNO's ability in frequency domain global modeling and enhances local spatial representation by introducing U-Net-style skip connections and multi-layer residual blocks; the final output is the temporal prediction result of the target physics field. The specific implementation process of step 3 includes: Step 3-1: Construct a multi-input neural operator fusion FNO module; The multi-input neural operator MIONet learns nonlinear operator mappings between multiple input functions and a target output function through a branch-trunk architecture; including: For a given set of input functions Branch network { B i Feature extraction is performed for each input function separately; The backbone network consists of a multi-layer fully connected network FNN. The input is the target space-time coordinates, which are transformed by linear transformation and non-linear activation mapping in the hidden layers into a high-dimensional basis function representation. The backbone network encodes time information into the output basis function vector, thereby extracting time step features. Finally, the basis function vectors and the input function features extracted by each branch network are fused through inner product to form a joint feature representation under time-space conditions; a basis function representation related to the target position x is generated, which is a vector generated by the backbone network for each target coordinate, including time step features and spatial location features; the final output is obtained through inner product fusion, as shown below: （7）； Where ⟨⋅ and ⋅> represent the inner product. b 0 For bias terms; u ( x The joint feature representation under time-space conditions represents the mapping output at the target location x, i.e., the predicted output function value; Building upon MIONet, an FNO layer is introduced. The FNO layer operates on the initial feature representation under time-space conditions, first performing a frequency-domain projection FFT, then executing pointwise multiplication convolution in the frequency domain to capture long-range dependencies, and finally mapping back to the spatial domain via inverse Fourier transform (IFFT) to obtain the enhanced feature representation. The final output is the prediction function value. A Fourier-MIONet module, i.e., a multi-input neural operator fusion FNO module, is constructed. The operations of the FNO layer include: 1) Perform a Fast Fourier Transform (FFT) on the input features to project them into the frequency domain representation; 2) Perform convolution operations, i.e., point-by-point multiplication, in the frequency domain to capture long-range dependencies across regions; 3) By using Inverse Fourier Transform (IFFT), long-range dependencies across regions are mapped back to the original space; Step 3-2: Construct the ResNet enhanced U-FNO module; A residual learning structure, namely the ResNet module, is introduced into the U-FNO decoder; The ResNet module is integrated into multiple scales of the U-FNO decoding path, with each scale containing multiple residual convolutional blocks. The core structure is as follows: （8）； in, x The input feature map is obtained from the output of the previous network layer in the U-FNO decoding path; This represents a nonlinear transformation consisting of two-dimensional convolution, batch normalization, nonlinear activation, and random deactivation operations. It is an activation function. y This represents the output feature map of the current residual convolutional module; In each U-Net downsampling layer, multiple consecutive residual blocks are added. Specifically, each residual block contains convolution operations with the same input and output channels; the input path is preserved and added to the convolution output; the output is activated by ReLU and then enters the next residual block or upsampling layer. In implementation, the output of each layer of the decoder is: （9）； In equation (9), This represents the hierarchical index in the network decoding path, where, Indicates the current layer. Indicates the previous layer; Indicates the first The intermediate feature map output by the layer decoding module is used as the first... Input to the layer decoding module; Indicates the first Fourier spectral convolution operator for layers; Indicates the first Layer-wise pointwise convolution operator; Indicates embedded in the first The residual U-Net module in the layer decoding path includes a downsampling encoding path, a residual convolution bottleneck, and an upsampling decoding path, and achieves multi-scale feature fusion through skip connections; This represents a linear activation function, used to introduce non-linear mapping capabilities into the feature fusion results; Indicates the first Output feature map of the layer decoding module.

2. The method for simulating and predicting the saturation field and pressure field during underground carbon dioxide injection based on deep learning as described in claim 1, characterized in that, The specific implementation process of step 1 includes: Step 1-1: Establish a multi-scale heterogeneous reservoir physical model; including: In a porous multiphase flow system, for any component α, the mass conservation equation is established as follows: (1)； In equation (1), M α Components α Accumulated amount over time; t For time, s; This refers to the convective term of mass transport caused by fluid flow; This is the diffusion term representing the molecular dispersion effect; q a Indicates the mass rate of external injection or extraction; Component cumulative term M α Further expressed as: (2)； in, Porosity is the ratio of pore volume to total volume in a porous medium. S p For the sake of the prime minister p Saturation, dimensionless, represents phase. p The proportion of pore volume occupied; ρ p For the sake of the prime minister p Density, kg / m³ 3 ; Mass fraction, dimensionless, represents the component α In phase p The proportion of quality in; Assuming in equation (1) diff If the diffusion term is 0, then after simplification... α The mass flux of the component is: (3)； in, For the sake of the prime minister p The flow velocity is calculated using Darcy's law: (4)； in, k Let mD be the absolute permeability tensor, which describes the ability of a porous medium to conduct fluids. For the sake of the prime minister p The fluid pressure gradient, Pa / m; g The acceleration vector due to gravity, in m / s² 2 ; k rp Relative permeability indicates the phase p Effective penetration rate; μ p Indicates phase p The dynamic viscosity, mPa·s; In a multiphase flow of CO2 and brine, the capillary pressure is determined by the interfacial effect between the gaseous CO2 and liquid brine phases, and the relationship is as follows: （5）； in, P c Capillary pressure, Pa, depends on interfacial tension and pore structure; P g , P w These represent the pressures of the CO2 phase and the brine phase, respectively, in Pa; Steps 1-2: Design a multi-condition injection scheme and simulate to generate a dataset; Input parameters are divided into two categories: spatially dependent field variables and global scalar parameters; Among them, the space-dependent field variables include: horizontal permeability field k x Material anisotropy ratio k x / k y Porosity field Dimensionless; global scalar parameter scalars include: initial reservoir pressure at the top of the subsurface reservoir. P init reservoir temperature T Injection rate Q Capillary pressure scaling factor 𝜆 Irreducible water saturation S wi and the top position of the perforation. Perf top and bottom position Perf bottom ; The output variables are the simulated gas phase saturation field and the cumulative pressure change field.

3. The method for simulating and predicting the saturation field and pressure field during underground carbon dioxide injection based on deep learning, as described in claim 1, is characterized in that... The specific implementation process of step 2 includes: Step 2-1: Data normalization processing; The max-min normalization method is used to map the characteristic values ​​of horizontal permeability, vertical permeability, porosity, injection location, and gas saturation to the [0,1] interval; the calculation formula is as follows: （6）； in, For the first i The result after standardization of each parameter For the first i The values ​​of the parameters, For the first i The maximum value of each parameter. For the first i The minimum value of each parameter; Step 2-2: Divide the training set and the test set; The samples are divided into training and test sets; each sample includes a set of randomly sampled partial differential equation control parameters and their corresponding numerical solutions.

4. The method for simulating and predicting the saturation field and pressure field during underground carbon dioxide injection based on deep learning according to any one of claims 1-3, characterized in that, The MRU-FNO model uses weighted averages. l p -loss function As shown below: ； Where y represents the true value. This represents the predicted value from the MRU-FNO model. p It is the order of the norm. β Hyperparameters for adjusting the weights of different error terms.

Images