Rock physical forward modeling method suitable for formation pore fluid of carbon dioxide

Abstract

The invention discloses a rock physical forward modeling method applicable to formation pore fluid of carbon dioxide, and relates to the technical field of petroleum and natural gas engineering. The method comprises the steps of firstly obtaining a carbon dioxide physical property constraint condition based on a state equation of carbon dioxide under a reservoir temperature and pressure condition, and then based on the constraint condition and according to logging information of a target well, core analysis and a laboratory test result, carrying out formation evaluation to calculate reservoir parameters. Then, based on reservoir parameters, the sandstone content is obtained through calculation, a rock physical volume model of a target well is constructed, and finally, based on the model, corresponding rock physical forward modeling is conducted according to different carbon dioxide saturations by means of the equivalent medium theory to generate stratum elastic attribute parameters. And the change relation between the corresponding saturation and the stratum elastic attribute parameters is quantitatively represented, so that the CO2 quantitative monitoring precision and reliability can be remarkably improved by comprehensively considering the special physical properties of CO2, the fluid replacement mechanism and the seismic wave propagation characteristics.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of petroleum and natural gas engineering, and particularly relates to a rock physics forward method suitable for a formation pore fluid being carbon dioxide. BACKGROUND

[0002] Under the global energy conservation and environmental protection background, carbon capture, utilization and storage (CCUS) as a key technology for reducing greenhouse gas emissions aims to capture carbon dioxide CO2 from the source and store it safely and effectively. A large number of previous studies have shown that CO2 geological storage is an important means to achieve the goal of energy conservation and environmental protection.

[0003] The safety and effectiveness of CO2 storage depend on the accurate monitoring of CO2 migration and distribution in the reservoir. Geophysical monitoring methods are considered to be the most effective means of CO2 monitoring, which can qualitatively describe the dynamic changes of the reservoir. However, the traditional seismic monitoring method (such as time-lapse seismic) still has certain limitations in the quantitative characterization of key parameters such as CO2 saturation. Meanwhile, in the field of CO2 geological storage, many works mainly focus on the optimization of seismic forward method, and few works focus on quantitative characterization research. In addition, the existing CO2 monitoring technology lacks systematic research on fluid displacement effect and seismic wave field response, which limits the accuracy and reliability of CO2 quantitative monitoring, and requires accurate rock physics forward method for quantitative support.

[0004] However, the existing rock physics forward method is mainly established for oil and gas exploration and development environment, and usually assumes that the formation pore is mainly filled with water, oil or gas, and the rock physics model and fluid replacement method used in the calculation process do not fully consider the unique physical characteristics of CO2. For example, CO2 may be in a supercritical state under formation conditions, and its density and elastic properties are highly sensitive to changes in temperature and pressure; at the same time, there is a complex interface effect between CO2 and formation water and rock matrix, which will significantly affect the wave impedance, velocity and reflection characteristics. If the traditional rock physics forward model in oil and gas exploration is followed, it is easy to cause distortion in the description of the seismic response of CO2 injection layer and surrounding rock layer.

[0005] Therefore, there is an urgent need for a rock physics forward method that can adapt to the characteristics of CO2 storage, so as to comprehensively consider the special physical properties of CO2, the fluid replacement mechanism and the seismic wave propagation characteristics, and thus provide more reliable theoretical and technical support for the seismic monitoring and reservoir evaluation of CO2 geological storage, and provide a more reliable quantitative interpretation framework for CO2 storage monitoring, and also has important theoretical significance and engineering value for promoting the industrialization application of CCUS technology. SUMMARY

[0006] The purpose of the present application is to provide a rock physics forward method suitable for a formation pore fluid being carbon dioxide, so as to solve the problem that the existing rock physics forward scheme is distorted in describing the seismic response of the CO2 injection layer and the surrounding rock layer, thereby affecting the accurate monitoring of the CO2 migration and distribution in the reservoir, so as to provide more reliable theoretical and technical support for the seismic monitoring and reservoir evaluation of CO2 geological storage, and to provide a more reliable quantitative interpretation framework for CO2 storage monitoring, and also has important theoretical significance and engineering value for promoting the industrial application of CCUS technology.

[0007] In order to achieve the above-mentioned purpose, the present application adopts the following technical solutions:

[0008] The present application provides a rock physics forward method suitable for a formation pore fluid being carbon dioxide, comprising:

[0009] The state equation of carbon dioxide under the temperature and pressure conditions of the reservoir is introduced, and the dynamic response relationship of the physical and chemical properties of pure carbon dioxide in the gaseous, liquid and supercritical states with temperature and pressure is obtained based on the state equation, so as to serve as a carbon dioxide physical property constraint condition;

[0010] Based on the carbon dioxide physical property constraint condition, and according to the logging data, core analysis and laboratory test results of the target well, the formation evaluation is carried out to calculate the reservoir parameters, wherein the reservoir parameters include shale content and comprehensive porosity;

[0011] The sandstone content is calculated based on the reservoir parameters, and the rock physics volume model of the target well is constructed using the sandstone content and the reservoir parameters;

[0012] Based on the rock physics volume model, for different carbon dioxide saturations, the corresponding rock physics forward is carried out using the equivalent medium theory to generate the formation elastic property parameters, and the change relationship between the corresponding saturation and the formation elastic property parameters is quantitatively characterized.

[0013] Based on the above invention content, a new scheme capable of adapting to the CO2 storage characteristics for rock physics forward is provided, that is, first, the CO2 physical property constraint condition is obtained based on the equation of state of CO2 under the temperature and pressure conditions of the reservoir, then the constraint condition is used, and the formation evaluation is carried out according to the well logging data, core analysis and laboratory test results of the target well to calculate the reservoir parameters, then the sandstone content is calculated based on the reservoir parameters and the rock physics volume model of the target well is constructed, finally, based on the model, for different CO2 saturations, the corresponding rock physics forward is carried out by using the equivalent medium theory to generate the formation elastic property parameters, and the change relationship between the corresponding saturation and the formation elastic property parameters is quantitatively characterized, so that by comprehensively considering the special physical properties of CO2, the fluid replacement mechanism and the seismic wave propagation characteristics in the rock physics forward process, the accuracy and reliability of the CO2 quantitative monitoring can be significantly improved, and then the scheme can be used as a key technical means for carbon capture, utilization and storage, and has a unique and irreplaceable role in the CO2 geological storage project, which is convenient for practical application and popularization.

[0014] In one possible design, the critical point of the supercritical state is a pressure of 7.38 MPa and a temperature of 31.1℃;

[0015] And / or, the equation of state adopts a Span-Wagner equation, wherein the Span-Wagner equation is a multi-parameter equation of state based on Helmholtz free energy, and is expressed as follows:

[0016]

[0017] In the formula, A(ρ, T) represents the Helmholtz free energy, represents the free energy part when the fluid is assumed to be an ideal gas, represents the difference free energy part between the real fluid and the ideal gas, T represents the absolute temperature, R represents the gas constant, δ represents the reduced density and has δ = ρ ÷ ρ c , ρ represents the density, ρ c = 467.6 kg·m -3 , τ represents the reduced inverse temperature and has τ = T c ÷ T, T c = 304.1282 K, ln(δ) represents the compressibility factor relationship of the ideal gas, represents the temperature dependence of the constant-pressure heat capacity of the ideal gas, represents a constant term, represents a linear temperature term, represents a logarithmic temperature term, represents the corresponding vibration temperature and is used to control the starting point of the change of the term with temperature, n irepresents the coefficient determining the weight of each term, a i represents the first power exponent and is used to control the dependence of the non-analytical term on the density, b i represents the second power exponent and is used to control the degree of singularity of the term on the critical point distance function, d i represents the third power exponent and is used to control the sensitivity of the term to the density, t i represents the fourth power exponent and is used to control the sensitivity of the term to the temperature, c i represents the decay factor, a i controls the width of the Gaussian function in the δ direction, b i controls the width of the Gaussian function in the τ direction, e i represents the position of the Gaussian center point in the δ direction, g i represents the position of the Gaussian center point in the τ direction, e represents the base of the natural logarithm, and Δ represents the critical zone distance function.

[0018] In one possible design, the calculation of the shale content includes the following:

[0019] The natural gamma logging method is used to determine the shale content of the rock, wherein the shale content is calculated by using the following linear interpolation formula:

[0020]

[0021] wherein V sh represents the shale content, GR log represents the gamma value of the measuring point, GR min represents the gamma value of pure sand, GR max represents the gamma value of pure mud.

[0022] In one possible design, the calculation of the comprehensive porosity includes the following:

[0023] The first porosity reflecting the difference between the total density of the rock and the matrix and the density of the pore fluid is calculated based on the density logging curve as the density porosity, and the density porosity is corrected by using the following formula:

[0024]

[0025] In the formula, PHIDC represents the corrected density porosity, PHID represents the density porosity under the corresponding matrix, V 2.76 represents the density corresponding to dry clay, p 2.65 represents the density of sandstone, and p1 represents the density of the formation pore fluid.

[0026] A second porosity for reflecting the influence of pores on the velocity of acoustic wave propagation is calculated based on acoustic logging curves as acoustic porosity, and the acoustic porosity is corrected by using the following formula:

[0027]

[0028] In the formula, PHIDTC represents the corrected acoustic porosity, PHIDT represents the acoustic porosity under the corresponding matrix, DT_VCL represents the acoustic time difference corresponding to dry clay, γ 56 represents the P-wave time difference of sandstone, γ 189 represents the P-wave time difference of formation pore fluid;

[0029] A third porosity for reflecting the content of hydrogen atoms in rock is calculated based on neutron logging curves as neutron porosity, and the neutron porosity is corrected by using the following formula:

[0030]

[0031] In the formula, PHINC represents the corrected neutron porosity, PHIN represents the neutron porosity under the corresponding matrix, δ -0.02 represents the neutron porosity of sandstone, δ 0.36 represents the neutron porosity corresponding to dry clay, and δ1 represents the hydrogen index of water;

[0032] The neutron-density porosity is corrected by using the following formula:

[0033] PHINDC=PHIND-VCLAY×δ 0.18

[0034] In the formula, PHINDC represents the corrected neutron-density porosity, PHIND represents the neutron-density porosity under the corresponding matrix and has PHIND=PHIN÷PHID, and δ 0.18 represents the neutron porosity of typical shale;

[0035] According to the density porosity, the acoustic porosity, the neutron porosity and the neutron-density porosity, the comprehensive porosity PHIT is calculated by using the following formula:

[0036] PHIT=PHIDC×ω1+PHINDC×ω2+PHINC×ω3+PHIDTC×ω4

[0037] In the formula, ω1, ω2, ω3 and ω4 respectively represent preset weight coefficients and have ω1+ω2+ω3+ω4=1.

[0038] In one possible design, the sandstone content is calculated based on the reservoir parameters, including:

[0039] Based on the clay content VCLAY and the comprehensive porosity PHIT in the reservoir parameters, the sandstone content VQUA = 1 - PHIT - VCLAY is calculated.

[0040] In one possible design, the carbon dioxide saturation ranges from [0% to 100%], with the different carbon dioxide saturations set at intervals of 0%, 20%, 40%, 60%, 80%, and 100%.

[0041] In one possible design, the formation elastic property parameters include density, P-wave transit time, S-wave transit time, and / or P-wave / S-wave velocity ratio.

[0042] In one possible design, the equivalent medium theory posits that the macroscopic properties of the rock are jointly determined by the properties of a skeleton model, a dry rock model with added pores, and a saturated rock model with added fluids. The skeleton model is constructed based on the arithmetic mean of the Voigt upper bound and the Reuss lower bound of the Voigt-Reuss-Hill model, and is expressed as follows:

[0043]

[0044] In the formula, M VRH M represents the average modulus of Hill. V Voigt represents the upper limit of the equivalent elastic modulus and has N represents the total number of components, i represents the component number, and f i M represents the volume component of the i-th component. i M represents the elastic modulus of the i-th component. R This represents the Reuss lower limit of the equivalent elastic modulus and has

[0045] In one possible design, the dry rock model is constructed based on the differential equivalent medium theory, which simulates a two-phase mixture by gradually adding inclusion phases into the solid mineral phase. Specifically, it employs an iterative calculation algorithm to incrementally embed pore structures of different geometries into the mineral framework phase, with an equivalent bulk modulus K. * and shear modulus μ * The coupled differential equations are expressed as follows:

[0046]

[0047] In the formula, the equivalent bulk modulus K * initial value K * (0) = K1, the shear modulus μ * initial value μ* (0) = μ1, K1 represents the bulk modulus of the initial principal material, μ1 represents the shear modulus of the initial principal material, K2 represents the bulk modulus of the gradually added inclusions, the gradually added inclusions being pores of different geometries, μ2 represents the shear modulus of the gradually added inclusions, y represents the porosity of the gradually added inclusions, K * (y) represents the equivalent bulk modulus K when the porosity is y. * μ * (y) represents the shear modulus μ when the porosity is y. * P (*2) (y) represents the porosity of y and for the condition having an equivalent bulk modulus K. * and shear modulus μ * The first geometric factor, Q, of the inclusion material in the background medium. (*2) (y) represents the porosity of y and for the condition having an equivalent bulk modulus K. * and shear modulus μ * The second geometric factor of the containing material in the background medium.

[0048] In one possible design, the saturated rock model is constructed based on the Gassmann fluid substitution equation, which describes the modulus change from a dry rock state to a saturated fluid porosity state, and is expressed as follows:

[0049]

[0050] In the formula, K sat1 K represents the effective bulk modulus of saturated rock for the first saturated fluid. sat2 K represents the effective bulk modulus of saturated rock for the second saturated fluid, K0 represents the bulk modulus of the minerals composing the rock, and K fl1 K represents the effective bulk modulus of the fluid in the first pore. fl2 This represents the effective bulk modulus of the fluid in the second pore. Indicates porosity.

[0051] The beneficial effects of the above scheme are:

[0052] (1) The application provides a new scheme for rock physics forward modeling capable of adapting to CO2 storage characteristics, that is, first, CO2 physical property constraint conditions are obtained based on a state equation of CO2 under temperature and pressure conditions of a reservoir, then, based on the constraint conditions and according to well logging data, core analysis and laboratory test results of a target well, formation evaluation is carried out to calculate reservoir parameters, then, sandstone content is calculated based on the reservoir parameters and a rock physics volume model of the target well is constructed, and finally, based on the model, corresponding rock physics forward modeling is carried out for different CO2 saturations by using equivalent medium theory to generate formation elastic property parameters and quantitatively represent the change relationship between the corresponding saturations and the formation elastic property parameters, so that the accuracy and reliability of CO2 quantitative monitoring can be significantly improved by comprehensively considering the special physical properties of CO2, fluid replacement mechanism and seismic wave propagation characteristics in the process of rock physics forward modeling, and the application can be used as a key technical means for carbon capture, utilization and storage and has a unique and irreplaceable role in CO2 geological storage engineering.

[0053] (2) Compared with the traditional method, the application can accurately reflect the sensitivity of different saturation CO2 to the change of formation elastic property parameters, and improves the limitations of geophysical monitoring methods in the quantitative representation of key parameters such as CO2 saturation.

[0054] (3) Compared with the traditional method, the application improves the accuracy and reliability of CO2 quantitative monitoring based on the quantitative support of CO2 rock physics forward modeling, and has good engineering adaptability, can be embedded into mainstream seismic forward and inversion processes, and is suitable for multiple links such as reservoir evaluation, injection monitoring and long-term safety evaluation, thereby providing reliable technical support for the whole life cycle of CO2 geological storage engineering, and has the advantages of high precision, strong applicability and significant engineering popularization value.

[0055] (4) The application can provide more reliable theoretical and technical support for seismic monitoring and reservoir evaluation of CO2 geological storage, and provide a more reliable quantitative interpretation framework for CO2 storage monitoring, and has important theoretical significance and engineering value for promoting the industrial application of CCUS technology, and is convenient for practical application and popularization. BRIEF DESCRIPTION OF DRAWINGS

[0056] In order to more clearly illustrate the technical solutions in the embodiments of the application or the prior art, the drawings needed to be used in the embodiments or prior art description will be briefly introduced. Obviously, the drawings in the following description only some embodiments of the application, and for those skilled in the art, other drawings can be obtained without creative labor on the basis of these drawings.

[0057] Figure 1A flowchart of a rock physics forward method suitable for a formation pore fluid of carbon dioxide is provided in the embodiments of the present application. DETAILED DESCRIPTION

[0058] In order to more clearly illustrate the technical solutions of the embodiments of the present application or the prior art, the present application will be briefly introduced below in combination with the drawings and the description of the embodiments or the prior art. Obviously, the following description of the drawings is only some embodiments of the present application, and for those skilled in the art, other embodiments can be obtained without creative labor on the basis of these embodiments. It should be noted that the description of these embodiments is used to help understand the present application, but does not constitute a limitation on the present application.

[0059] It should be understood that although the terms first and second, etc. may be used herein to describe various objects, these objects should not be limited by these terms. These terms are only used to distinguish one object from another object. For example, the first object can be called the second object, and similarly, the second object can be called the first object, without departing from the scope of the example embodiments of the present application.

[0060] It should be understood that for the term "and / or" which may appear in the present application, it is only a description of the association relationship of the associated objects, which means that there can be three relationships, for example, A and / or B can mean that A exists alone, B exists alone, or A and B exist together, etc. For example, A, B and / or C can mean that any one of A, B and C exists or any combination thereof; for the term " / and" which may appear in the present application, it is another description of the association relationship of another associated object, which means that there can be two relationships, for example, A / and B can mean that A exists alone or A and B exist together; in addition, for the character " / " which may appear in the present application, it generally means that the associated objects before and after are in an "or" relationship.

[0061] EMBODIMENT

[0062] As shown in the embodiment, the rock physics forward method suitable for a formation pore fluid of carbon dioxide comprises but is not limited to the following steps S1-S4. Figure 1

[0063] S1. Introducing a state equation of carbon dioxide under the temperature and pressure conditions of a reservoir, and obtaining a dynamic response relationship of physical and chemical properties of pure carbon dioxide in gaseous, liquid and supercritical states with temperature and pressure changes based on the state equation as a carbon dioxide physical property constraint condition.

[0064] ​In the step S1, the equation of state can specifically but not limited to adopt the Span-Wagner equation, wherein the Span-Wagner equation is a multi-parameter equation of state based on Helmholtz free energy, and is expressed as follows:

[0065]

[0066] In the formula, A(ρ, T) represents the Helmholtz free energy, represents the free energy part when the fluid is assumed to be an ideal gas, represents the difference free energy part between the real fluid and the ideal gas, T represents the absolute temperature, R represents the gas constant, δ represents the reduced density and has δ = ρ ÷ ρ c , ρ represents the density, ρ c = 467.6 kg·m -3 , τ represents the reduced inverse temperature and has τ = T c ÷ T, T c = 304.1282 K, ln(δ) represents the compressibility factor relationship of the ideal gas, represents the temperature dependence corresponding to the constant volume heat capacity of the ideal gas, represents the constant term, represents the linear temperature term, represents the logarithmic temperature term, represents the vibration temperature corresponding to and is used to control the starting point of the term with the change of temperature, n i represents the coefficient to determine the weight of each term, a i represents the first power exponent and is used to control the dependence of the non-analytic term in the density direction, b i represents the second power exponent and is used to control the degree of singularity of the term on the critical point distance function, d i represents the third power exponent and is used to control the sensitivity of the term to the density, t i represents the fourth power exponent and is used to control the sensitivity of the term to the temperature, c i represents the attenuation factor, α i controls the width of the Gaussian function in the δ direction, β i controls the width of the Gaussian function in the τ direction, ∈ i represents the position of the Gaussian center point in the δ direction, γ i represents the position of the Gaussian center point in the τ direction, e represents the base of the natural logarithm, and Δ represents the critical region distance function. The gas constant R can generally take 0.1889241 kJ·kg -1 · K -1, the critical region distance function Δ is used to ensure the correct critical singularity, the Helmholtz free energy (i.e. Helmholtz free energy) is a state function in thermodynamics, defined as the internal energy of the system minus the product of absolute temperature and entropy, i.e. A = U - TS (where U is the internal energy, T is the absolute temperature, and S is the entropy), which is usually decomposed into the ideal gas part of Helmholtz energy (i.e. ) and the residual part of Helmholtz energy (i.e. ), which is an energy index describing the ability of the system to do work under constant temperature conditions, commonly used to judge the thermodynamic equilibrium state of the isothermal process. In addition, the critical point of the supercritical state is a pressure of 7.38 MPa and a temperature of 31.1℃, and the supercritical state is usually met when the formation depth exceeds 800-1000 meters; and the physical and chemical properties of pure carbon dioxide include but are not limited to density and elasticity.

[0067] S2. Based on the carbon dioxide property constraint condition, and according to the logging data, core analysis and laboratory test results of the target well, carry out formation evaluation to calculate the reservoir parameters, wherein the reservoir parameters include but are not limited to shale content and comprehensive porosity.

[0068] In the step S2, there are many methods to calculate the shale content, such as methods based on natural gamma, natural potential, resistivity, neutron and neutron-density crossplot, etc. These methods have their own advantages and disadvantages. In order to ensure the accuracy of the calculation results, the natural gamma logging method is preferred to determine the shale content of the rock, wherein the technical principle of the natural gamma logging method is: first, by measuring the natural gamma intensity of the formation and comparing with the pure sand and pure mud as the benchmark, the content of radioactive elements such as potassium, uranium and / or thorium in shale is used to quantitatively estimate the shale content, i.e. the natural gamma logging method can be used to calculate the shale content by using the following linear interpolation formula:

[0069]

[0070] Wherein, V sh represents the shale content, GR log represents the gamma value of the measuring point, GR min represents the gamma value of pure sand, and GR max represents the gamma value of pure mud.

[0071] In the step S2, the comprehensive porosity refers to the weighted average calculation result of the density porosity, acoustic porosity, neutron porosity and neutron-density porosity. Specifically, the calculation method of the comprehensive porosity includes but is not limited to the following steps S21-S25.

[0072] S21. Based on the density logging curve, a first porosity is calculated to reflect the difference between the overall rock density and the matrix and pore fluid density, which is then used as the density porosity. The density porosity is corrected using the following formula:

[0073]

[0074] In the formula, PHIDC represents the corrected density porosity, PHID represents the density porosity under the corresponding skeleton, VCLAY represents the clay content, and ρ 2.76 ρ represents the density corresponding to dry clay. 2.65 ρ represents the density of sandstone, and ρ1 represents the density of pore fluid in the formation.

[0075] In step S21, the density logging curve is one of the results obtained from formation evaluation based on logging data, core analysis, and laboratory test results. The first porosity, which reflects the difference between the overall rock density and the matrix and pore fluid density, can be conventionally calculated. Specifically, in the above formula, density ρ... 2.76 For example, the dimensionless value of ρ can be taken as 2.76, and the density ρ 2.65 For example, the dimensionless value of ρ can be taken as 2.65, and the dimensionless value of density ρ1 can be taken as 1.

[0076] S22. A second porosity reflecting the influence of pores on the propagation velocity of sound waves is calculated based on the sonic logging curve and used as the sonic porosity, and the sonic porosity is corrected using the formula:

[0077]

[0078] In the formula, PHIDTC represents the corrected acoustic porosity, PHIDT represents the acoustic porosity under the corresponding skeleton, DT_VCL represents the acoustic transit time corresponding to dry clay, and γ 56 The longitudinal wave transit time of sandstone, γ 189 This represents the P-wave travel time of the formation pore fluid.

[0079] In step S22, the sonic logging curve is the second result obtained from formation evaluation based on logging data, core analysis, and laboratory test results. The second porosity, reflecting the influence of pores on the sonic wave propagation velocity, can be conventionally calculated. Specifically, in the above formula, the P-wave travel time γ... 56 For example, the dimensionless value of γ can be taken as 56, representing the longitudinal wave time difference. 189 The dimensionless value can be taken as 189 for example. Furthermore, the sonic logging curves are primarily P-wave time-of-flight curves.

[0080] S23. Calculate a third porosity for reflecting the content of hydrogen atoms in the rock based on the neutron logging curve as a neutron porosity, and correct the neutron porosity by using the following formula:

[0081]

[0082] In the formula, PHINC represents the corrected neutron porosity, PHIN represents the neutron porosity under the corresponding matrix, and δ -0.02 represents the neutron porosity of sandstone, δ 0.36 represents the neutron porosity corresponding to dry clay, and δ1 represents the hydrogen index of water.

[0083] In the step S23, the neutron logging curve is one of the results obtained from the formation evaluation based on logging data, core analysis and laboratory test results, and the third porosity for reflecting the content of hydrogen atoms in the rock can be calculated conventionally. In detail, in the above formula, the dimensionless value of the neutron porosity δ -0.02 may be taken as -0.02, the dimensionless value of the neutron porosity δ 0.36 may be taken as 0.36, and the dimensionless value of the hydrogen index δ1 may be taken as 1.

[0084] S24. Correct to obtain a neutron-density porosity according to the following formula:

[0085] PHINDC = PHIND - VCLAY x δ 0.18

[0086] In the formula, PHINDC represents the corrected neutron-density porosity, PHIND represents the neutron-density porosity under the corresponding matrix and has PHIND = PHIN ÷ PHID, and δ 0.18 represents the neutron porosity of typical shale.

[0087] In the step S24, in detail, the dimensionless value of the neutron porosity δ 0.18 may be taken as 0.18.

[0088] S25. Calculate a comprehensive porosity PHIT according to the following formula based on the density porosity, the acoustic porosity, the neutron porosity and the neutron-density porosity:

[0089] PHIT = PHIDC x ω1 + PHINDC x ω2 + PHINC x ω3 + PHIDTC x ω4

[0090] In the formula, ω1, ω2, ω3 and ω4 respectively represent preset weight coefficients and have ω1 + ω2 + ω3 + ω4 = 1.

[0091] In the step S25, in detail, the weight coefficients ω1, ω2, ω3 and ω4 can be exemplified as ω1=0.35, ω2=0.35, ω3=0.2 and ω4=0.1 respectively.

[0092] S3. Calculate the sand content based on the reservoir parameters, and construct a rock physics volume model of the target well by using the sand content and the reservoir parameters.

[0093] In the step S3, in detail, the sand content calculated based on the reservoir parameters includes but is not limited to: based on the shale content VCLAY and the comprehensive porosity PHIT in the reservoir parameters, the sand content VQUA=1-PHIT-VCLAY is calculated.

[0094] S4. Based on the rock physics volume model, for different respective carbon dioxide saturations, corresponding rock physics forward is performed by using the equivalent medium theory to generate formation elastic attribute parameters, and a change relationship between the corresponding saturation and the formation elastic attribute parameters is quantitatively characterized.

[0095] In the step S4, in detail, the value range of the carbon dioxide saturation is [0%, 100%], and the different respective carbon dioxide saturations can be exemplified but are not limited to 0%, 20%, 40%, 60%, 80% and 100% and the like which are respectively and interval set. The formation elastic attribute parameters include but are not limited to density, P-wave time difference, S-wave time difference and / or P-S wave velocity ratio and the like, which are needed to be constructed based on the equivalent medium theory. In detail, the equivalent medium theory is a method of regarding a multi-phase heterogeneous medium as a homogeneous equivalent medium, and the overall physical property parameters of the equivalent medium are calculated through a certain averaging or approximation process; the equivalent medium theory considers that the macroscopic properties of rock are jointly determined by the properties of a skeleton model, a dry rock model with added pores and a saturated rock model with added fluid. The skeleton model is constructed by the arithmetic average of the Voigt upper limit and the Reuss lower limit based on the Voigt-Reuss-Hill model (which is a theoretical model for estimating the elastic modulus of a composite material, and the effective elastic modulus is obtained by calculating the arithmetic average of the upper limit Voigt model and the lower limit Reuss model, wherein Voigt, Reuss and Hill are existing terms, the Voigt model is assumed to be independent of strain, and the upper limit is calculated by superimposing the stress of each phase, while the Reuss model is assumed to be independent of stress, and the lower limit is calculated by superimposing the strain of each phase), and is expressed as follows:

[0096]

[0097] In the formula, M VRH represents the Hill average modulus, M VVoigt upper bound of the equivalent elastic modulus and has N represents the total number of components (i.e. constituent components), i represents the component number, f i represents the volume fraction of the i-th component, M i represents the elastic modulus of the i-th component, M R Reuss lower bound of the equivalent elastic modulus and has

[0098] In the step S4, in detail, the dry rock model is constructed based on the differential effective medium theory, wherein the differential effective medium theory is to simulate the two-phase mixture by gradually adding the inclusions phase into the solid mineral phase, i.e. using the iterative calculation algorithm, the different geometric pore structures are gradually embedded into the mineral skeleton phase in an incremental manner, and the coupled differential equations of the equivalent bulk modulus K * and shear modulus μ * are respectively expressed as follows:

[0099]

[0100] In the formula, the initial value K * of the equivalent bulk modulus K * (0) = K1, the initial value μ * of the shear modulus μ * (0) = μ1, K1 represents the bulk modulus of the initial host material, μ1 represents the shear modulus of the initial host material, K2 represents the bulk modulus of the gradually added inclusions, the gradually added inclusions are different geometric pores, μ2 represents the shear modulus of the gradually added inclusions, y represents the porosity of the gradually added inclusions, K * (y) represents the equivalent bulk modulus K * at the porosity y, μ * (y) represents the shear modulus μ * at the porosity y, P (*2) (y) represents the first geometric factor at the porosity y and for the inclusion material in the background medium with the equivalent bulk modulus K * and the shear modulus μ * , Q (*2) (y) represents the second geometric factor at the porosity y and for the inclusion material in the background medium with the equivalent bulk modulus K * and the shear modulus μ * .

[0101] In the step S4, in detail, the saturated rock model is constructed based on a Gassmann fluid substitution equation used for describing modulus change from a dry rock state to a saturated fluid pore state and expressed as follows:

[0102]

[0103] wherein K sat1 represents an effective bulk modulus of the first saturated fluid, K sat2 represents an effective bulk modulus of the second saturated fluid, K0 represents a bulk modulus of a mineral constituting the rock, K fl1 represents an effective bulk modulus of the first pore fluid, K fl2 represents an effective bulk modulus of the second pore fluid, represents a porosity.

[0104] In summary, by using the method provided in the embodiment, the following technical effects are achieved:

[0105] (1) The embodiment provides a new scheme capable of adapting to CO2 storage features for rock physics forward modeling, i.e., first obtaining CO2 physical property constraint conditions based on a state equation of CO2 under temperature and pressure conditions of a reservoir, then based on the constraint conditions and according to well logging data of a target well, core analysis and laboratory test results, carrying out formation evaluation to calculate reservoir parameters, then based on the reservoir parameters, calculating sandstone content and constructing a rock physics volume model of the target well, and finally based on the model, for different CO2 saturations, using equivalent medium theory to carry out corresponding rock physics forward modeling to generate formation elastic property parameters and quantitatively characterize a change relationship between the corresponding saturation and the formation elastic property parameters, thereby by comprehensively considering special physical properties of CO2, fluid replacement mechanism and seismic wave propagation characteristics in the rock physics forward modeling process, the accuracy and reliability of CO2 quantitative monitoring can be significantly improved, and the method can be used as a key technical means for carbon capture, utilization and storage, and has a unique and irreplaceable role in CO2 geological storage engineering;

[0106] (2) Compared with a traditional method, the scheme can accurately reflect sensitivity of different saturation CO2 to change of formation elastic property parameters, and improves limitations of a geophysical monitoring method in quantitative characterization of key parameters such as CO2 saturation;

[0107] (3) Compared with the traditional method, the scheme is quantitatively supported by CO2 rock physics forward, improves the accuracy and reliability of CO2 quantitative monitoring, has good engineering adaptability, can be embedded in the mainstream seismic forward and inversion process, is suitable for multiple links such as reservoir evaluation, injection monitoring and long-term safety evaluation, thereby providing reliable technical support for the whole life cycle of CO2 geological storage engineering, and has the advantages of high precision, strong applicability and significant engineering popularization value, and beneficial effects;

[0108] (4) The scheme can provide more reliable theoretical and technical support for seismic monitoring and reservoir evaluation of CO2 geological storage, provide a more reliable quantitative interpretation framework for CO2 storage monitoring, and has important theoretical significance and engineering value for promoting the industrial application of CCUS technology, and is convenient for practical application and popularization.

[0109] Finally, it should be pointed out that: the above only for the preferred embodiments of the present application, and not for limiting the protection scope of the present application. Any modification, equivalent replacement, improvement, etc. within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A rock physics forward modeling method suitable for formations where the pore fluid is carbon dioxide, characterized in that, include: An equation of state for carbon dioxide under reservoir temperature and pressure conditions is introduced, and based on the equation of state, the dynamic response relationships of the physicochemical properties of pure carbon dioxide in gaseous, liquid and supercritical states as a function of temperature and pressure are obtained as constraints on the physical properties of carbon dioxide. Based on the aforementioned carbon dioxide physical property constraints, and according to the logging data, core analysis, and laboratory test results of the target well, formation evaluation is carried out to calculate reservoir parameters, which include clay content and overall porosity. The sandstone content is calculated based on the reservoir parameters, and a rock physical volume model of the target well is constructed using the sandstone content and the reservoir parameters. Based on the aforementioned rock physical volume model, for different carbon dioxide saturations, the equivalent medium theory is used to perform corresponding rock physical forward modeling to generate formation elastic property parameters, and the relationship between the corresponding saturation and the formation elastic property parameters is quantitatively characterized.

2. The rock physics forward modeling method as described in claim 1, characterized in that, The critical point of the supercritical state is a pressure of 7.38 MPa and a temperature of 31.1 °C. And / or, the equation of state adopts the Span–Wagner equation, wherein the Span–Wagner equation is a multi-parameter equation of state based on the Helmholtz free energy, and is expressed as follows: In the formula, A(ρ,T) represents the Helmholtz free energy. This represents the free energy component when the fluid is assumed to be an ideal gas. This represents the difference in free energy between a real fluid and an ideal gas, where T represents absolute temperature, R represents the gas constant, and δ represents the reduced density, and δ = ρ ÷ ρ c ρ represents density, ρ c = 467.6 kg·m -3 τ represents the reduced reciprocal temperature and τ = T c ÷T,T c = 304.1282 K, ln(δ) represents the compressibility factor relationship for an ideal gas. This represents the temperature dependence of the isobaric heat capacity of an ideal gas. Represents a constant term. Represents the linear temperature term. Represents the logarithmic temperature term. Indicates the corresponding vibration temperature and is used to control the starting point for this term's variation with temperature, n i The coefficients determine the weight of each item, a. i b represents the first power exponent and is used to control the dependence of the non-analytic term in the density direction. i d represents the second power exponent and is used to control the degree of singularity of this term on the critical point distance function. i t represents the third power exponent and is used to control the sensitivity of this term to density. i c represents the fourth power exponent and is used to control the sensitivity of this term to temperature. i α represents the attenuation factor. i Controlling the width of the Gaussian function in the δ direction, β i Controlling the width of the Gaussian function in the τ direction, ∈ i This indicates the position of the Gaussian center point in the δ direction, γ i Let represent the position of the Gaussian center point in the τ direction, e represent the base of the natural logarithm, and Δ represent the critical zone distance function.

3. The rock physics forward modeling method as described in claim 1, characterized in that, The calculation methods for the mud content include the following: Natural gamma ray logging is used to determine the clay content of rocks. Specifically, the natural gamma ray logging method employs the following linear interpolation formula to calculate the clay content: Among them, V sh The GR indicates the clay content. log GR represents the gamma value of the measuring point. min GR represents the gamma value of pure sand. max This indicates the gamma value of pure mud.

4. The rock physics forward modeling method as described in claim 1, characterized in that, The calculation method for the overall porosity includes the following: The first porosity, which reflects the difference between the total rock density and the matrix and pore fluid density, is calculated based on the density logging curve and is used as the density porosity. The density porosity is then corrected using the following formula: In the formula, PHIDC represents the corrected density porosity, PHID represents the density porosity under the corresponding skeleton, VCLAY represents the clay content, and ρ 2.76 ρ represents the density corresponding to dry clay. 2.65 ρ represents the density of sandstone, and ρ1 represents the density of pore fluid in the formation; A second porosity, reflecting the influence of pore size on the propagation velocity of sound waves, is calculated based on the acoustic logging curve and used as the acoustic porosity. This acoustic porosity is then corrected using a formula: In the formula, PHIDTC represents the corrected acoustic porosity, PHIDT represents the acoustic porosity under the corresponding skeleton, DT_VCL represents the acoustic transit time corresponding to dry clay, and γ 56 The longitudinal wave transit time of sandstone, γ 189 This represents the P-wave transit time of pore fluids in the formation. A third porosity, reflecting the hydrogen atom content in the rock, is calculated based on neutron logging curves and used as neutron porosity. This neutron porosity is then corrected using the following formula: In the formula, PHINC represents the corrected neutron porosity, PHIN represents the neutron porosity under the corresponding framework, and δ -0.02 δ represents the neutron porosity of sandstone. 0.36 δ1 represents the neutron porosity corresponding to dry clay, and δ1 represents the hydrogen index of water. Neutron-density porosity is obtained by correction using the following formula: PHINDC=PHIND-VCLAY×δ 0.18 In the formula, PHINDC represents the corrected neutron-density porosity, PHIND represents the neutron-density porosity under the corresponding framework, and PHIND = PHIN ÷ PHID, δ 0.18 This represents the neutron porosity of typical clay. The comprehensive porosity PHIT is calculated based on the density porosity, the acoustic porosity, the neutron porosity, and the neutron-density porosity according to the following formula: PHIT=PHIDC×ω1+PHINDC×ω2+PHINC×ω3+PHIDTC×ω4 In the formula, ω1, ω2, ω3 and ω4 represent the preset weight coefficients and ω1+ω2+ω3+ω4=1.

5. The rock physics forward modeling method as described in claim 1, characterized in that, The sandstone content calculated based on the reservoir parameters includes: Based on the clay content VCLAY and the comprehensive porosity PHIT in the reservoir parameters, the sandstone content VQUA = 1 - PHIT - VCLAY is calculated.

6. The rock physics forward modeling method as described in claim 1, characterized in that, The carbon dioxide saturation value ranges from [0% to 100%], and the different carbon dioxide saturations are set at intervals of 0%, 20%, 40%, 60%, 80%, and 100%.

7. The rock physics forward modeling method as described in claim 1, characterized in that, The formation elastic property parameters include density, P-wave transit time, S-wave transit time, and / or P-wave / S-wave velocity ratio.

8. The rock physics forward modeling method as described in claim 1, characterized in that, The equivalent medium theory posits that the macroscopic properties of rocks are jointly determined by the properties of a skeleton model, a dry rock model incorporating pores, and a saturated rock model incorporating fluids. The skeleton model is constructed based on the arithmetic mean of the Voigt upper bound and the Reuss lower bound of the Voigt-Reuss-Hill model, and is expressed as follows: In the formula, M VRH M represents the average modulus of Hill. V Voigt represents the upper limit of the equivalent elastic modulus and has N represents the total number of components, i represents the component number, and f i M represents the volume component of the i-th component. i M represents the elastic modulus of the i-th component. R This represents the Reuss lower limit of the equivalent elastic modulus and has 9. The rock physics forward modeling method as described in claim 8, characterized in that, The dry rock model is constructed based on the differential equivalent medium theory. This theory simulates a two-phase mixture by gradually adding inclusion phases into the solid mineral phase. Specifically, it employs an iterative calculation algorithm to incrementally embed pore structures of different geometries into the mineral framework phase. The equivalent bulk modulus K... * and shear modulus μ * The coupled differential equations are expressed as follows: In the formula, the equivalent bulk modulus K * initial value K * (0) = K1, the shear modulus μ * initial value μ * (0) = μ1, K1 represents the bulk modulus of the initial principal material, μ1 represents the shear modulus of the initial principal material, K2 represents the bulk modulus of the gradually added inclusions, the gradually added inclusions being pores of different geometries, μ2 represents the shear modulus of the gradually added inclusions, y represents the porosity of the gradually added inclusions, K * (y) represents the equivalent bulk modulus K when the porosity is y. * μ * (y) represents the shear modulus μ when the porosity is y. * P (*2) (y) represents the porosity of y and for the condition having an equivalent bulk modulus K. * and shear modulus μ * The first geometric factor, Q, of the inclusion material in the background medium. (*2) (y) represents the porosity of y and for the condition having an equivalent bulk modulus K. * and shear modulus μ * The second geometric factor of the containing material in the background medium.

10. The rock physics forward modeling method as described in claim 8, characterized in that, The saturated rock model is constructed based on the Gassmann fluid substitution equation, which describes the modulus change from a dry rock state to a saturated fluid porosity state, and is expressed as follows: In the formula, K sat1 K represents the effective bulk modulus of saturated rock for the first saturated fluid. sat2 K represents the effective bulk modulus of saturated rock for the second saturated fluid, K0 represents the bulk modulus of the minerals composing the rock, and K fl1 K represents the effective bulk modulus of the fluid in the first pore. fl2 This represents the effective bulk modulus of the fluid in the second pore. Indicates porosity.

Images