A method for identifying displacement fronts in aquifer gas storage construction

Abstract

The present invention belongs to the technical field of gas storage construction, and specifically relates to a method for identifying the displacement front of an aquifer gas storage construction. The present invention constructs an aquifer gas storage model grid system, locates the displacement front of the aquifer gas storage model based on interface energy gradient analysis, calculates the interface energy of each grid and its three-dimensional gradient, and calibrates the dynamic front boundary with a maximum gradient threshold of 10%. The "coarse-fine" grid fractal dimension method is further used to quantify the complexity of the front, and the front is divided into a stable area, a transition area, and a risk area by quartiles. The present invention directly characterizes the dynamic evolution characteristics of the gas-water interface through the interface energy gradient mechanism, combines the fractal dimension to quantitatively identify high-risk gas channeling areas, and provides technical support for the construction optimization and gas channeling prevention and control of heterogeneous aquifer gas storage.

Description

Technical Field

[0001] The present invention belongs to the field of gas storage, and in particular relates to a method for identifying the displacement front of an aquifer gas storage reservoir. Background Art

[0002] Aquifer gas storage is governed by a combination of sedimentary environment, tectonic evolution, and hydrogeological conditions. Reservoirs have complex pore structures (dominated by a dual medium of pores and fractures), require high caprock sealing, and exhibit significant heterogeneity (permeability coefficient of variation > 0.8). Due to the long-term natural hydrodynamic balance and perturbations of injection-production cycles, the dynamic evolution of the gas-water interface is complex, making precise control of the gas-water front a key technical challenge for the efficient construction and safe operation of gas storage.

[0003] Artificial gas storage spaces are created through gas injection and drainage, and its core mechanisms include immiscible displacement, gravity-induced displacement, and coordinated pressure regulation. Currently, domestic and international scholars have conducted in-depth research on gas-water front morphology, critical gas crossflow pressure, and seal evaluation based on core displacement experiments and three-dimensional geomechanical simulations. However, traditional methods rely on interwell pressure monitoring or saturation interpolation, making it difficult to capture dynamic heterogeneous phenomena such as fingering and bypass of the gas-water front in highly heterogeneous reservoirs in real time. This results in a gas storage reservoir's effective working gas utilization rate of less than 65%, increasing the risk of seal failure.

[0004] Chinese patent number CN107515246B, "A Method for Detecting the Displacement Front of Carbon Dioxide Flooding," proposes a non-intrusive method for detecting the displacement front of carbon dioxide flooding. This method uses a curve of average sound velocity to identify the abrupt change, thereby determining the time when the gas displacement front reaches the position measured by the ultrasonic probe, and thus the location of the displacement front. However, this method can ambiguously identify the abrupt change point in highly heterogeneous reservoirs (fractures cause acoustic scattering), resulting in large positioning errors.

[0005] Chinese patent number CN106761613B, "A Method for Determining the CO2 Displacement Front by Well Testing," discloses a method for determining the CO2 displacement front by using relevant parameters and solving a relationship equation for the CO2 displacement front. This method suffers from nonlinear pressure gradients in highly heterogeneous reservoirs (fracture conductivity distorts the pressure field), model failure, and long well testing cycles (requiring several weeks of pressure buildup testing), making it impossible to track the front in real time.

[0006] Chinese patent number CN112182897B, "Method and Apparatus for Determining the Displacement Front in Interwell Water Channels in Double-High Reservoirs," provides a method for determining the displacement front by determining the length of the oil-water two-phase region and the location of the equivalent piston displacement point based on the displacement front shape equation. This method assumes that the influence of capillary forces on the stability of the displacement front is ignored.

[0007] Chinese patent number CN118187783A, "A Method for Predicting the CO2 Miscible Displacement Front in Low-Permeability Reservoirs Based on Bottomhole Flowing Pressure," discloses a method for predicting the CO2 miscible displacement front based on bottomhole flowing pressure, using basic reservoir parameters to generate different miscible front prediction models. However, this method cannot quantify the complexity of the front.

[0008] The above patents all propose a method for identifying the displacement front, but all of them are based on CO2-crude oil miscible / immiscible displacement, which is essentially different from gas-water two-phase immiscible displacement. The existing technology does not involve the identification of the gas-water displacement front during the construction of gas storage facilities, and none of them proposes a method for using the fractal dimension method to identify the interfacial energy area to determine the displacement front. Summary of the Invention

[0009] The present invention aims to provide a method for identifying the displacement front during the construction of an aquifer gas storage reservoir. The method identifies the position of the displacement front during the gas flooding process by using the interfacial energy generated at the contact surface between the injected gas and formation water. This method aims to solve the problem of unclear displacement front position during the construction of an aquifer gas storage reservoir, and provides guidance for improving the effect of gas flooding and controlling and delaying the formation of gas channeling.

[0010] Interfacial energy refers to the additional energy generated by interfacial tension at the interface between two immiscible phases. Its magnitude is related to the gas-liquid properties, interfacial tension, and reservoir microstructural parameters. During the process of injecting gas to displace formation water, interfacial energy is generated at the point where the gas contacts the water in the formation, as the two phases are immiscible. Studying this interfacial energy can help determine the location of the displacement front.

[0011] The present invention adopts the interfacial energy gradient mechanism to integrate capillary force and three-dimensional saturation gradient to locate the displacement front of the aquifer gas storage model, adopts the "coarse-fine" grid fractal dimension method to quantify the complexity of the displacement front, and divides the stable zone, transition zone, and risk zone into quartiles. The interfacial energy calculation formula can calculate the dynamic changes of grid interface energy. No additional monitoring equipment is required, and low-cost front prediction can be achieved based on geological modeling and numerical simulation.

[0012] The specific steps include:

[0013] (1) Based on the geological data of the gas storage aquifer obtained from geophysical logging, empirical formulas for capillary force, permeability, and porosity were obtained through indoor experiments;

[0014] Specifically, the relevant geological data of the target layer for database construction is obtained based on geophysical logging, including porosity , permeability , water saturation , reservoir pressure , formation water viscosity , formation water density Indoor experiments were conducted to obtain gas-water interfacial tension, gas viscosity, irreducible water saturation, and porosity distribution index. Mercury injection tests were then conducted on n cores of the target reservoir. By fitting the relationship between permeability, porosity, and capillary force for the n cores, the empirical formula for capillary force, permeability, and porosity was obtained as follows:

[0015] ;

[0016] ;

[0017] Where, P c is the capillary force, Pa; is the air-water interfacial tension, ; is the porosity of the water layer; is the water layer permeability, ; is the dimensionless function of water saturation; is water saturation; is the irreducible water saturation; is the porosity distribution index.

[0018] (2) Construct an aquifer gas storage model and discretize it into orthogonal grids, with each grid storing the permeability k i , porosity φ i , water saturation S wi and capillary forces P ci The aquifer gas storage model includes an injection well model, a drainage well model, a layer velocity model, a layer model, a porosity model and a permeability model; the construction process is as follows:

[0019] (2-1) Import the wellhead coordinate data, well trajectory and logging data of gas injection well a and drainage well b into Petrel software to build models of gas injection well and drainage well;

[0020] (2-2) Import 3D seismic volume and well-seismic calibration data, combine them with seismic data, establish interval velocity model, and constrain wellbore time-depth conversion;

[0021] (2-3) Use seismic interpretation to build a model and perform fault interpretation using ant volume tracking; use seismic volume tracking controlled by seed points to interpret horizons, generate a layer model, and convert the interpreted horizons into the depth domain;

[0022] (2-4) Seismic attribute modeling establishes porosity model and permeability model; the reservoir is discretized into orthogonal grids, with grid numbers i=1, 2, ..., N, and each grid stores permeability , porosity , water saturation and capillary forces .

[0023] (3) Based on the aquifer gas storage model, the working data of the well is imported, the simulation time parameters are set, and the water saturation value of each grid at different locations and times of the target gas storage model is obtained based on the differential solution of the control equation. , capillary force , permeability ;

[0024] Specifically, the injection rate of gas injection well a is X million cubic meters per day, and the drainage rate of drainage well b is Y million cubic meters per day. The calculation time module sets a time step of m days, and the simulation cycle is n days, for a total of time steps. In the numerical calculation of the entire data set, the gas-displacement water process satisfies the following governing equation:

[0025] Equations of motion:

[0026] ;

[0027] ;

[0028] Where, v gi is the gas phase velocity, m / s ; k gi is the gas phase permeability, mD ; μ gi is the gas phase viscosity, mPa·S ; v wi is the liquid phase velocity, m / s ; k wi is the liquid permeability, mD ; μ wi is the liquid viscosity, mPa·S ; and are the pressure gradients of the gas and liquid phases, respectively, .

[0029] The mass conservation equation:

[0030] ;

[0031] ;

[0032] ;

[0033] Where, is the porosity; 、 are the gas and liquid phase densities, ; 、 are gas saturation and water saturation, respectively; is the gas phase source term, ; is the liquid phase sink term, .

[0034] Equation of state:

[0035] ;

[0036] Where, is the gas phase pressure, Pa; is the molar mass of the gas, ; Z is the gas compressibility factor; R is the gas constant, ; T is temperature, K.

[0037] The above equations are differentially solved using the simulation module of Petrel software to obtain the formation pressure, gas saturation, and water saturation values for each grid at different locations and times of the target gas storage.

[0038] (4) Call the aquifer gas storage model in step (3) to calculate the gas-water interfacial tension in each grid at time t σ i , capillary force P ci , water saturation S wi , permeability k i , calculate the interface energy of each grid E i , the calculation formula is as follows:

[0039] ;

[0040] ;

[0041] ;

[0042] ;

[0043] Where, E iis the interfacial energy of the ith grid, J; is the air-water interfacial tension, mN / m; P ci is the capillary pressure in the i-th grid, Pa; is the dimensionless function of water saturation of the i-th grid; S wi is the water saturation of the i-th grid, 0≤ S wi ≤1; is the irreducible water saturation; is the porosity distribution index, which is taken as 2.1 here. is the water saturation gradient amplitude, m -1 ; x, y, z correspond to the three orthogonal directions of the grid; is the grid volume, m³; X, Y, Z are the length, width, and height of the grid, respectively, in meters, all of which are 5 meters here; is the correction factor of the i-th grid; For the i The permeability of each grid, mD; is the average permeability, mD; 、 are the viscosities of water and gas, mPa·S, respectively.

[0044] (5) After obtaining the interface energy of each grid Then, the interface energy gradient amplitude of each grid is calculated :

[0045] ;

[0046] ;

[0047] ;

[0048] ;

[0049] Where (x, y, z) is the grid index; x 、 y 、 z Three orthogonal directions of the grid respectively; is the grid spacing, m; 、 are the interface energies at the grid (x+1, y, z) and the grid (x-1, y, z), J; is the interface energy gradient amplitude of the ith grid, J / m.

[0050] (6) The maximum value of the calculated interface energy gradient amplitude is recorded as ,Will 10% as the threshold , the interface energy gradient amplitude is greater than or equal to the threshold The grid is marked as the displacement front grid.

[0051] (7) Calculate the fractal dimension of the displacement front, evaluate the front complexity based on the fractal dimension, and control the gas injection parameters in different zones;

[0052] (7-1) The displacement front grid obtained by numerical simulation is discretized into a point matrix and covered with two-scale boxes of "coarse-fine"; the number of boxes required to cover the front using the coarse grid is , using a fine grid to count the number of boxes needed to cover the leading edge is After obtaining the number of boxes required for the two-scale grid coverage, the fractal dimension method is used to reflect the complexity growth rate under the scale change of the displacement front, which is calculated as follows:

[0053] ;

[0054] Where, For the i The fractal dimension of a region, i =1, 2, 3…100; For the i The number of boxes under the coarse grid statistics of each region; For the i The number of boxes under the fine grid statistics of each area; R is the side length of the coarse square grid, m; r is the side length of the fine square grid, m;

[0055] (7-2) Count the values of all leading edge points D f , calculate the quartiles; Q 1 is the 25th percentile, which is the boundary of the more regular area of the front edge; Q 3 is the 75th percentile, which is the boundary of the more fragmented area at the front edge;

[0056] Stable zone: D f ≤Q 1;

[0057] Transition Zone: Q 1< D f ＜Q 3;

[0058] Risk areas: D f ≥Q 3;

[0059] Comparing the calculated quantiles with the boundaries of each zone reveals the regularity of the gas-water front at that moment. For the front located in transition and risk zones, technical policy adjustments are implemented for the corresponding gas injection wells to achieve a regular displacement front, improve the sweep coefficient, and enhance the efficiency of gas storage construction.

[0060] Compared with the prior art, the present invention has the following advantages:

[0061] (1) Accurate dynamic identification of strong heterogeneous fronts

[0062] This method breaks through the limitations of traditional interwell pressure monitoring and saturation interpolation. Through interfacial energy gradient analysis (integrating a capillary force prediction model with 3D saturation gradient calculation), it directly quantifies the dynamics of the gas-water interface using physical and chemical mechanisms, reducing positioning errors. This method is particularly suitable for highly heterogeneous reservoirs with high permeability variation coefficients.

[0063] (2) Quantitative risk classification using fractal dimension

[0064] The fractal dimension of the "coarse-fine" grid is used to evaluate the complexity of the front edge, and the stable zone, transition zone, and risk zone are divided into quartiles to quantitatively identify high-risk gas channeling areas. Combined with dynamic gas injection adjustment, the sweep coefficient is improved and the gas channeling risk is reduced.

[0065] (3) Strong adaptability and low cost advantages

[0066] By using the correction factor α (coupling of permeability heterogeneity and gas-water viscosity ratio) and three-dimensional gradient calculation, it is adapted to highly heterogeneous reservoirs, reduces monitoring costs compared to traditional well testing methods, supports real-time decision-making, and is suitable for the construction of gas storage reservoirs in my country's terrestrial highly heterogeneous aquifers. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 This is a schematic diagram of the fractal calculation of the gas drive front morphology, where black represents the gas phase, white represents the water phase, and the grid lines represent the grid divisions;

[0068] Figure 2 Schematic diagram of the interface energy model, where red represents Ei≠0 and blue represents Ei=0;

[0069] Figure 3 Schematic diagram of the gas-driven water front. DETAILED DESCRIPTION

[0070] The exploration data of Qianjiang aquifer gas storage revealed that the aquifer is mainly composed of sandstone, buried at a depth of 1000m, and the water layer pressure is P The average permeability is 10MPa. k 300mD, average porosity is 0.2, the formation water viscosity The average water saturation is 0.9 is 0.8.

[0071] (1) The air-water interfacial tension can be obtained by conducting indoor experiments is 50 mN / m, gas viscosity is 0.025, the bound water saturation The porosity and permeability of each core and the overall porosity distribution index of the 100 cores were obtained by conducting indoor experiments on 100 cores of the target layer. ,in The relationship between permeability, porosity and capillary force of 100 cores is fitted, and the empirical formula of capillary force, permeability and porosity is obtained as follows:

[0072] ;

[0073] ;

[0074] Where, P c is the capillary force, Pa; is the air-water interfacial tension, Here we take 50 mN / m; is the porosity of the water layer; is the water layer permeability, ; is the dimensionless function of water saturation; is the water saturation, ; is the irreducible water saturation, which is taken as 0.3 here; is the porosity distribution index, which is taken as 2.1 here.

[0075] (2) Construct an aquifer gas storage model and discretize it into orthogonal grids, with each grid storing the permeability , porosity , water saturation and capillary forces The aquifer gas storage model includes an injection well model, a drainage well model, a layer velocity model, a layer model, a porosity model and a permeability model; the construction process is as follows:

[0076] (2-1) The wellhead coordinate data of one gas injection well and four drainage wells were imported into the Petrel software. The wellhead coordinates of the gas injection well are (1124, 3607, 98), and the wellhead coordinates of the four drainage wells are (997, 3811, 28), (1288, 3765, 32), (1010, 3402, 17), and (1315, 3465, 40), respectively. The well trajectory and logging data were used to establish the models of the gas injection well and drainage well.

[0077] (2-2) Import 3D seismic volume and well-seismic calibration data, combine them with seismic data, establish an interval velocity model, and constrain the wellbore time-depth conversion.

[0078] (2-3) Use seismic interpretation to build a model and use ant volume tracking to interpret faults; use seismic volume tracking controlled by seed points to interpret horizons, generate a layer model, and convert the interpreted horizons into the depth domain.

[0079] (2-4) Seismic attribute modeling: Porosity model, permeability model and water saturation model were established with a model size of 5000m × 5000m × 1500m.

[0080] The aquifer gas storage model is discretized into orthogonal grids with grid numbers i=1, 2, ..., N. Each grid stores the permeability , porosity , water saturation and capillary forces .

[0081] (3) Based on the aquifer gas storage model, the working data of the well is imported, the simulation time parameters are set, and the water saturation value of each grid at different locations and times of the target gas storage model is obtained based on the differential solution of the control equation. , capillary force , permeability ;

[0082] One gas injection well has an injection rate of 600,000 cubic meters per day, and four drainage wells have a drainage rate of 200,000 cubic meters per day. The computational time module is set to a 10-day time step, with a full simulation period of 12 months, for a total of 36 time steps. In the numerical calculation of the entire data set, the gas-displacement water process satisfies the following governing equations:

[0083] Equations of motion:

[0084] ;

[0085] ;

[0086] Where, v gi is the gas phase velocity, m / s ; k gi is the gas phase permeability, mD ; μ gi is the gas phase viscosity, mPa·S , here we take 0.025 mPa·S ; v wi is the liquid phase velocity, here it is taken as 0.9mPa·S ; k wi is the liquid permeability, mD ; μ wi is the liquid viscosity, mPa·S ; and are the pressure gradients of the gas and liquid phases, respectively, .

[0087] The mass conservation equation:

[0088] ;

[0089] ;

[0090] ;

[0091] Where, is the porosity; 、 are the gas and liquid phase densities, , here we take 、 ; 、 are gas saturation and water saturation, respectively; is the gas phase source term, ; is the liquid phase sink term, .

[0092] Equation of state:

[0093] ;

[0094] Where, is the gas phase pressure, Pa; is the molar mass of the gas, ; Z is the gas compressibility factor; R is the gas constant, , here we take 8.314 ; T is temperature, K.

[0095] The above equations are solved by differential method using the simulation module of Petrel software to obtain the water saturation value of each grid at different locations and times of the target gas storage. , capillary force , permeability .

[0096] (4) Call the aquifer gas storage model in step (3) to calculate the gas-water interfacial tension in each grid at time t , capillary force , water saturation , permeability , calculate the interface energy of each grid ;

[0097] ;

[0098] Where, is the interface energy of the grid (15,25,8), J; is the air-water interfacial tension, mN / m, here it is taken as 50 mN / m; For the i Capillary pressure within each grid, Pa; For the i A dimensionless function of water saturation of each grid; is the water saturation of the grid (15,25,8), which is 0.45 here; is the bound water saturation, which is taken as 0.3 here; is the porosity distribution index, which is taken as 2.1 here. is the water saturation gradient amplitude, m -1 ; x 、 y 、 z They correspond to the three orthogonal directions of the grid; is the grid volume, m³; X, Y, Z are the length, width, and height of the grid, respectively, in meters, all of which are 5 meters here; is the correction factor for the grid (15,25,8); is the permeability of the grid (15, 25, 8), mD, which is 600mD here; is the average permeability, mD, here it is taken as 300mD; 、 are the viscosities of water and air, which are taken as 0.9 and 0.025 respectively.

[0099] (5) After obtaining the interface energy of each grid Then, the interface energy gradient amplitude of each grid is calculated :

[0100] ;

[0101] Where (x, y, z) is the grid index; x, y, and z are the three orthogonal directions of the grid respectively; is the grid spacing, m, here it is 5m; 、 are the interface energies at the grids (16, 25, 8) and (14, 25, 8), J, which are taken as 200 J and 4000 J respectively (the same applies to the y and z directions); is the interface energy gradient amplitude of the grid (15, 25, 8), J / m.

[0102] (6) The maximum value of the calculated interface energy gradient amplitude , here 10% of 3100 J / m is used as the threshold , the interface energy gradient amplitude is greater than or equal to the threshold The grid is marked as the displacement front grid. Greater than threshold , so the grid (15, 25, 8) is marked as the displacement front grid.

[0103] (7) Calculate the fractal dimension of the displacement front, evaluate the front complexity based on the fractal dimension, and control the gas injection parameters in different zones.

[0104] (7-1) The displacement front grid obtained by numerical simulation is discretized into a point matrix and divided into 100 regions, which are covered by two scales of "coarse-fine" boxes. The number of boxes required to cover the front using the coarse grid is , using a fine grid to count the number of boxes needed to cover the leading edge is After obtaining the number of boxes required for the two-scale grid coverage, the fractal dimension method is used to reflect the complexity growth rate under the scale change of the displacement front, which is calculated as follows:

[0105] ;

[0106] Where, For the i The fractal dimension of a region, i =1, 2, 3…100; for the i The number of boxes under the coarse grid statistics of the region, Here we take 1850; For the i The number of boxes under the fine grid statistics of each area is 8900 here; R is the side length of the coarse square grid, m, here it is 5m; r is the side length of the fine square grid, m, here it is taken as 1m.

[0107] (7-2) Count the values of all leading edge points D f , calculate the quartiles; Q 1 is the 25th percentile, which is 1.12 here, and is the boundary of the more regular area of the front edge; Q 3 is the 75th percentile, which is 1.58 here, marking the boundary of the more fragmented area at the front edge;

[0108] Stable zone: D f ≤Q 1;

[0109] Transition Zone: Q 1< D f ＜Q 3;

[0110] Risk areas: D f ≥Q 3;

[0111] Comparing the calculated quantiles with the boundaries of each zone reveals the regularity of the gas-water front at that moment. For the front located in transition and risk zones, technical policy adjustments are implemented for the corresponding gas injection wells to achieve a regular displacement front, improve the sweep coefficient, and enhance the efficiency of gas storage construction.

Claims

1. A method for identifying the displacement front of an aquifer gas storage reservoir, characterized in that: Interfacial energy gradient analysis is used to locate the displacement front of the aquifer gas storage model, and the "coarse-fine" grid fractal dimension method is used to quantify the complexity of the displacement front. The specific steps are as follows: (1) Based on the geological data of the gas storage aquifer obtained from geophysical logging, empirical formulas for capillary force, permeability, and porosity were obtained through indoor experiments; (2) Construct an aquifer gas storage model and discretize it into orthogonal grids, with each grid storing the permeability k i , porosity φ i , water saturation S wi and capillary forces P ci ; (3) Based on the aquifer gas storage model, the well working data is imported, the simulation time parameters are set, and the water saturation of each grid at different locations and times of the target gas storage model is obtained based on the differential solution of the control equation. S wi , capillary force P ci , permeability k i ; (4) Call the aquifer gas storage model in step (3) to calculate the gas-water interfacial tension in each grid at time t σ i , capillary force P ci , water saturation S wi , permeability k i , calculate the interface energy of each grid E i , the calculation formula is as follows: ; ； ； ； Where, E i For the i The interface energy of the grid, J; σ i is the air-water interfacial tension, mN / m; P ci is the capillary pressure in the i-th grid, Pa; S wi is the water saturation of the i-th grid, 0≤ S wi ≤1; is the water saturation gradient amplitude, m -1 ; x 、 y 、 z They correspond to the three orthogonal directions of the grid; is the volume of the grid, m³; X, Y, Z are the length, width, and height of the grid, respectively, in m; is the correction factor of the i-th grid; is the permeability of the i-th grid, mD; is the average permeability, mD; 、 are the viscosities of water and gas, mPa·S; (5) After obtaining the interface energy of each grid Then, the interface energy gradient amplitude of each grid is calculated : ； ； ； ； Where (x, y, z) is the grid index; x 、 y 、 z Three orthogonal directions of the grid respectively; is the grid spacing, m; 、 are the interface energies at the grid (x+1, y, z) and the grid (x-1, y, z), J; is the interface energy gradient amplitude of the ith grid, J / m; (6) The maximum value of the calculated interface energy gradient amplitude is recorded as ,Will 10% as the threshold , the interface energy gradient amplitude is greater than or equal to the threshold The grid is marked as the displacement front grid; (7) Calculate the fractal dimension of the displacement front, evaluate the front complexity based on the fractal dimension, and control the gas injection parameters in different zones.

2. The method for identifying the displacement front of an aquifer gas storage reservoir according to claim 1, characterized in that: In step (1), the empirical formula for capillary force, permeability and porosity is: ； ； Where, is the capillary force, Pa; is the air-water interfacial tension, ; is the porosity of the water layer; is the water layer permeability, ; is the dimensionless function of water saturation; is water saturation; is the irreducible water saturation; is the porosity distribution index.

3. The method for identifying the displacement front of an aquifer gas storage reservoir according to claim 1, characterized in that: In step (2), the aquifer gas storage model includes an injection well model, a drainage well model, a layer velocity model, a layer model, a porosity model and a permeability model; the construction process is as follows: (2-1) Import the wellhead coordinate data, well trajectory and logging data of gas injection well a and drainage well b into Petrel software to build models of gas injection well and drainage well; (2-2) Import 3D seismic volume and well-seismic calibration data, combine them with seismic data, and establish an interval velocity model; (2-3) Use seismic interpretation to build a model, use ant volume tracking for fault interpretation, use seed point controlled seismic volume tracking for horizon interpretation, and generate a layer model; (2-4) Establish porosity model and permeability model based on seismic attribute modeling.

4. The method for identifying displacement fronts during construction of an aquifer gas storage facility according to claim 1, wherein: The control equations include motion equations, mass conservation equations and state equations.

5. The method for identifying displacement fronts during construction of an aquifer gas storage facility according to claim 1, characterized in that: The specific steps of step (7) are: (7-1) The displacement front grid obtained by numerical simulation is discretized into a point matrix and covered with "coarse-fine" two-scale boxes; the number of boxes required to cover the front using the coarse grid is , using a fine grid to count the number of boxes needed to cover the leading edge is After obtaining the number of boxes required for the two-scale grid coverage, the fractal dimension method is used to reflect the complexity growth rate under the scale change of the displacement front, which is calculated as follows: ； Where, For the i The fractal dimension of a region, i =1, 2, 3…100; For the i The number of boxes under the coarse grid statistics of each region; For the i The number of boxes under the fine grid statistics of each area; R is the side length of the coarse square grid, m; r is the side length of the fine square grid, m; (7-2) Count the values of all leading edge points D f , calculate the quartiles; Q 1 is the 25th percentile, which is the boundary of the more regular area of the front edge; Q 3 is the 75th percentile, which is the boundary of the more fragmented area at the front edge; Stable zone: D f ≤Q 1; Transition Zone: Q 1< D f ＜Q 3; Risk areas: D f ≥Q 3; By comparing the calculated quantiles with the boundaries of each partition, the regularity of the gas-water front at that moment can be obtained.

Images