Wave velocity-attenuation full-flow calculation method based on Lo two-phase fluid model

Abstract

The invention belongs to the technical field of rock physics and multiphase flow coupling seismic wave propagation simulation, and discloses a wave velocity-attenuation whole-process calculation method based on an Lo two-phase fluid model. According to the method, thermodynamic characteristics of CO2-saline two-phase fluid are obtained by utilizing temperature and pressure conditions of a target reservoir, seepage characteristics are fitted according to experimental data, and a unified medium physical property layer is constructed by combining fluid parameters and rock skeleton parameters; after an elastic matrix, an inertia matrix and a viscosity matrix of a multiphase medium are calculated according to physical property parameters, phase velocity and attenuation characteristics of three types of longitudinal waves and one type of transverse waves are obtained through a Lo two-phase fluid model control equation. According to the method, a complete technical process from physical property construction, control equation solution to sensitivity judgment and mechanism analysis is formed, the method can be used for seismic monitoring, reservoir feature recognition and relevant parameter inversion in the CO2 geological storage process, and the method has good engineering applicability and popularization value.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of rock physics and coupled seismic wave propagation simulation of multiphase flow, and particularly relates to a wave velocity-attenuation full-process calculation method based on a Lo two-phase fluid model, a multiphase medium seismic wave propagation modeling method under the action of supercritical CO2-saline water two-phase fluid in a saline layer, and a numerical simulation and parameter sensitivity analysis method of seismic attributes such as wave velocity and attenuation. BACKGROUND

[0002] Geological carbon storage (GCS) is a key technology path to achieve global deep decarbonization. Large-scale deployment of GCS is crucial to achieve the urgent climate goal of "greenhouse gas emissions peaking in 2025 and reducing by 43% by 2030". Among various potential reservoirs, deep saline aquifers account for more than 90% of the total global geological storage potential. Due to its huge storage potential, it is considered as the primary target reservoir for implementing GCS. In saline aquifers, injected scCO2 will coexist with in-situ saline water, forming a dynamically evolving saline water-scCO2 two-phase fluid system. The spatial distribution and saturation change of this two-phase fluid will change the elastic properties of the reservoir rock (such as wave impedance, velocity), which will directly affect the amplitude, phase and propagation characteristics of the seismic reflection wave.

[0003] Biot theory laid the foundation for the theory of seismic wave propagation in saturated porous media. This theory couples the motion of the solid skeleton with the pore fluid and predicts the existence of type II compression waves. However, since the GCS target layer is considered a two-fluid system with saline water and scCO2 coexisting in the pores, the wave propagation characteristics within it are more complex: the difference in compressibility between the different fluid phases can induce local pore pressure gradients, thereby causing local flow and diffusion of the fluid in the pores. This phenomenon of waves causing local fluid flow is considered an important mechanism for seismic wave energy loss and velocity dispersion in porous two-fluid media. Subsequently, many models were developed to describe porous media saturated with two immiscible fluids. The Berryman-Thigpen-Chin (BTC) model, within the Biot framework, introduces a second type of fluid, considering the low-frequency limiting response and inertial coupling of the two-fluid system. However, this model neglects capillary pressure and severely underestimates the wave velocity dispersion and attenuation caused by capillary effects in the CO2-saltwater system. The Santos-Carcione model, based on the BTC model, introduces capillary pressure and improves the inertial coupling term between the two phases, enabling it to predict the third type of longitudinal wave (P3) generated by capillary pressure. However, the parameters of this model are mostly empirical or generalized assumptions, lacking sensitivity to pore geometry, and the inertial coupling coefficients are difficult to directly correspond to measured physical quantities. The Coussy model, through coupled pore elasticity theory, effective stress principle, and capillary effect, is mainly suitable for describing quasi-static seepage-deformation processes (such as soil consolidation and settlement), but is difficult to directly apply to consolidated rock systems. The Lo two-phase fluid theory is based on the mixture theory within the Euler framework. It establishes mass and momentum conservation equations for each phase (skeleton, wetting phase, and non-wetting phase) and clarifies the dynamic responses between the coupled phases through viscous drag and inertial interaction coefficients. Compared to the simplified treatment of capillary action and coupling parameters in previous BTC and Santos-Carcione models, the Lo model introduces a capillary pressure-saturation relationship and defines physically interpretable flow-flow and fluid-structure coupling parameters, enabling it to accurately predict three-branch body wave modes (including fast waves and two slow waves) and their corresponding dispersion and attenuation characteristics under low-frequency conditions. The Lo model provides crucial theoretical support for seismic wave monitoring of the saline-scCO2 two-phase fluid system in the GCS. By quantifying the wave velocity-attenuation coupling mechanism between water and scCO2, it can accurately analyze the wavefield response dynamics during CO2 seismic storage, providing a theoretical basis for elucidating the physical mechanism of seismic wave response changes induced by the two-phase fluid in saline CO2 seismic storage.

[0004] Based on the above analysis, the problems and shortcomings of the existing technology are as follows: (1) Current researches mostly use traditional van Genuchten (vG) model to represent the capillary pressure-saturation relationship, and the fitting accuracy is insufficient in the medium-low CO2 saturation interval; at the same time, the pore fluid parameters and rock skeleton parameters are often not obtained or constructed under the same temperature and pressure conditions, resulting in a lack of consistency and internal coupling relationship of medium physical property layers, which is difficult to meet the requirements of high consistency parameter system for two-phase fluid seismic simulation.

[0005] (2) Although the existing mixed or extended Biot model (such as BTC model and Santos-Carcione model) can describe the wave response of the two-fluid system, the coupling terms are mostly set empirically or simply, and it is difficult to realize the parameterization of physical interpretable parameters in terms of inertia coupling, viscous damping and capillary pressure, which cannot fully reflect the real interaction between CO2-saline two-phase fluid and solid skeleton.

[0006] (3) The existing methods mostly only output the calculation results of wave velocity and attenuation, and lack a special sensitivity analysis framework for the influence degree of different physical property parameter changes on the calculation results, which cannot clearly distinguish the relative influence of each parameter in the change of wave velocity and attenuation.

[0007] In summary, the existing technology lacks a complete seismic wave modeling and analysis method for the deep saline layer CO2 storage scene, which can integrate high-precision physical property input, clearly based on physical mechanism two-phase wave propagation model, and has systematic parameter influence analysis and mechanism interpretation ability. SUMMARY

[0008] In order to solve the problems of existing two-phase saturated rock seismic model parameters not easy to correspond physical quantities, capillary pressure curve error leading to unstable wave velocity prediction, and unclear different wave type response mechanism, the present disclosure provides a seismic wave propagation overall modeling scheme suitable for deep saline layer CO2 storage situation, which can form a continuous, calculable and derivable seismic physical simulation process from model parameter construction, two-phase fluid model wave velocity and attenuation solution, to key parameter sensitivity quantification. The present disclosure is implemented as follows based on the wave velocity-attenuation full process calculation method of Lo two-phase fluid model, including the following steps: S1, medium physical property construction: according to the temperature and pressure conditions of the target reservoir, the thermodynamic parameters of CO2-saline are obtained, and the capillary pressure-saturation relationship seepage characteristics are obtained by fitting formula according to experimental data; the two-phase fluid parameters and rock skeleton parameters are combined to construct a unified medium physical property layer; wherein the capillary pressure-saturation relationship is represented by an extended van Genuchten model to correct the capillary pressure fitting error in the medium-low saturation interval; S2, wave velocity and attenuation calculation: the elastic, inertial, viscous and shear modulus matrix is calculated by the medium property layer parameters, and substituted into the control equation of Lo two-phase fluid model to solve, to obtain the phase velocity and attenuation coefficient of three types of longitudinal waves and one type of transverse wave in CO2-saline saturated porous medium; S3, sensitivity and mechanism analysis: based on the wave velocity and attenuation calculation results, local and global sensitivity analysis is performed on each type of physical property parameter to identify the key parameters affecting the wave velocity; wherein, P1 wave is implemented three-channel decomposition, the influence of parameters on wave velocity is split into three parts of inertial channel, viscous channel and elastic channel, to clarify the sensitivity composition under different physical mechanisms.

[0009] In step S1, the thermodynamic parameters of CO2-saline include density, sound speed and dynamic viscosity; wherein, the density and sound speed of CO2 are calculated by Span-Wagner state equation, the expression is: ; In the formula, is pressure, unit is MPa; is CO2 density, is gas constant; is temperature, unit is K; is reduced density, , is critical density, ; is CO2 sound speed, is reduced temperature, , is critical temperature; is dimensionless Helmholtz free energy of ideal gas part, is dimensionless Helmholtz free energy of residual part, is first-order partial derivative of dimensionless Helmholtz free energy of residual part with respect to reduced density, ; is second-order mixed partial derivative of dimensionless Helmholtz free energy of residual part with respect to reduced density and reduced temperature, ; is second-order partial derivative of dimensionless Helmholtz free energy of ideal gas part with respect to reduced temperature, ; is second-order partial derivative of dimensionless Helmholtz free energy of residual part with respect to reduced density, ; The dynamic viscosity of CO2 is calculated by Laesecke-Muzny model: ; In the formula, viscosity of CO2, zero-density limit viscosity, linear density term, residual viscosity term, critical enhancement term; The density, sound velocity and dynamic viscosity of the brine are calculated using the Batzle-Wang empirical formula, which is expressed as: ; In the formula, is the density of the brine, is the sound velocity of the brine, is the sound velocity of pure water, , is an experimental fitting parameter; is the mass fraction of NaCl, is the viscosity of the brine.

[0010] In step S1, after obtaining the thermodynamic parameters of CO2-brine, the relative permeability model is used to describe the seepage capacity of CO2-brine two-phase in the pore medium, which is expressed as: ; In the formula, is the relative permeability of CO2, is the relative permeability of the brine, is the CO2 saturation, are all experimental fitting parameters.

[0011] In step S1, the capillary pressure-saturation relationship is characterized by the extended van Genuchten model to correct the capillary pressure fitting error in the medium-low saturation interval, including: The correction term takes as the core shape function, combined with the linear weight coefficient Without changing the physical meaning of the original vG parameters, the curve shape in the medium saturation interval is locally fine-tuned, and the expression of the extended van Genuchten capillary pressure model is: ; In the formula, is the capillary pressure, is the density of the brine, is the acceleration of gravity, are all experimental data fitting parameters; The skeleton parameters are composed of bulk modulus, shear modulus, porosity and permeability, which are used to represent the basic physical properties of the solid framework of the reservoir. The skeleton and fluid parameters are used as the total physical property input for seismic wave propagation calculation.

[0012] In step S2, after the medium property layer is constructed and all parameters of the rock skeleton and pore fluid are obtained, the Lo two-phase fluid model is used to calculate the propagation characteristics of different wave modes in the CO2-salt water system; the model unifies the elastic response of the solid framework, the inertial effect and viscous dissipation of the pore fluid, and the interface coupling effect caused by the capillary pressure gradient in the wave equation, and the phase velocity, dispersion relationship and attenuation coefficient of three types of longitudinal waves P1, P2, P3 and one type of transverse wave S wave are calculated by solving the propagation matrix in the form of complex eigenvalue, and the response law varying with frequency, saturation and pore structure is described.

[0013] Further, the Lo two-phase fluid model is based on Euler description, and the solid skeleton, non-wetting phase fluid and wetting phase fluid are regarded as a three-phase coexistence system, the motion behavior of each phase is described by mass conservation and momentum conservation equations, and the elastic response of the solid phase is constructed by combining the linear stress-strain relationship; the inertial coupling between the solid phase and the flow phase, the cross-viscous resistance between the flow phases and the capillary pressure gradient term are introduced, so that the deformation of the solid skeleton, the inertial effect of the pore fluid, the viscous dissipation and the interface tension process are uniformly represented; The three-phase medium constitutive equation is: ; In the formula, is the medium displacement, , is the solid phase, is the non-wetting flow phase, is the wetting flow phase, is the vector differential operator, which can be expressed as , is the first-order time derivative of the medium displacement, is the second-order time derivative of the medium displacement, is the inertial coupling matrix, which describes the virtual mass effect caused by the density of each phase and the acceleration of different phases, including the inertial coupling between the solid phase and each flow phase and the cross-inertial coupling between the flow phases; is the viscous coupling matrix, which describes the Darcy resistance and Yuster type cross-viscous resistance; is the elastic coefficient matrix, is the shear modulus coefficient matrix; The expressions of the inertial coupling matrix, the viscous coupling matrix and the elastic coefficient matrix are: ; ; ; ; In the formula, is Density of phase, ; is Volume fraction of phase, ; are constitutive coefficients related to inertial coupling between solid and fluid phases, respectively, are inertial cross-coupling coefficients between two fluids, respectively, are constitutive coefficients related to viscous coupling between solid and fluid phases, respectively, are cross-coupling coefficients between two fluids, respectively, is elastic coefficient, cross terms are symmetric, , ; is shear modulus of porous skeleton; Elastic coefficient is expressed as: ; ; ; ; wherein, is Bulk modulus of phase, ; is porosity of rock skeleton, is a dimensionless parameter describing porosity closure, obtained by uncased and cased experiments, ; are first and second water storage coefficients, respectively, is capillary pressure, is CO2 saturation; then we have: ; ; ; wherein, is bulk modulus of rock skeleton; is obtained according to the extended vG model: ; wherein, is rate of change of capillary pressure with respect to saturation; Inertial coefficient expression is: ; wherein, is pore tortuosity factor, for a random system with uniform circular pores of all orientations, the theoretical value is 3; when the pores are uniform and parallel to the axis and pressure gradient, the value is 1; when the solid particles of the porous medium are spherical, ; The viscous parameter is expressed as: ; ; is expressed as: ; In the formula, is The dynamic viscosity of the phase, ; is The relative permeability of the phase, ; is the intrinsic permeability; The divergence and the curl are obtained by using the three-phase medium constitutive equation, and the control equations of the longitudinal wave and the transverse wave are obtained: ; ; ; ; The general form of the solution of the three-phase medium constitutive equation is introduced: ; ; In the formula, is The longitudinal displacement of the phase, is The longitudinal wave amplitude of the phase, is The transverse displacement of the phase, is The transverse wave amplitude of the phase, is the complex wave number, is the propagation direction, is the angular frequency, is the time variable; After the general form of the solution of the three-phase medium constitutive equation is substituted into the control equations of the longitudinal wave and the transverse wave, the final control equations of the longitudinal wave and the transverse wave are obtained; The coefficient determinant is 0, so the control equation has a trivial solution: ; ; After expansion, it is a polynomial of The longitudinal wave control equation has six roots, and the transverse wave control equation has two roots. Under the physical constraint that the wave always decreases along the propagation direction, three types of longitudinal waves and one type of transverse wave are obtained.

[0014] In step S3, based on the wave velocity and attenuation calculation results, local and global sensitivity analyses are performed on various physical property parameters to identify key parameters affecting wave velocity, including: Sensitivity analysis was performed at the global level using the PAWN method; for each output parameter... Divide the entire interval into There are 10 sub-intervals, each sub-interval has a fixed number of sub-intervals. The remaining parameters are sampled in the entire space using the Latin hypercube sampling method to obtain the output. Conditional distribution The overall influence of each input parameter is quantified by calculating the Kolmogorov-Smirnov statistic between the conditional and unconditional distributions of the output. ,but: ; In the formula, In the first Given a given interval, the maximum absolute difference between the conditional cumulative distribution function and the unconditional cumulative distribution function. To take the upper bound for all possible inputs, To divide according to the range of values ​​of the input parameters After the first sub-interval, the second... A range; In order to meet the conditions Falling in each interval The conditional cumulative distribution function under the given conditions, It is the unconditional accumulation distribution function; For each parameter By combining the KS statistics of all sub-intervals, a global sensitivity index of PAWN is defined. : ; In the formula, To all By performing summary statistics, the final PAWN sensitivity index is obtained. .

[0015] Local sensitivity analysis is performed on the wave velocity and attenuation results to characterize the instantaneous response of the model to physical property disturbances under typical operating conditions. A disturbance is applied to each input parameter near a given set of reference parameters, and the difference between the output before and after the disturbance is calculated to obtain the first-order sensitivity of the output with respect to the parameter. For output variables With input parameters The local sensitivity coefficient is expressed as: ; In the formula, is a local sensitivity parameter, is a medium parameter for the local sensitivity parameter. In step S3, three-channel decomposition is performed on the P1 wave, including: By using the three-channel decomposition method, the change of the eigenvalue corresponding to the P1 wave is divided into the contributions of the elastic-capillary, inertial and viscous three physical channels, so as to identify the dominant path controlling the wave velocity change; The critical saturation at which the P1 wave changes from descending to ascending is determined by the velocity turning point solving method, and the dependence of the turning point on each physical parameter is analyzed; Under the frequency domain eigenvalue problem of the Lo model, the P1 wave velocity The first-order sensitivity of is written as the sum of the inertial, elastic-capillary and viscous three channels; the characteristic equation is rewritten as a generalized eigenvalue problem, and the expression is: ; The left eigenvalue is introduced, and the normalization condition is applied, and the eigenvalue is written as: ; Under the above normalization condition, the derivative of the mechanism variable is obtained, and taking the porosity as an example, the left multiplication is performed to obtain the first-order derivative of the eigenvalue: ; For the body wave, at a fixed frequency, there is: ; The first-order sensitivity of the velocity with respect to the porosity is obtained, and is decomposed into the sum of the elastic-inertial-viscous three channels: ; Wherein: ; ; ; In the formula, is the elastic-capillary channel, , is the inertial channel, is the viscous channel.

[0016] In combination with all the above technical solutions, the application has the beneficial effects of: First, the present application is around the problem of modeling and interpretation of seismic response in CO2-saline two-phase medium, and a complete technical process covering physical property construction, wave velocity-attenuation solution and sensitivity mechanism analysis is constructed and implemented. The Lo two-phase fluid model has the advantages of strong consistency and strict physical constraints in describing the coupling dynamics of multi-phase fluid-solid, but it usually needs to be processed step by step for parameter changes, wave field response and mechanism interpretation in the application process, which is easy to cause information fragmentation. Based on this situation, the present application expands the traditional modeling framework, unifies the temperature-pressure related fluid properties, relative permeability model and improved capillary pressure model into the medium physical layer, and realizes the accurate calculation of phase velocity, dispersion and attenuation of various wave types under two-phase conditions on this basis. Further, the present application carries out sensitivity analysis on key physical parameters, and based on the inertial, viscous and elastic-capillary three-channel decomposition method, it determines the overall modeling process from parameter perturbation to wave field response to physical mechanism, thereby solving the problem of information fragmentation and mechanism quantification in traditional methods.

[0017] Second, the present application verifies the process by using typical reservoir properties and CO2-saline data, and the results show that the method can stably calculate the responses of three types of P waves and one type of S wave under different conditions, correctly depict the key features in the velocity-saturation relationship, and effectively identify the influence direction and contribution mechanism of main control parameter changes on seismic attributes. The whole process has clear structure, strong physical interpretation and high calculation efficiency, and can be directly used for seismic monitoring, reservoir state identification and parameter inversion in the process of CO2 storage, which has good engineering applicability and popularization value.

[0018] Third, the CO2-saline two-phase medium modeling process proposed by the present application can be used in CO2 geological storage (GCS) engineering to quantitatively analyze the velocity changes caused by fluid saturation changes and fluid phase state changes in the process of CO2 injection, so as to help engineers master the disturbance amplitude and evolution characteristics of the reservoir velocity structure in the injection process. The process can be used as a physical modeling tool in carbon storage monitoring and evaluation, and can provide reference information with physical basis for injection safety evaluation and monitoring scheme optimization.

[0019] The existing multi-phase medium seismic simulation is mostly dispersed processing, such as setting physical properties, solving velocity or local parameter analysis, and lacks a complete technical route that can penetrate modeling, solving and sensitivity analysis. The present application unifies the temperature-pressure related fluid parameters, relative permeability and capillary pressure into a multi-phase physical layer under the Lo multi-phase coupling model framework, and solves the velocity and attenuation of multiple wave types according to this; on this basis, further sensitivity analysis and corresponding inertial-viscous-elastic capillary mechanism decomposition are carried out, so that the parameter influence path is revealed systematically. The whole chain of "physical property -> solving -> sensitivity -> mechanism" constructed by the present application makes up for the lack of systematic process of multi-phase medium seismic simulation.

[0020] In the multi-phase medium, the coupling relationship of physical parameters is complex, the wave type is multiple and interferes with each other, so that it is difficult to realize from parameter change to velocity and attenuation response and then to the division of influencing factors. The present application makes the action path of each parameter on velocity and attenuation under the condition of multi-phase coupling clear through unified physical modeling, stable solving of multiple wave types and supporting sensitivity analysis and mechanism decomposition, so as to solve the problem of "difficult to quantitatively distinguish the source of parameter influence" in the multi-phase system.

[0021] The existing multi-phase medium seismic simulation method is usually difficult to finely distinguish the contribution of each physical parameter and its coupling term when analyzing the influencing factors of velocity and attenuation, especially when the capillary pressure and the joint action of multiple parameters are involved, the influence path is often difficult to be clearly presented in the same framework. The present application makes the influence source of each parameter on velocity and attenuation be able to be identified in a clearer and quantitative way through systematic sensitivity analysis and mechanism decomposition of inertial, viscous and elastic-capillary action paths, so as to overcome the lack of refinement of traditional methods in parameter influence analysis. BRIEF DESCRIPTION OF DRAWINGS

[0022] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the present disclosure and serve to explain the principles of the present disclosure, together with the description; Figure 1 is a flow chart of the wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model provided by the embodiment of the present application; Figure 2 is a route map of the wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model provided by the embodiment of the present application; Figure 3FIG. 1 is a schematic view of the variation of the CO2-saline water physical property parameters under different temperature and pressure conditions provided by the embodiments of the present application; wherein, (a) is the density variation with temperature (P = 15 MPa), (b) is the longitudinal wave velocity variation with temperature (P = 15 MPa), (c) is the viscosity variation with temperature (P = 15 MPa), (d) is the density variation with pressure (T = 40℃), (e) is the longitudinal wave velocity variation with pressure (T = 40℃), and (f) is the viscosity variation with pressure (T = 40℃); Figure 4 FIG. 2 is a graph of experimental data and fitting curves provided by the embodiments of the present application; wherein, (a) is the Pini (2012) experimental data, the traditional vG model and the Ext-vG model fitting curves, and (b) is the Wang (2016) experimental data and fitting curves; Figure 5 FIG. 3 is a schematic view of the relationship between the wave velocity and attenuation and the CO2 saturation and frequency provided by the embodiments of the present application; wherein, (a) is the P1 wave velocity, (b) is the P1 wave attenuation, (c) is the P2 wave velocity, (d) is the P2 wave attenuation, (e) is the P3 wave velocity, (f) is the P3 wave attenuation, (g) is the S wave velocity, and (h) is the S wave attenuation; Figure 6 FIG. 4 is a schematic view of the global sensitivity analysis of the wave velocity provided by the embodiments of the present application; wherein, (a) is the P1 wave, (b) is the P2 wave, (c) is the P3 wave, and (d) is the S wave; Figure 7 FIG. 5 is a schematic view of the local sensitivity analysis of the wave velocity provided by the embodiments of the present application; wherein, (a) is the P1 wave, (b) is the P2 wave, (c) is the P3 wave, and (d) is the S wave; Figure 8 FIG. 6 is a sensitivity decomposition graph of the P1 wave to the porosity and the two-fluid density provided by the embodiments of the present application; wherein, (a) is the porosity, (b) is the CO2 density, and (c) is the salt water density; Figure 9 FIG. 7 is a saturation dependence schematic view of the P1 wave velocity and the elastic modulus and the density sensitivity in the CO2-saline water-sandstone system provided by the embodiments of the present application; wherein, (a) is the P1 wave velocity and the turning point, and (b) is the modulus sensitivity and the density sensitivity intersection graph. DETAILED DESCRIPTION

[0023] In order to make the above objectives, characteristics and advantages of the present application more obvious and easy to understand, the specific embodiments of the present application are described in detail below with reference to the drawings. In the following description, a large number of specific details are set forth in order to provide a thorough understanding of the present application. However, the present application can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without departing from the spirit of the present application, so the present application is not limited to the specific implementations disclosed below.

[0024] The innovation of the present application is: (1) The extended van Genuchten capillary pressure model is proposed, which improves the fitting accuracy of the capillary pressure-saturation relationship in the medium-low saturation section, and provides a more reliable seepage parameter basis for multiphase medium modeling.

[0025] (2) The CO2-saline two-phase system physical property model is constructed to meet the temperature and pressure conditions of the reservoir, and the fluid thermodynamic parameters, relative permeability and capillary pressure are integrated into the multiphase physical layer to realize the physical consistency of the physical property input under the two-fluid condition.

[0026] (3) The turning point feature of P1 wave velocity changing with saturation is revealed, and it is clear that the turning point feature is caused by the intersection of effective bulk modulus and effective density, which provides a physical basis for explaining the velocity inversion or abnormal change during injection.

[0027] (4) A multi-wave type sensitivity analysis system is established, and parameter sensitivity analysis and mechanism decomposition are carried out on P1, P2 and P3 waves respectively, so as to clearly distinguish the influence of inertia, viscosity and elastic-capillary on velocity and attenuation, and to provide a systematic method for mechanism interpretation of multiphase wave response.

[0028] In the embodiments of the present application, the wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model includes the following steps: Figure 1 and Figure 2 In the embodiments of the present application, the wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model includes the following steps: S1, medium property construction: according to the temperature and pressure conditions of the target reservoir, the thermodynamic parameters of CO2-saline are obtained, and the capillary pressure-saturation relationship seepage characteristics are obtained by fitting formula according to experimental data; the two-phase fluid parameters and rock skeleton parameters are combined to construct a unified medium property layer; wherein the capillary pressure-saturation relationship is characterized by the extended van Genuchten model to correct the capillary pressure fitting error in the medium-low saturation interval; (1) After inputting given temperature and pressure, the density and sound velocity of CO2 are calculated by using Span and Wagner generalized equation of state (formula 1), and the dynamic viscosity of CO2 is calculated by using Laesecke and Muzny model (formula 2). The density, longitudinal wave sound velocity and dynamic viscosity of brine (NaCl-H2O system) are calculated by using the empirical relationship formula proposed by Batzle and Wang (formula 3).

[0029] The thermodynamic parameters of the CO2-brine include density, sound velocity and dynamic viscosity; wherein the density and sound velocity of CO2 are calculated by using Span-Wagner equation of state, and the expression is as follows: ; In the formula, is pressure, unit: MPa; is CO2 density, is gas constant; is temperature, unit: K; is reduced density, , is critical density, ; is CO2 sound velocity, is reduced temperature, , is critical temperature; is dimensionless Helmholtz free energy of ideal gas part, is dimensionless Helmholtz free energy of residual part, is first-order partial derivative of dimensionless Helmholtz free energy of residual part with respect to reduced density, ; is second-order mixed partial derivative of dimensionless Helmholtz free energy of residual part with respect to reduced density and reduced temperature, ; is second-order partial derivative of dimensionless Helmholtz free energy of ideal gas part with respect to reduced temperature, ; is second-order partial derivative of dimensionless Helmholtz free energy of residual part with respect to reduced density, ; The dynamic viscosity of CO2 is calculated by using Laesecke-Muzny model: ; In the formula, is CO2 viscosity, is zero-density limit viscosity, is linear density term, is remaining viscosity term, is the critical enhancement term; The density, sound speed and dynamic viscosity of brine are calculated by Batzle-Wang empirical formula, which is expressed as: (3); wherein, is the density of brine, is the sound speed of brine, is the sound speed of pure water, , is the experimental fitting parameter; is the mass fraction of NaCl, is the viscosity of brine.

[0030] Figure 3 The above formula is obtained under different temperature (30-60℃) and pressure (5-35MPa) conditions. Among them, Figure 3 (a) in the figure is the density change with temperature (P=15MPa), Figure 3 (b) in the figure is the longitudinal wave velocity change with temperature (P=15MPa), Figure 3 (c) in the figure is the viscosity change with temperature (P=15MPa), Figure 3 (d) in the figure is the density change with pressure (T=40℃), Figure 3 (e) in the figure is the longitudinal wave velocity change with pressure (T=40℃), Figure 3 (f) in the figure is the viscosity change with pressure (T=40℃).

[0031] (2) Capillary pressure and relative permeability parameter calibration: input capillary pressure and relative permeability experimental data, and fit the model parameters. Capillary pressure uses Pini (2012) Berea sandstone experimental data at 50℃, 9MPa, which is respectively substituted into the traditional vG capillary pressure equation and Ext-vG capillary pressure equation for fitting. The vG model fitting parameters are: and . The Ext-vG model fitting parameters are: , , , ; the relative permeability uses the CO2-brine two-phase relative permeability experimental data given by Wang (2016), which is substituted into the relative permeability equation for fitting. The fitting parameters are , and .

[0032] After obtaining the thermodynamic parameters of CO2-brine, the relative permeability model is used to describe the seepage capacity of CO2-brine two-phase in porous media, which is expressed as: ; where, is the relative permeability of CO2, is the relative permeability of brine, is the CO2 saturation, are experimental fitting parameters.

[0033] The capillary pressure-saturation relationship is characterized by an extended van Genuchten model to correct the capillary pressure fitting error in the medium-low saturation interval, including: The correction term is with the core shape function combined with the linear weight coefficient Without changing the physical meaning of the original vG parameters, the curve shape in the medium saturation interval is locally fine-tuned, and the expression of the extended van Genuchten capillary pressure model is: ; ; where, is the capillary pressure, is the brine density, is the gravitational acceleration, are experimental data fitting parameters; The skeleton parameters are composed of bulk modulus, shear modulus, porosity and permeability, which are used to characterize the basic physical properties of the solid framework of the reservoir, and the skeleton and fluid parameters are all physical properties input for seismic wave propagation calculation.

[0034] The experimental data and fitting results are shown in Figure 4 , wherein, Figure 4 (a) of the figure is the Pini (2012) experimental data, the traditional vG model and the Ext-vG model fitting curve, Figure 4 (b) of the figure is the Wang (2016) experimental data and the fitting curve.

[0035] (3) Rock parameter input: The present application uses the typical rock physical property parameters of Berea sandstone for numerical experiment, selects the mineral matrix bulk modulus 39.69 GPa, shear modulus 37.22 GPa, dry rock frame bulk modulus 11.6 GPa, shear modulus 11.1 GPa, porosity 0.19, permeability 2.72×10 -13 m 2 As the input parameters of the rock skeleton in the two-phase fluid-solid system.

[0036] S2, wave velocity and attenuation calculation: the elastic, inertial and viscous matrices are calculated by the medium property layer parameters, and substituted into the control equations of Lo two-phase fluid model to solve, and the phase velocity and attenuation coefficient of three types of longitudinal waves and one type of transverse wave in CO2-saline saturated porous medium are obtained; (1) define the excitation frequency and calculate the elastic, viscous and inertial matrices: at the set excitation frequency, the rock skeleton parameters and pore fluid properties provided by the medium property layer are directly called to construct the elastic matrix, viscous matrix and inertial matrix of the multi-phase coupled system.

[0037] (2) solve the control equation to obtain the wave velocity and attenuation of four types: input the above matrix into the Lo model control equation, solve the frequency domain eigenvalue problem, and obtain the phase velocity, dispersion characteristics and attenuation coefficient of three types of compression waves and one type of shear wave.

[0038] The relationship between wave velocity and attenuation with CO2 saturation and frequency is shown in Figure 5 . Among them, Figure 5 (a) of (a) in (a) is the P1 wave velocity, Figure 5 (b) of (a) in (a) is the P1 wave attenuation, Figure 5 (c) of (c) in (c) is the P2 wave velocity, Figure 5 (d) of (c) in (c) is the P2 wave attenuation, Figure 5 (e) of (e) in (e) is the P3 wave velocity, Figure 5 (f) of (e) in (e) is the P3 wave attenuation, Figure 5 (g) of (g) in (g) is the S wave velocity, Figure 5 (h) of (g) in (g) is the S wave attenuation.

[0039] After the medium property layer is constructed and all the parameters of rock skeleton and pore fluid are obtained, the Lo two-phase fluid model is used to calculate the propagation characteristics of different wave types in CO2-saline system; the model unifies the elastic response of solid framework, the inertial effect and viscous dissipation of pore fluid, and the interface coupling effect caused by capillary pressure gradient in the wave equation, and the phase velocity, dispersion relationship and attenuation coefficient of three types of longitudinal waves P1, P2, P3 and one type of transverse wave S wave are calculated by solving the complex eigenvalue form of propagation matrix, and the response law with frequency, saturation and pore structure is described.

[0040] The Lo two-phase fluid model is based on Euler description, and the solid skeleton, non-wetting phase fluid and wetting phase fluid are regarded as a three-phase coexistence system, the motion behavior of each phase is described by mass conservation and momentum conservation equation, and the elastic response of solid phase is constructed by linear stress-strain relationship; the inertial coupling between solid phase and flow phase, cross-viscous resistance between flow phase and flow phase, and capillary pressure gradient term are introduced, so that the deformation of solid skeleton, inertial effect of pore fluid, viscous dissipation and interface tension process can be expressed uniformly; The constitutive equation of three-phase medium is: ; where, is the medium displacement, , is the solid phase, is the non-wetting fluid phase, is the wetting fluid phase, is the vector differential operator, which in Cartesian coordinates can be expressed as , is the first time derivative of the medium displacement, is the second time derivative of the medium displacement, is the inertial coupling matrix, which characterizes the inertial coupling between the solid and the fluids and the cross-inertial coupling between the fluids; is the viscous coupling matrix, which describes the Darcy resistance and the Yuster-type cross-viscous resistance; is the elastic coefficient matrix, is the shear modulus coefficient matrix; The expressions of the inertial coupling matrix, the viscous coupling matrix and the elastic coefficient matrix are: ; ; ; ; where, is the density of the phase, ; is the volume fraction of the phase, ; are the constitutive coefficients related to the inertial coupling between the solid and the fluids, are the cross-inertial coupling coefficients between the two fluids, are the constitutive coefficients related to the viscous coupling between the solid and the fluids, are the cross-viscous coupling coefficients between the two fluids, are the elastic coefficients, the cross terms are symmetric, , ; is the shear modulus of the porous skeleton; The elastic coefficients are expressed as: ; ; ; ; where, for The bulk modulus of the phase, ; For rock skeleton porosity, Dimensionless parameters describing porosity closure were obtained through both sleeveless and sleeved experiments. ; These are the first and second water storage coefficients, respectively. For capillary pressure, CO2 saturation; Then we have: ; ; ; In the formula, The bulk modulus of the rock skeleton; Obtained from the extended vG model: ; In the formula, The rate of change of capillary pressure with respect to saturation; The expression for the inertia coefficient is: ; In the formula, The pore tortuosity factor has a theoretical value of 3 for a random system with uniform circular pores of all orientations; a value of 1 when the pores are uniform and parallel to the axis and pressure gradient; and a value of 1 when the solid particles in the porous medium are spherical. ; The viscosity parameter is expressed as: ; ; Represented as: ; In the formula, for The dynamic viscosity of the phase, ; for The relative permeability of the phase, ; This is the inherent penetration rate; By using the constitutive equations of a three-phase medium to obtain the divergence and curl, the governing equations for longitudinal and transverse waves are derived: ; ; Introducing the general form of the solution to the constitutive equation for a three-phase medium: ; where, is the longitudinal displacement of the phase, is the longitudinal wave amplitude of the phase, is the transverse displacement of the phase, is the transverse wave amplitude of the phase, is the complex wave number, is the propagation direction, is the angular frequency, is the time variable; Substitute the general form of the three-phase medium constitutive equation solution into the control equations of longitudinal and transverse waves to obtain the final control equations of longitudinal and transverse waves; Set the coefficient determinant to 0, then the control equation has a trivial solution: ; After expansion, it is a polynomial of , the longitudinal wave control equation has six roots, and the transverse wave control equation has two roots. Under the physical constraint that the wave always decreases along the propagation direction, there are three types of longitudinal waves and one type of transverse wave.

[0041] S3, sensitivity and mechanism analysis: based on the wave velocity and attenuation calculation results, local and global sensitivity analysis is performed on each type of physical parameter to identify the key parameters that affect the wave velocity; among them, three-channel decomposition is implemented for P1 wave, the influence of parameters on wave velocity is divided into three parts of inertia channel, viscous channel and elastic channel, and the sensitivity under different physical mechanisms is clarified.

[0042] (1) Use the global sensitivity method (such as PAWN) to calculate the sensitivity index of three types of P waves and S waves excited at 100 Hz under different parameter perturbations, obtain the overall influence degree of wave velocity on each input parameter, and form Figure 6 . Among them, Figure 6 the (a) figure in (a) is P1 wave, Figure 6 the (b) figure in (b) is P2 wave, Figure 6 the (c) figure in (c) is P3 wave, Figure 6 and the (d) figure in (d) is S wave.

[0043] (2) Under the condition of 100 Hz excitation, a 1% amplitude perturbation is applied to each input parameter centered on the reference state, the local partial derivative sensitivity of wave velocity with respect to porosity, permeability and fluid properties is calculated, the local response characteristics of three types of P waves and S waves are obtained, and Figure 7 is generated. Among them, Figure 7 the (a) figure in (a) is P1 wave, Figure 7 the (b) figure in (b) is P2 wave, Figure 7 the (c) figure in (c) is P3 wave,Figure 7 Figure (d) in the diagram represents the S-wave.

[0044] (3) Under 100Hz excitation conditions, porosity, CO2 density, and saline water density were selected as the three main control parameters for the P1 wave. The three-channel decomposition method was used to separate and compare the sensitivity contributions of the P1 wave to these three inputs, analyze the relative strength of the three channels under different saturation levels, and form Figure 8 .in, Figure 8 Figure (a) shows porosity. Figure 8 Figure (b) shows the CO2 density. Figure 8 Figure (c) shows the density of saline water.

[0045] (4) Under 100Hz excitation conditions, the sensitivity sequences of the equivalent bulk modulus channel and the effective density channel as a function of saturation are calculated respectively; then, a cross-hatching operation is performed on the two sensitivity sequences. When the sensitivity values ​​of the two channels are equal, the corresponding saturation value is recorded; this saturation is the inflection point saturation of the P1 wave velocity. The saturation dependence of P1 wave velocity and elastic modulus with density sensitivity in the CO2-saltwater-sandstone system is as follows: Figure 9 As shown. Among them, Figure 9 Figure (a) shows the P1 wave velocity and inflection point. Figure 9 Figure (b) in the figure is a cross plot of modulus sensitivity and density sensitivity.

[0046] Based on the wave velocity and attenuation calculation results, local and global sensitivity analyses were performed on various physical property parameters to identify key parameters affecting wave velocity, including: Sensitivity analysis was performed at the global level using the PAWN method; for each output parameter... Divide the entire interval into There are 10 sub-intervals, each sub-interval has a fixed number of sub-intervals. The remaining parameters are sampled in the entire space using the Latin hypercube sampling method to obtain the output. Conditional distribution The overall influence of each input parameter is quantified by calculating the Kolmogorov-Smirnov statistic between the conditional and unconditional distributions of the output. ,but: ; In the formula, In the first Given a given interval, the maximum absolute difference between the conditional cumulative distribution function and the unconditional cumulative distribution function. To take the upper bound for all possible inputs, To divide according to the range of values ​​of the input parameters After the first sub-interval, the second... interval; is the conditional cumulative distribution function of the parameter ; is the unconditional cumulative distribution function of the parameter ; is the unconditional cumulative distribution function of the parameter ; , the PAWN global sensitivity index is defined as : ; where is the summary statistic of all ; . The local sensitivity analysis of wave velocity and attenuation results is performed to characterize the instantaneous response of the model to the physical property perturbation under typical operating conditions. In the vicinity of a given set of reference parameters, a perturbation is imposed on each input parameter, and the difference between the outputs before and after the perturbation is calculated to obtain the first-order sensitivity of the output with respect to the parameter; For the output variable and the input parameter , the local sensitivity coefficient is expressed as: ; where is the local sensitivity parameter, is the medium parameter for local sensitivity analysis. Three-channel decomposition is performed on P1 wave, including: Using the three-channel decomposition method, the change in the eigenvalue corresponding to P1 wave is divided into the contributions of the elastic-capillary, inertial, and viscous three physical channels, to identify the dominant path that controls the wave velocity change; The critical saturation at which P1 wave changes from descending to ascending is determined by the velocity turning point method, and the dependence of this turning point on each physical property parameter is analyzed; Under the frequency-domain eigenproblem of the Lo model, the first-order sensitivity of P1 wave velocity with respect to is written as the sum of the inertial, elastic-capillary, and viscous three channels; the characteristic equation is rewritten as a generalized eigenvalue problem, expressed as: ; The left eigenvalue is introduced, and the normalization condition is applied, and the eigenvalue is written as: ; Under the above normalization condition, the derivative of the mechanism variable is taken, and taking porosity as an example, left-multiplying , the first-order derivative of the eigenvalue is obtained: ; For body waves, at a fixed frequency, we have: ; Obtain the velocity with respect to porosity The first-order sensitivity is decomposed into the sum of the elastic, inertial, and viscous channels: ; in: ; ; ; In the formula, It is an elastic capillary channel. , For inertial channels, It is a viscous channel.

[0047] After quantifying the contributions of the elastic-capillary, inertial, and viscous channels to the P1 wave velocity sensitivity using a three-channel decomposition method, the inflection point of the P1 wave velocity was determined, thereby identifying the critical saturation and controlling factors at which the velocity changes from decreasing to increasing. The approximate expression for the P1 wave velocity is as follows: ; In the formula, It is the equivalent bulk modulus of the medium, which is determined by the solid skeleton, fluid compressibility, and coupling effects. It is the equivalent density of the medium, obtained by the linear superposition of the contributions of each phase; the relationship between the equivalent bulk modulus and the equivalent density is analyzed. Relative sensitivity: ; In the formula, It is the equivalent bulk modulus pair Sensitivity It is an equivalent density pair The sensitivity of equivalent bulk modulus and equivalent density to saturation; when the relative rates of change of equivalent bulk modulus and equivalent density are equal, the sensitivity of equivalent bulk modulus to equivalent density is... The relative sensitivity expression gives a sensitivity term of zero, and the P1 wave velocity reaches a local minimum. This condition is the critical saturation criterion for the velocity to change from decreasing to increasing.

[0048] Example 2: The wave velocity-attenuation full-process calculation system based on the Lo two-phase fluid model provided in this embodiment of the invention includes: A medium property construction module is configured to construct a unified medium property layer according to temperature, pressure and rock skeleton parameters of a target reservoir, wherein the percolation parameters include a capillary pressure-saturation relationship, and the capillary pressure-saturation relationship is characterized by an extended van Genuchten capillary pressure model. A wave velocity and attenuation calculation module is in communication connection with the medium property construction module, configured to, after the medium property layer parameters are constructed into an elastic, viscous, inertial and shear modulus matrix, substitute the parameters into control equations of a Lo two-phase fluid model to solve the equations and calculate phase velocities and attenuation coefficients of three types of longitudinal waves and one type of transverse wave in a CO2-saline-saturated porous medium. A sensitivity and mechanism analysis module is in communication connection with the wave velocity and attenuation calculation module, configured to, based on the wave velocity and attenuation calculation results, perform sensitivity analysis on physical property parameters affecting the wave velocity and attenuation, wherein the sensitivity analysis includes three-channel decomposition analysis for P1 wave for physical mechanism interpretation.

[0049] Preferably, the medium property construction module includes: A fluid thermodynamic parameter calculation unit is configured to, by inputting only temperature and pressure of a work area, calculate density and sound velocity of CO2 based on a Span-Wagner state equation, calculate dynamic viscosity of CO2 based on a Laesecke-Muzny model, and calculate density, sound velocity and dynamic viscosity of saline water based on a Batzle-Wang empirical formula. A percolation parameter characterization unit is configured to determine a relative permeability-saturation relationship based on a Mualem-van Genuchten model and determine a capillary pressure-saturation relationship based on an extended van Genuchten capillary pressure model.

[0050] Preferably, the wave velocity and attenuation calculation module is configured to calculate phase velocities, attenuation coefficients and dispersion characteristics of three types of longitudinal waves (P1, P2 and P3) and one type of transverse wave (S wave) by the Lo two-phase fluid model.

[0051] Preferably, the sensitivity and mechanism analysis module includes: A global sensitivity analysis unit is configured to quantify overall influence degrees of each input parameter on wave velocity and attenuation results by using a PAWN method; A local sensitivity analysis unit is configured to calculate first-order sensitivities of an output quantity with respect to each input parameter in a vicinity of a reference parameter group; A special mechanism analysis unit is configured to perform three-channel decomposition analysis of P1 wave and velocity turning point calculation.

[0052] Example 3: Comparison of effects of a traditional model and an extended van Genuchten capillary pressure model. Figure 4(a) The least square fitting of the traditional vG model and the extended vG model is performed using the Pc-S experimental data of Berea sandstone based on Pini's report. The optimal parameters of the traditional vG model are , . The fitting parameters of the Ext-vG model are , , , . On this basis, the root mean square error (RMSE) of the traditional vG model is 0.812 kPa, and the determination coefficient R 2 = 0.9759; the RMSE of the extended vG model is 0.343 kPa, and the R 2 is 0.9957, and the fitting residual of the Pc-S curve in the low to moderate CO2 saturation interval is reduced as a whole. For this Berea sandstone data set, the extended vG model can thus more smoothly reproduce the trend of experimental capillary pressure with saturation, while still maintaining the asymptotic behavior consistent with the traditional vG model at the limits of CO2→0 and SCO2→1. In the above embodiments, the description of each embodiment is focused on, and the parts not described or recorded in a certain embodiment can be referred to the related description of other embodiments.

[0053] To further prove the positive effect of the above embodiments, the present application based on the above technical solutions carries out the following experiments.

[0054] (1) Model fitting effect: Ext-vG has obvious accuracy improvement compared with vG; According to the experimental capillary pressure-saturation data of Berea sandstone, as shown in Figure 4 , the fitting accuracy of the Ext-vG model proposed by the present application in the 0.1-0.6 saturation interval is obviously improved. The root mean square error (RMSE) of the traditional van Genuchten model is 0.812 kPa, while the Ext-vG model of the present application reduces the RMSE to 0.343 kPa; the determination coefficient R 2 from 0.9759 to 0.9957. From the fitting curve, it can be seen that the Ext-vG model is well fitted with the experimental data in the low to moderate saturation section, effectively avoiding the systematic deviation of the traditional vG model.

[0055] (2) Wave speed-attenuation variation law based on multiphase model; After calculating the speed and attenuation of P1, P2, P3 and S waves under multiphase conditions, it can be observed that various wave types present clear and coherent physical variation laws with saturation and frequency, as shown in Figure 5P1 wave speed presents a rapid decay after a small amount of CO2 injection followed by a slow decay with saturation increase under most temperature-pressure conditions; under partial temperature-pressure conditions, it shows a turning trend from decline to rise with saturation increase, and its minimum point can be clearly identified on the calculated speed-saturation curve; P2 wave presents a trend of rapid decay followed by slow decay; P3 wave appears slow wave characteristics under low frequency conditions, and its speed and attenuation have a typical diffusion-type response with saturation change; S wave is not sensitive to saturation change, and its speed curve is relatively stable. The solving process of the application can present a reasonable evolution relationship between multiple wave types, so that the change trend, relative difference and characteristic interval of P1, P2, P3 and S wave can be fully embodied in the calculation results.

[0056] (3) Sensitivity analysis: the main control pattern of four types of wave parameters is clear; In the given saturation range, the sensitivity of each wave type presents a clear main control pattern, as shown in Figure 6 , 7 : P1 wave is most sensitive to porosity and the densities of two types of fluids, and the sensitivity gradually increases with the increase of saturation; although the densities of CO2 and brine are opposite in direction, they both change linearly with saturation when the saturation is more than 10%, and the influence of brine density is more significant; the speed of sound only plays a role in the low saturation interval, and the influences of fluid viscosity and skeleton permeability are close to negligible. P2 wave is positively sensitive to skeleton permeability and CO2 speed, and negatively sensitive to porosity, and the sensitivity of each type converges rapidly at low saturation and at medium and high saturation; the sensitivity of CO2 and brine viscosity respectively presents a trend of negative enhancement and continuous decay with saturation, and tends to be stable at high saturation. P3 wave is positively sensitive to skeleton permeability and CO2 density, and negatively sensitive to porosity, and its response to fluid viscosity is different from that of P2 wave: the sensitivity of CO2 viscosity decreases with the increase of saturation, while the negative sensitivity of brine viscosity increases with the increase of saturation. The sensitivity trend of S wave to porosity and fluid density is similar to that of P1 wave, and it is not sensitive to other parameters. Global PAWN analysis further shows that the main control factors of each wave type have significant distinguishability: P1 wave is mainly dominated by porosity; P2 wave is mainly controlled by CO2 speed; P3 wave is most sensitive to saturation, and its dependence on CO2 density, porosity and brine viscosity changes at different saturation intervals; S wave is almost completely dominated by porosity. The above results collectively show that different wave types have stable and differentiated main control parameter patterns under multi-phase conditions. Based on three-channel decomposition of the perturbed parameters, the dominant mechanism at different saturation stages can be clearly observed to change systematically with saturation, as shown in Figure 8The total sensitivity is determined by the elastic-capillary channel and the inertial channel in the low saturation range, in which the elastic-capillary term has a significant negative contribution due to the large derivative of the capillary pressure; the channel rapidly decays and tends to zero after entering the transition zone, and the inertial channel begins to dominate, so that the change of the total sensitivity mainly comes from the inertial term; in the medium and high saturation ranges, the contributions of the elastic-capillary term and the inertial term gradually approach, and the total sensitivity presents a single-channel characteristic dominated by the inertial channel. The three-channel decomposition of CO2 and salt water density also shows a similar structure: the CO2 density mainly produces a negative contribution through the inertial channel, while the elastic-capillary term has a limited positive impact at low saturation and rapidly decays as the saturation increases; the salt water density presents a positive contribution under very low saturation conditions, and then its elastic-capillary term gradually weakens and tends to be constant, and finally the sensitivity is dominated by the inertial term. Overall, the three-channel decomposition not only ensures the convergence and self-consistency of the total sensitivity, but also reveals the continuous evolution law of the dominant path under low, medium and high saturation, providing a clear mechanism for understanding the physical response of wave speed to parameter perturbation.

[0057] The above merely describes preferred specific embodiments of the present application, but the protection scope of the present application is not limited thereto, and any modification, equivalent replacement and improvement made by any person skilled in the art within the technical range disclosed by the present application, as long as it is within the spirit and principle of the present application, should be covered within the protection scope of the present application.

Claims

1. A method for calculating wave velocity-attenuation throughout the entire process based on the Lo two-phase fluid model, characterized in that, The method includes the following steps: S1, Construction of medium properties: Based on the temperature and pressure conditions of the target reservoir, the thermodynamic parameters of CO2-brine are obtained. Based on the experimental data, the capillary pressure-saturation relationship seepage characteristics are obtained through fitting formulas. The two-phase fluid parameters and rock skeleton parameters are combined to construct a unified medium property layer. Among them, the capillary pressure-saturation relationship is characterized by an extended van Genuchten model to correct the capillary pressure fitting error in the low and medium saturation range. S2, Wave velocity and attenuation calculation: The elastic, inertial, viscous and shear modulus matrices are calculated by the parameters of the medium physical layer, and then substituted into the control equation of the Lo two-phase fluid model for solution, to obtain the phase velocity and attenuation coefficient of three types of longitudinal waves and one type of transverse wave in CO2-salt water saturated porous medium. S3, Sensitivity and Mechanism Analysis: Based on the wave velocity and attenuation calculation results, local and global sensitivity analyses are performed on various physical property parameters to identify key parameters affecting wave velocity; among them, a three-channel decomposition is performed on the P1 wave, and the influence of parameters on wave velocity is divided into three parts: inertial channel, viscous channel and elastic channel, clarifying the sensitivity composition under different physical mechanisms.

2. The wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model according to claim 1, characterized in that, In step S1, the thermodynamic parameters of the CO2-salt water include density, sound velocity, and dynamic viscosity; wherein, the density and sound velocity of CO2 are calculated using the Span-Wagner equation of state, and the expression is: In the formula, Pressure, unit is MPa; CO2 density It is the gas constant; Temperature, in Kelvin (K). To reduce density, , Critical density, ; The speed of sound of CO2 To reduce the temperature, , It is the critical temperature; For the dimensionless Helmholtz free energy of the ideal gas part, The dimensionless Helmholtz free energy of the residual part. The first-order partial derivative of the dimensionless Helmholtz free energy of the residual part with respect to the reduced density. ; The second-order mixed partial derivative of the residual dimensionless Helmholtz free energy with respect to the reduced density and reduced temperature is given. ; For the dimensionless Helmholtz free energy of the ideal gas part, the second partial derivative with respect to the reduction temperature is given. ; The second partial derivative of the dimensionless Helmholtz free energy of the residual part with respect to the reduced density is given by . ; The dynamic viscosity of CO2 was calculated using the Laesecke-Muzny model: In the formula, For CO2 viscosity, Zero density limiting viscosity, For linear density terms, This is the remaining viscosity term. This is a critical enhancing term; The density, sound velocity, and dynamic viscosity of the salt water were calculated using the Batzle-Wang empirical formula, expressed as follows: In the formula, The density of the salt water is... For the speed of sound in salt water, For the speed of sound in pure water, , These are the parameters for the experimental fitting. This represents the mass fraction of NaCl. This refers to the viscosity of the salt water.

3. The wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model according to claim 1, characterized in that, In step S1, after obtaining the thermodynamic parameters of CO2-salt water, the relative permeability model is used to describe the flow capacity of the two phases of CO2-salt water in the porous medium, and the expression is: In the formula, The relative permeability of CO2. The relative permeability of the saline solution. CO2 saturation All of these are experimental fitting parameters.

4. The wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model according to claim 3, characterized in that, In step S1, the capillary pressure-saturation relationship is characterized using an extended van Genuchten model to correct the capillary pressure fitting error in the low-to-medium saturation range, including: The correction item is The core shape function, combined with linear weighting coefficients Without altering the original physical meaning of the vG parameters, by making local fine-tuning to the curve shape in the moderate saturation range, the expression for the extended van Genuchten capillary pressure model becomes: In the formula, For capillary pressure, The density of the salt water is... It is the acceleration due to gravity. All parameters are fitted to experimental data; The skeleton parameters consist of bulk modulus, shear modulus, porosity, and permeability, and are used to characterize the basic physical properties of the reservoir's solid framework. The skeleton and fluid parameters serve as the complete physical property inputs for seismic wave propagation calculations.

5. The wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model according to claim 1, characterized in that, In step S2, after the medium property layer is constructed and all parameters of the rock skeleton and pore fluid are obtained, the Lo two-phase fluid model is used to obtain the propagation characteristics of different wave types in the CO2-salt water system. This model unifies the elastic response of the solid framework, the inertial effect and viscous dissipation of the pore fluid, and the interfacial coupling effect generated by the capillary pressure gradient in the wave equation. By solving the propagation matrix in the form of complex eigenvalues, the phase velocity, dispersion relationship and attenuation coefficient of three types of longitudinal waves P1, P2, P3 and one type of transverse wave S wave are calculated, and the response law with frequency, saturation and pore structure changes is described.

6. The wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model according to claim 5, characterized in that, The Lo two-phase fluid model, based on the Euler description, treats the solid skeleton, the non-wetting fluid phase, and the wetting fluid phase as a three-phase coexisting system. It describes the motion behavior of each phase through mass and momentum conservation equations and constructs the elastic response of the solid phase by combining linear stress-strain relationships. It introduces inertial coupling between the solid phase and the fluid phase, cross viscous drag, and capillary pressure gradient terms to unify the deformation of the solid skeleton, the inertial effect of the pore fluid, viscous dissipation, and interfacial tension processes. The constitutive equation for a three-phase medium is: In the formula, For medium displacement, , As a solid phase, It is a non-wetting flow phase. For wetting the flow phase, As a vector differential operator, it can be represented in Cartesian coordinates as: , The first time derivative of the medium displacement. The second time derivative of the medium displacement. The inertial coupling matrix characterizes the virtual mass effect caused by the density of each phase and the acceleration of different phases, including the inertial coupling between the solid phase and each fluid phase, as well as the cross-inertial coupling between fluids; This is a viscous coupling matrix that describes Darcy drag and Yuster-type cross viscous drag. The elasticity coefficient matrix, This is the shear modulus coefficient matrix; The expressions for the inertial coupling matrix, viscous coupling matrix, and elastic coefficient matrix are as follows: In the formula, for Phase density, ; for Phase volume fraction, ; These are the constitutive coefficients related to the inertial coupling between the solid and fluid phases, respectively. Let be the inertial cross-coupling coefficients between the two fluids, respectively. , are constitutive coefficients related to the viscous coupling between the solid and the fluid, respectively. Let be the cross-coupling coefficients between the two fluids, respectively. The elasticity coefficient is symmetric with cross terms. , ; The shear modulus of the porous skeleton; The elastic coefficient is expressed as: In the formula, for The bulk modulus of the phase, ; For rock skeleton porosity, Dimensionless parameters describing porosity closure were obtained through both sleeveless and sleeved experiments. ; These are the first and second water storage coefficients, respectively. For capillary pressure, CO2 saturation; Then we have: In the formula, The bulk modulus of the rock skeleton; Obtained from the extended vG model: In the formula, The rate of change of capillary pressure with respect to saturation; The expression for the inertia coefficient is: In the formula, The pore tortuosity factor has a theoretical value of 3 for a random system with uniform circular pores of all orientations; a value of 1 when the pores are uniform and parallel to the axis and pressure gradient; and a value of 1 when the solid particles in the porous medium are spherical. ; The viscosity parameter is expressed as: Represented as: In the formula, for The dynamic viscosity of the phase, ; for The relative permeability of the phase, ; This is the inherent penetration rate; By using the constitutive equations of a three-phase medium to obtain the divergence and curl, the governing equations for longitudinal and transverse waves are derived: Introducing the general form of the solution to the constitutive equation for a three-phase medium: In the formula, for longitudinal displacement of the phase, for Longitudinal wave amplitude, for Lateral displacement of the phase, for Phase transverse wave amplitude, For complex wave number, For the direction of dissemination, Angular frequency, It is a time variable; Substituting the general form of the solution to the constitutive equation of the three-phase medium into the control equations of the longitudinal and transverse waves, we obtain the final control equations of the longitudinal and transverse waves. If the determinant of the coefficient matrix is ​​set to 0, then the governing equation has a trivial solution: After unfolding, it becomes The polynomials of the equations have six roots for the longitudinal wave control equations and two roots for the transverse wave control equations. Under the physical constraint that the wave always decreases along the direction of propagation, there are three types of longitudinal waves and one type of transverse wave.

7. The wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model according to claim 1, characterized in that, In step S3, based on the wave velocity and attenuation calculation results, local and global sensitivity analyses are performed on various physical property parameters to identify key parameters affecting wave velocity, including: Sensitivity analysis was performed at the global level using the PAWN method; for each output parameter... Divide the entire interval into There are 10 sub-intervals, each sub-interval has a fixed number of sub-intervals. The remaining parameters are sampled in the entire space using the Latin hypercube sampling method to obtain the output. Conditional distribution The overall influence of each input parameter is quantified by calculating the Kolmogorov-Smirnov statistic between the conditional and unconditional distributions of the output. ,but: In the formula, In the first Given a given interval, the maximum absolute difference between the conditional cumulative distribution function and the unconditional cumulative distribution function. To take the upper bound for all possible inputs, To divide according to the range of values ​​of the input parameters After the first sub-interval, the second... A range; In order to meet the conditions Falling in each interval The conditional cumulative distribution function under the given conditions, It is the unconditional accumulation distribution function; For each parameter By combining the KS statistics of all sub-intervals, a global sensitivity index of PAWN is defined. : In the formula, To all By performing summary statistics, the final PAWN sensitivity index is obtained. ; Local sensitivity analysis is performed on the wave velocity and attenuation results to characterize the instantaneous response of the model to physical property disturbances under typical operating conditions. A disturbance is applied to each input parameter near a given set of reference parameters, and the difference between the output before and after the disturbance is calculated to obtain the first-order sensitivity of the output with respect to the parameter. For output variables With input parameters The local sensitivity coefficient is expressed as: In the formula, For local sensitivity parameters, The ring parameters are used to determine the local sensitivity parameters.

8. The wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model according to claim 6, characterized in that, In step S3, the P1 wave is decomposed into three channels, including: A three-channel decomposition method was used to break down the changes in the corresponding eigenvalues ​​of the P1 wave into contributions from three physical channels: elastic-capillary, inertial, and viscous, in order to identify the dominant path controlling the change in wave velocity. The critical saturation degree at which the P1 wave changes from decreasing to increasing is determined by the velocity inflection point determination method, and the dependence of this inflection point on various physical property parameters is analyzed. Under the frequency domain eigenvalue problem of the Lo model, the velocity of the P1 wave is... about The first-order sensitivity can be written as the sum of the inertial, elastic-capillary, and viscous channels; the characteristic equation can be rewritten as a generalized eigenvalue problem, expressed as: Introducing left eigenvalues and apply normalization conditions. The eigenvalues ​​are written as: Under the above normalization conditions, the derivative of the mechanism variable is taken with respect to porosity. For example, left multiplication The first derivatives of the eigenvalues ​​are obtained: For body waves, at a fixed frequency, we have: Obtain the velocity with respect to porosity The first-order sensitivity is decomposed into the sum of the elastic, inertial, and viscous channels: in: In the formula, It is an elastic capillary channel. , For inertial channels, It is a viscous channel.

9. The wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model according to claim 8, characterized in that, After quantifying the contributions of the elastic-capillary, inertial, and viscous channels to the P1 wave velocity sensitivity using a three-channel decomposition method, the inflection point of the P1 wave velocity was determined, thereby identifying the critical saturation and controlling factors at which the velocity changes from decreasing to increasing. The approximate expression for the P1 wave velocity is as follows: In the formula, It is the equivalent bulk modulus of the medium, which is determined by the solid skeleton, fluid compressibility, and coupling effects. It is the equivalent density of the medium, obtained by the linear superposition of the contributions of each phase; the relationship between the equivalent bulk modulus and the equivalent density is analyzed. Relative sensitivity: In the formula, It is the equivalent bulk modulus pair Sensitivity It is an equivalent density pair The sensitivity of equivalent bulk modulus and equivalent density to saturation; when the relative rates of change of equivalent bulk modulus and equivalent density are equal, the sensitivity of equivalent bulk modulus to equivalent density is... The relative sensitivity expression gives a sensitivity term of zero, and the P1 wave velocity reaches a local minimum. This condition is the critical saturation criterion for the velocity to change from decreasing to increasing.

10. The application of the wave velocity-attenuation full-process calculation method based on the Lo two-phase fluid model as described in any one of claims 1-9 in the seismic wave response simulation, monitoring data interpretation, or reservoir parameter inversion of CO2 transport during the geological carbon sequestration process in deep saline aquifers.

Images