Method for calculating hydrate saturation in a well based on a multi-scale rock physics model

Abstract

The application belongs to the field of applied geophysical logging, and relates to a hydrate saturation calculation method in a well based on a multi-scale rock physics model, which considers multiple scale attenuation mechanisms and constructs a multi-scale rock physics model; according to logging data of a work area, the volume modulus and shear modulus of a solid phase are calculated, the volume modulus and shear modulus of the solid phase are substituted into the multi-scale rock physics model, the longitudinal wave velocity and attenuation coefficient of the reservoir in the work area are obtained, and compared with the longitudinal wave velocity and attenuation coefficient calculated by acoustic logging, the hydrate saturation of the reservoir in the work area is determined. The multi-scale rock physics model constructed by the application can finely depict the longitudinal wave velocity and attenuation characteristics of hydrate-bearing sediments, the calculated velocity and attenuation are accurate, and the saturation is calculated based on acoustic parameters and the multi-scale rock physics model, which provides a theoretical basis for quantitative interpretation of the saturation of hydrate-bearing sediments.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of applied geophysical well logging, and particularly relates to a hydrate saturation calculation method in a well based on a multi-scale rock physics model. BACKGROUND

[0002] Natural gas hydrate in nature widely exists in permafrost zones and marine sediments in the outer periphery of land. Because natural gas hydrate has the characteristics of wide distribution, huge resource quantity and high density, it has been widely concerned by many countries since the end of the 1960s and has been praised as a new type of potential energy with clean, efficient and abundant reserves. In the current hydrate reservoir logging interpretation and evaluation, the elastic wave rock physics model has an important position. However, the current existing elastic wave rock physics model cannot accurately describe the attenuation mechanism of the actual hydrate reservoir, which affects the calculation of hydrate saturation and brings great challenges to the quantitative geophysical analysis of hydrate-bearing sediments.

[0003] Wave-induced fluid flow (hereinafter referred to as: WIFF) is the main reason for the dispersion and attenuation of elastic waves in hydrate-bearing sediments, and mesoscopic and microscopic heterogeneity is the main mechanism causing WIFF. In hydrate-bearing sediments, mesoscopic and microscopic heterogeneity exists at the same time and can cause a significant change in the velocity of longitudinal waves, which means that it is necessary to consider the influence of the two mechanisms on dispersion and attenuation at the same time. Although the existing models (for example: BISQ model) have achieved certain application effect in simulating the velocity and attenuation characteristics of hydrate-bearing sediments, these models cannot simultaneously describe the attenuation mechanisms of the three scales of micro, meso and macro, which leads to the inability to accurately depict the longitudinal wave velocity and attenuation characteristics of hydrate-bearing sediments and seriously affects the calculation of hydrate saturation in the well. SUMMARY

[0004] The application provides a hydrate saturation calculation method in a well based on a multi-scale rock physics model, which solves the problem that the prior art cannot accurately describe the longitudinal wave velocity and attenuation of hydrate-bearing sediments. The method constructs a multi-scale rock physics model for describing the longitudinal wave velocity and attenuation of actual natural gas hydrate reservoirs, forms a corresponding modeling method, and calculates the saturation based on acoustic parameters on this basis. The multi-scale rock physics model constructed by the method can accurately describe the longitudinal wave velocity and attenuation characteristics of hydrate-bearing sediments in different occurrence forms, and provides a theoretical basis for the quantitative interpretation of hydrate saturation of hydrate-bearing sediments.

[0005] In order to achieve the above purpose, the application provides a hydrate saturation calculation method in a well based on a multi-scale rock physics model, and the steps are as follows:

[0006] S1, a multi-scale rock physics model is established, and the specific steps are as follows:

[0007] S11, the jet flow model is added to the Kuster-Toksuez model to calculate the hydrate equivalent bulk modulus and the hydrate equivalent shear modulus of the gas / water inclusion;

[0008] S12, the bulk modulus of the pore water and the methane gas and the equivalent bulk modulus of the hydrate obtained in step S11 are substituted into the wood formula to obtain the bulk modulus of the first fluid phase, which is expressed as:

[0009]

[0010] wherein, is the bulk modulus of the first fluid phase, is the bulk modulus of the pore water, is the bulk modulus of the methane gas, is the equivalent bulk modulus of the hydrate, is the saturation of the pore water, is the saturation of the methane gas, is the saturation of the hydrate;

[0011] The bulk modulus of the pore water and the methane gas is substituted into the wood formula to obtain the bulk modulus of the second fluid phase, which is expressed as:

[0012]

[0013] The bulk modulus and the shear modulus of the quartz particles and the equivalent bulk modulus and the equivalent shear modulus of the hydrate obtained in step S11 are substituted into the Hill average equation to obtain the bulk modulus and the shear modulus of the solid phase, which are expressed as:

[0014]

[0015]

[0016] wherein, K is the bulk modulus of the solid phase, G is the shear modulus of the solid phase; K q is the bulk modulus of the quartz particles, G q is the shear modulus of the quartz particles, f H is the volume fraction of the hydrate in the solid phase, f q is the volume fraction of the hydrate in the solid phase;

[0017] S13, the equivalent volume modulus and equivalent shear modulus of the hydrate obtained in step S11 are substituted into the cemented sandstone model to obtain the volume modulus and shear modulus of the dry rock skeleton in contact with the cemented hydrate; the volume modulus and shear modulus of the quartz particles are substituted into the cemented sandstone model to obtain the volume modulus and shear modulus of the dry rock skeleton in contact with the particle-coated hydrate; the volume modulus and shear modulus of the solid phase obtained in step S12 are substituted into the Hashin-Shtrikman model to obtain the volume modulus and shear modulus of the dry rock skeleton in contact with the particle-supported hydrate; and the volume modulus and shear modulus of the quartz particles are substituted into the Hashin-Shtrikman model to obtain the volume modulus and shear modulus of the dry rock skeleton in contact with the pore-filled hydrate;

[0018] S14, the volume modulus and shear modulus of the fluid phase obtained in step S12 and the volume modulus and shear modulus of the dry rock skeleton obtained in step S13 are combined to obtain the longitudinal wave velocity and attenuation coefficient of the hydrate-bearing sediment in four different occurrence forms, i.e., the contact cemented hydrate, the particle-coated hydrate, the particle-supported hydrate, and the pore-filled hydrate, according to the Biot-Rayleigh theoretical model, i.e., the multi-scale rock physics model;

[0019] S2, the volume fraction of the rock composition is obtained according to the natural gamma ray logging, porosity logging, and lithology density logging data of the target layer of the work area, the volume modulus and shear modulus of the solid phase are calculated according to the Hill average equation, the saturation of the hydrate is obtained according to the resistivity logging data of the work area, then the modulus of the solid phase and the saturation of the hydrate are substituted into the multi-scale rock physics model constructed in step S1, and finally the longitudinal wave velocity and attenuation coefficient of the reservoir of the work area are obtained. Compared with the longitudinal wave velocity and attenuation coefficient obtained by the measured acoustic logging, the saturation of the hydrate is updated until the predicted and measured longitudinal wave velocity and attenuation are within the set error range, and the predicted saturation of the hydrate of the reservoir of the work area is obtained.

[0020] Preferably, in step S11, the jet flow model is represented as:

[0021]

[0022] In the formula, is the volume modulus of the inclusion, is the shear modulus of the inclusion, is the equivalent volume modulus of the fluid in the hydrate inclusion and the pore, and the subscript , N is the different types of inclusions; is the volume modulus of the hydrate, and the subscript H is the hydrate, is the relaxation time, ωfor the angular frequency, for the viscosity of the inclusion in the hydrate;

[0023] The equivalent bulk modulus and the equivalent shear modulus of the hydrate are calculated by adding the jet flow model to the Kuster-Toksöz model, and are expressed as:

[0024]

[0025] wherein, is the equivalent bulk modulus of the hydrate, is the equivalent shear modulus of the hydrate, is the shear modulus of the hydrate, is the first i volume fraction of the inclusion in the hydrate; , is a geometric factor, , α is the radius of the water and methane gas inclusion.

[0026] Preferably, in step S13, the cemented sandstone model is expressed as:

[0027]

[0028]

[0029] wherein, is the bulk modulus of the dry rock, is the shear modulus of the dry rock, is the bulk modulus of the cement or the hydrate, is the shear modulus of the quartz grain or the hydrate; the parameter is proportional to the positive direction of the combination of two grains with cement, is proportional to the shear stiffness of the combination of two grains with cement; the parameter and depend on the content of the cement and the characteristics of the cement and the skeleton grain; n is the coordination number of the grain, and in the hydrate-bearing sediment n = 8.5; is the initial porosity of the sediment;

[0030] When the porosity is less than the critical porosity, the Hashin-Shtrikman model is expressed as:

[0031]

[0032]

[0033] When the porosity is greater than or equal to the critical porosity, the Hashin-Shtrikman model is expressed as:

[0034]

[0035]

[0036] wherein, is the porosity; is the critical porosity, is the bulk modulus at the critical porosity, is the shear modulus at the critical porosity, is the Poisson's ratio of the mineral phase calculated from K and G K is the bulk modulus of the solid phase or quartz grains, G is the shear modulus of the solid phase or quartz grains, P is the equivalent pressure, .

[0037] Preferably, in step S14:

[0038] the bulk modulus and the shear modulus of the second fluid phase obtained in step S12 and the bulk modulus and the shear modulus of the dry rock matrix of the contact cement hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and the attenuation coefficient of the contact cement hydrate deposit;

[0039] the bulk modulus and the shear modulus of the second fluid phase obtained in step S12 and the bulk modulus and the shear modulus of the dry rock matrix of the grain-coated hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and the attenuation coefficient of the grain-coated hydrate deposit;

[0040] the bulk modulus and the shear modulus of the second fluid phase obtained in step S12 and the bulk modulus and the shear modulus of the dry rock matrix of the grain-supported hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and the attenuation coefficient of the grain-supported hydrate deposit;

[0041] the bulk modulus and the shear modulus of the first fluid phase obtained in step S12 and the bulk modulus and the shear modulus of the dry rock matrix of the pore-filling hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and the attenuation coefficient of the pore-filling hydrate deposit.

[0042] ​Preferably, the wave equation of the Biot-Rayleigh double-pore theoretical model is:

[0043]

[0044]

[0045]

[0046]

[0047] wherein, N , A , , , i =1,2 is a rigidity coefficient, u is the average particle displacement of the solid, , are the first and second derivatives of u respectively, , , denote the divergence of the displacement field of the solid, the fluid in pore 1 and the fluid in pore 2 respectively; , i =1,2,3, j =1,2,3 is a density parameter, is the density of the host fluid, , are the porosities of the water-saturated pore and the hydrate-saturated pore respectively, , are the local porosities of the two regions, , i =1,2 is a dissipation parameter, is the increment of fluid strain caused by the local flow process, , are the first and second derivatives of respectively, is the viscosity of the host fluid, is the permeability of the host fluid; if the hydrate is the pore fluid, is the displacement when the first fluid phase is the formation water, is the displacement when the second fluid phase is the natural gas hydrate / free gas; if the hydrate is the rock skeleton, is the displacement when the first fluid phase is the fluid in the host skeleton, is the displacement when the second fluid phase is the fluid in the hydrate inclusion;

[0048] According to the plane wave analysis method, the above equation is solved by bringing in the analytical solution of the plane longitudinal wave, to obtain:

[0049]

[0050] wherein, represents the wave number of the contact cement hydrate or the particle coating hydrate or the particle support hydrate or the pore filling hydrate; 、 represents the coefficient matrix of the equation, i =1,2,3, j =1,2,3;

[0051] Then the calculation formula of the velocity and the attenuation is represented as:

[0052]

[0053]

[0054] wherein, is the attenuation coefficient; v is the P-wave velocity of the hydrate-bearing sediment; and respectively take the real part and the imaginary part of k .

[0055] Compared with the prior art, the present application has the beneficial effects that:

[0056] (1) The hydrate saturation calculation method in the well based on the multi-scale rock physics model of the present application, the equivalent elastic modulus of the hydrate is calculated based on the jet flow model according to the Kuster-Toksöz model; the elastic modulus of the fluid phase and the solid phase is respectively calculated according to the wood formula and the Hill average equation; the dry rock skeleton elastic modulus of the contact cement hydrate and the particle coating hydrate is calculated by using the cemented sandstone model, and the dry rock skeleton elastic modulus of the pore filling hydrate and the particle support hydrate is calculated by using the Hashin-Shtrikman model; the elastic modulus of the fluid phase and the dry rock skeleton is combined, and the P-wave velocity and the attenuation coefficient of the four kinds of hydrate-bearing sediment in the four kinds of occurrence modes are obtained according to the Biot-Rayleigh theoretical model, that is, the multi-scale rock physics model. The multi-scale rock physics model constructed by the present application considers the attenuation mechanism of multiple scales, compared with the existing theoretical model, is more consistent with the actual situation of the hydrate-bearing sediment, can depict the acoustic response law of the hydrate-bearing sediment reservoir, and then the hydrate-bearing sediment can be detected and identified by using the acoustic logging.

[0057] (2) The hydrate saturation calculation method in the well based on the multi-scale rock physics model of the present application considers four different hydrate occurrence modes, which is consistent with the actual hydrate-bearing sediment reservoir occurrence mode, can accurately describe the P-wave velocity and the attenuation coefficient of the hydrate-bearing sediment in different occurrence modes, and provides a theoretical basis for the quantitative interpretation and evaluation of the hydrate-bearing sediment.

[0058] (3) The multi-scale rock physical model based on the hydrate saturation calculation method in the well provided by the application can describe the P-wave velocity and attenuation coefficient of hydrate-bearing sediments in different occurrence forms, and the multi-scale rock physical model constructed by the method is applied to hydrate detection, thereby providing a new method for hydrate detection by using acoustic parameters (i.e., P-wave velocity and attenuation coefficient), and can provide a measurement method for improvement and perfection of acoustic logging.

[0059] (4) The multi-scale rock physical model based on the hydrate saturation calculation method in the well provided by the application uses the multi-scale rock physical model constructed by the application, and uses P-wave velocity and attenuation to determine hydrate saturation, compared with the Archie formula or the resistivity saturation calculation method, the new saturation calculation method based on acoustic parameters provided by the application avoids the inapplicable problem of the Archie formula in argillaceous formations or weakly lithified formations, greatly expands the saturation calculation method of hydrate reservoirs, and to some extent, solves the evaluation and interpretation of hydrate saturation in weakly lithified formations. BRIEF DESCRIPTION OF DRAWINGS

[0060] Figure 1 The flowchart of the multi-scale rock physical model based on the hydrate saturation calculation method in the well is described in the embodiments of the application;

[0061] Figure 2a The P-wave velocity dispersion curve of the pore-filling hydrate deposit is shown in the schematic diagram;

[0062] Figure 2b The attenuation coefficient dispersion curve of the pore-filling hydrate deposit is shown in the schematic diagram;

[0063] Figure 3a The P-wave velocity curve of different hydrate occurrence forms with hydrate saturation is shown in the diagram;

[0064] Figure 3b The attenuation coefficient curve of different hydrate occurrence forms with hydrate saturation is shown in the diagram. DETAILED DESCRIPTION

[0065] In the following, the application will be specifically described through exemplary embodiments. However, it should be understood that the elements, structures and features in one embodiment can also be beneficially combined into other embodiments without further description.

[0066] Reference Figure 1 The embodiments of the application provide a multi-scale rock physical model based hydrate saturation calculation method in a well, and the specific steps are as follows:

[0067] S1, a multi-scale rock physics model is established, and the specific steps are as follows:

[0068] S11, the jet flow model is added to the Kuster-Toksuez model to calculate the equivalent bulk modulus and equivalent shear modulus of the hydrate containing gas / water inclusions.

[0069] Specifically, the jet flow model is represented as:

[0070]

[0071] In the formula, is the bulk modulus of the inclusion, is the shear modulus of the inclusion, is the equivalent bulk modulus of the fluid in the hydrate inclusion and the pore, subscript , N is the different kinds of inclusions; is the bulk modulus of the hydrate, subscript H is the hydrate, is the relaxation time, ω is the angular frequency, is the viscosity of the inclusion in the hydrate;

[0072] The jet flow model is added to the Kuster-Toksuez model to calculate the equivalent bulk modulus and equivalent shear modulus of the hydrate, which is represented as:

[0073]

[0074] In the formula, is the equivalent bulk modulus of the hydrate, is the equivalent shear modulus of the hydrate, is the shear modulus of the hydrate, is the volume fraction of the first i inclusion in the hydrate; , is the geometric factor, , α is the radius of the water and methane gas inclusions.

[0075] S12, the bulk modulus of the pore water and the methane gas and the equivalent bulk modulus of the hydrate obtained in step S11 are substituted into the wood formula to obtain the bulk modulus of the first fluid phase, which is represented as:

[0076]

[0077] In the formula, is the bulk modulus of the first fluid phase, is the bulk modulus of the pore water, is the bulk modulus of the methane gas, the equivalent bulk modulus of the hydrate, the saturation of the pore water, the saturation of the methane gas, the saturation of the hydrate;

[0078] Substituting the bulk modulus and the shear modulus of the quartz grains and the equivalent bulk modulus and the equivalent shear modulus of the hydrate obtained in step S11 into the Hill averaging equation, the bulk modulus and the shear modulus of the solid phase are expressed as:

[0079]

[0080] Substituting the bulk modulus and the shear modulus of the quartz grains and the equivalent bulk modulus and the equivalent shear modulus of the hydrate obtained in step S11 into the Hill averaging equation, the bulk modulus and the shear modulus of the solid phase are expressed as:

[0081]

[0082]

[0083] wherein, K the bulk modulus of the solid phase, G the shear modulus of the solid phase; K q the bulk modulus of the quartz grains, G q the shear modulus of the quartz grains, f H the volume fraction of the hydrate in the solid phase, f q the volume fraction of the hydrate in the solid phase.

[0084] S13, substituting the equivalent bulk modulus and the equivalent shear modulus of the hydrate obtained in step S11 into the cemented sandstone model, the bulk modulus and the shear modulus of the dry rock skeleton contacting the cemented hydrate are obtained; substituting the bulk modulus and the shear modulus of the quartz grains into the cemented sandstone model, the bulk modulus and the shear modulus of the dry rock skeleton of the grain-coated hydrate are obtained; substituting the bulk modulus and the shear modulus of the solid phase obtained in step S12 into the Hashin-Shtrikman model, the bulk modulus and the shear modulus of the dry rock skeleton of the grain-supported hydrate are obtained; substituting the bulk modulus and the shear modulus of the quartz grains into the Hashin-Shtrikman model, the bulk modulus and the shear modulus of the dry rock skeleton of the pore-filled hydrate are obtained.

[0085] It is to be noted that when the hydrate is part of the pore fluid, it is added to the fluid phase to form a pore-filling hydrate. When the hydrate is a solid particle, it is added to the solid matrix to form a particle-supported hydrate. When the hydrate is a cement, a contact-cemented hydrate is formed when the cementation mode is that all the hydrate precipitates at the particle contacts, and a particle-coating hydrate is formed when the cementation mode is that the hydrate precipitates uniformly on the particle surfaces.

[0086] Two cementation modes and parameters α When the cementation mode is that all the hydrate precipitates at the particle contacts, When the cementation mode is that the hydrate precipitates uniformly on the particle surfaces, wherein, S is the saturation of the cement (hydrate) in the pore space, n is the coordination number of the particles, taken in the hydrate-bearing sediment, n = 8.5; is the initial porosity of the sediment.

[0087] In particular, the cemented sandstone model is expressed as:

[0088]

[0089]

[0090] wherein, is the bulk modulus of the dry rock, is the shear modulus of the dry rock, is the bulk modulus of the cement or hydrate, is the shear modulus of the quartz particles or hydrate; the parameters are proportional to the normal stiffness of the two-particle assembly with cement, are proportional to the shear stiffness of the two-particle assembly with cement, the parameters and depend on the content of the cement and the properties of the cement and the backbone particles; n is the coordination number of the particles, taken in the hydrate-bearing sediment, n = 8.5; is the initial porosity of the sediment.

[0091] It is to be noted that when the equivalent bulk modulus and the equivalent shear modulus of the hydrate obtained in step S11 are substituted into the cemented sandstone model, the bulk modulus of the dry rock skeleton and the shear modulus of the dry rock skeleton of the contact-cemented hydrate are obtained, is the equivalent bulk modulus of the hydrate , is the equivalent shear modulus of the hydrate The bulk modulus and the shear modulus of the quartz grains are substituted into the cemented sandstone model to obtain the dry rock skeleton bulk modulus and the dry rock skeleton shear modulus of the grain-coated hydrate, is the bulk modulus of the quartz grains K q , is the shear modulus of the quartz grains G q .

[0092] Specifically, when the porosity is less than the critical porosity, the Hashin-Shtrikman model is expressed as:

[0093]

[0094]

[0095] When the porosity is greater than or equal to the critical porosity, the Hashin-Shtrikman model is expressed as:

[0096]

[0097]

[0098] wherein, is the porosity; is the critical porosity, is the bulk modulus at the critical porosity, is the shear modulus at the critical porosity, is the Poisson's ratio of the mineral phase calculated from K and G is the Poisson's ratio of the mineral phase calculated from K is the bulk modulus of the solid phase or the quartz grains, G is the shear modulus of the solid phase or the quartz grains, P is the equivalent pressure, .

[0099] It should be noted that when the bulk modulus and the shear modulus of the solid phase obtained in step S12 are substituted into the Hashin-Shtrikman model to obtain the dry rock skeleton bulk modulus and the dry rock skeleton shear modulus of the grain-supported hydrate, the parameters in the Hashin-Shtrikman model are K is the bulk modulus of the solid phase K , G is the shear modulus of the solid phase GThe volume modulus and shear modulus of the quartz particles are substituted into the Hashin-Shtrikman model to obtain the volume modulus and shear modulus of the dry rock skeleton of the pore-filling hydrate, and the parameter of the Hashin-Shtrikman model is K the volume modulus of the quartz particles K q , G the shear modulus of the quartz particles G q .

[0100] S14, the volume modulus and shear modulus of the fluid phase obtained in step S12 and the volume modulus and shear modulus of the dry rock skeleton obtained in step S13 are combined to obtain the longitudinal wave velocity and attenuation coefficient of the hydrate-bearing sediment in four different occurrence forms of contact-cemented hydrate, particle-coating hydrate, particle-supporting hydrate, and pore-filling hydrate according to the Biot-Rayleigh theoretical model, that is, a multi-scale rock physics model.

[0101] Specifically, the volume modulus and shear modulus of the second fluid phase obtained in step S12 and the volume modulus and shear modulus of the dry rock skeleton of the contact-cemented hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and attenuation coefficient of the contact-cemented hydrate sediment.

[0102] Specifically, the volume modulus and shear modulus of the second fluid phase obtained in step S12 and the volume modulus and shear modulus of the dry rock skeleton of the particle-coating hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and attenuation coefficient of the particle-coating hydrate sediment.

[0103] Specifically, the volume modulus and shear modulus of the second fluid phase obtained in step S12 and the volume modulus and shear modulus of the dry rock skeleton of the particle-supporting hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and attenuation coefficient of the particle-supporting hydrate sediment.

[0104] Specifically, the volume modulus and shear modulus of the first fluid phase obtained in step S12 and the volume modulus and shear modulus of the dry rock skeleton of the pore-filling hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and attenuation coefficient of the pore-filling hydrate sediment.

[0105] Specifically, the wave equation of the Biot-Rayleigh double-pore theoretical model is:

[0106]

[0107]

[0108]

[0109]

[0110] where, N , A , , , i =1,2 are stiffness coefficients, u is the solid average particle displacement, , are the first and second derivatives of u respectively, , , denote the divergence of the displacement field of the solid, pore 1 fluid and pore 2 fluid respectively; , i =1,2,3, j =1,2,3 are density parameters, is the density of the host phase fluid, , are the water-saturated and hydrate-saturated porosities respectively, , are the local porosities of the two regions, , i =1,2 are dissipation parameters, is the fluid strain increment due to local flow processes, , are the first and second derivatives of respectively, is the viscosity of the host phase fluid, is the permeability of the host phase fluid; if hydrate is the pore fluid, is the displacement when the first fluid phase is formation water, is the displacement when the second fluid phase is natural gas hydrate / free gas; if hydrate is the rock skeleton, is the displacement when the first fluid phase is fluid in the host skeleton, is the displacement when the second fluid phase is fluid in the hydrate inclusions;

[0111] The above equations are solved according to the plane wave analysis method, with the analytical solution of the plane longitudinal wave, to obtain:

[0112]

[0113] where, kω represents the wave number of the contact cementation hydrate or the particle coating hydrate or the particle support hydrate or the pore filling hydrate; , C represents the coefficient matrix of the equation, i = 1, 2, 3, j = 1, 2, 3;

[0114] When ω = 1, k ω represents the wave number of the contact cementation hydrate k 1, the calculation formula of the velocity and the attenuation is represented as:

[0115]

[0116]

[0117] wherein, is the attenuation coefficient; v 1 is the longitudinal wave velocity of the contact cementation hydrate deposit; and respectively take the real part and the imaginary part of k 1.

[0118] When ω = 2, k ω represents the wave number of the particle coating hydrate k 2, the calculation formula of the velocity and the attenuation is represented as:

[0119]

[0120]

[0121] wherein, is the attenuation coefficient; v 2 is the longitudinal wave velocity of the particle coating hydrate deposit; and respectively take the real part and the imaginary part of k 2.

[0122] When ω = 3, k ω represents the wave number of the particle support hydrate k 3, the calculation formula of the velocity and the attenuation is represented as:

[0123]

[0124]

[0125] wherein, is the attenuation coefficient; v 3 is the longitudinal wave velocity of the particle support hydrate deposit; and respectively take the real part and the imaginary part of k 3.

[0126] When k Wave number representing pore-filling hydrate k 4, the formula for calculating velocity and attenuation is expressed as:

[0127]

[0128]

[0129] wherein, is the attenuation coefficient; v 4 is the P-wave velocity of the pore-filling hydrate deposit; and respectively take k the real part and the imaginary part of

[0130] S2, the volume fractions of rock components (calcite, illite and quartz) are obtained according to the gamma ray logging, porosity logging and density logging data of the work area, the bulk modulus and shear modulus of the solid phase are calculated according to the Hill average equation, the saturation of hydrate is obtained according to the resistivity logging data of the work area, then the modulus of the solid phase and the saturation of hydrate are substituted into the multi-scale rock physics model constructed in the S1 step, and finally the P-wave velocity and the attenuation coefficient of the reservoir in the work area are obtained. Compared with the P-wave velocity and the attenuation coefficient obtained from the measured acoustic logging data, the saturation of hydrate is updated until the predicted and measured P-wave velocity and attenuation are within the set error range, and the predicted saturation of hydrate in the reservoir in the work area is obtained.

[0131] Specifically, the specific method for updating the saturation of hydrate is: a least square objective function is constructed by the velocity and attenuation obtained from the multi-scale rock physics model and the velocity and attenuation coefficient calculated from the measured acoustic logging data, and the water saturation corresponding to the minimum least square of the constructed objective function is the saturation of hydrate in the target layer when the least square of the constructed objective function reaches the minimum by continuously adjusting each parameter in the multi-scale rock physics model. It should be noted that some parameters are obtained by experiments on regional strata and substituted into the multi-scale rock physics model, and some are different, such as the saturation of hydrate.

[0132] In the above method, the multi-scale petrophysical model of the hydrate-bearing sediment constructed considers multiple scale attenuation mechanisms, is more in line with the actual situation of the hydrate-bearing sediment compared with the existing theoretical model, can more accurately describe the P-wave velocity and attenuation coefficient of the hydrate-bearing sediment in different occurrence forms, and the calculated velocity and attenuation are more accurate, which has important significance and value for the acoustic detection and identification of the hydrate-bearing sediment, and provides a theoretical basis for the saturation quantitative interpretation of the hydrate-bearing sediment. The multi-scale petrophysical model is applied to hydrate detection, and then a new method for hydrate detection by using acoustic parameters (i.e. P-wave velocity and attenuation coefficient) is provided, which can provide a measurement method for the improvement and perfection of acoustic logging.

[0133] In order to illustrate the effect of the multi-scale petrophysical model constructed by the above method of the present application. Figure 2a 、 Figure 2b 、 Figure 3a 、 Figure 3b The frequency dispersion curves of the P-wave velocity and attenuation coefficient in the hydrate-bearing sediment and the variation trend of the P-wave velocity and attenuation coefficient with the hydrate saturation in the four different hydrate occurrence forms are given, and are compared with the experimental data.

[0134] Figure 2a 、 Figure 2b The P-wave velocity and P-wave attenuation coefficient as functions of frequency for different inclusion aspect ratios are shown respectively. The multi-scale petrophysical model constructed by the above method of the present application predicts three relaxation peaks, i.e. local flow (mesoscale), global flow (macro scale) and jet flow (micro scale). By analyzing Figure 2a 、 Figure 2b The following conclusions can be drawn: with the increase of the aspect ratio of the water (methane gas) inclusion, the third peak moves to a higher frequency, the aspect ratio controls the jet flow relaxation time, and then affects the position of the attenuation peak.

[0135] Figure 3a 、 Figure 3b The variation of the P-wave velocity and attenuation coefficient with the hydrate saturation in the four different hydrate occurrence forms is shown respectively. It can be seen from Figure 3a 、 Figure 3b that the "excess water" method mainly generates pore-filling hydrate, and the "excess gas" method mainly generates cemented hydrate, and the simulation result is more consistent with the particle coating hydrate. When the hydrate saturation is below 0%-40%, the multi-scale petrophysical model constructed by the above method of the present application can capture and measure the variation trend of the P-wave velocity and attenuation coefficient with the hydrate saturation.

[0136] The above examples are used to explain the present application, but not to limit the present application, any modification and change made to the present application within the spirit and protection scope of the claims, fall into the protection scope of the present application.

Claims

1. A method for calculating hydrate saturation in a well based on a multi-scale rock physics model, characterized in that, The steps are: S1, establishing a multi-scale rock physics model, the specific steps are: S11, adding the jet flow model to the Kuster-Toksouz model to calculate the equivalent bulk modulus and equivalent shear modulus of hydrate containing gas / water inclusions; S12, substituting the bulk modulus of pore water and methane gas and the equivalent bulk modulus of hydrate obtained in step S11 into the wood formula to obtain the bulk modulus of the first fluid phase, which is expressed as: wherein is the bulk modulus of the first fluid phase, is the bulk modulus of the pore water, is the bulk modulus of the methane gas, is the equivalent bulk modulus of the hydrate, is the saturation of the pore water, is the saturation of the methane gas, is the saturation of the hydrate; Substitute the bulk modulus of pore water and methane gas into the wood formula to obtain the bulk modulus of the second fluid phase, which is expressed as: Substitute the bulk modulus and shear modulus of quartz particles and the equivalent bulk modulus and equivalent shear modulus of hydrate obtained in step S11 into the Hill average equation to obtain the bulk modulus and shear modulus of the solid phase, which is expressed as: wherein K is the bulk modulus of the solid phase, G is the shear modulus of the solid phase; K q is the bulk modulus of the quartz grains, G q is the shear modulus of the quartz grains, f H is the volume fraction of hydrates in the solid phase, is the equivalent shear modulus of the hydrates, f q is the volume fraction of hydrates in the solid phase; S13, substituting the equivalent bulk modulus and equivalent shear modulus of hydrate obtained in step S11 into the cemented sandstone model to obtain the bulk modulus and shear modulus of the dry rock skeleton of contact cemented hydrate; substituting the bulk modulus and shear modulus of quartz particles into the cemented sandstone model to obtain the bulk modulus and shear modulus of the dry rock skeleton of particle coating hydrate; substituting the bulk modulus and shear modulus of the solid phase obtained in step S12 into the Hashin-Shtrikman model to obtain the bulk modulus and shear modulus of the dry rock skeleton of particle supported hydrate; and substituting the bulk modulus and shear modulus of quartz particles into the Hashin-Shtrikman model to obtain the bulk modulus and shear modulus of the dry rock skeleton of pore filling hydrate; S14, combining the bulk modulus and shear modulus of the fluid phase obtained in step S12 and the bulk modulus and shear modulus of the dry rock skeleton obtained in step S13, and according to the Biot-Rayleigh theoretical model, the longitudinal wave velocity and attenuation coefficient of hydrate-containing sediment in four different occurrence forms of contact cemented hydrate, particle coating hydrate, particle supported hydrate and pore filling hydrate are obtained, that is, a multi-scale rock physics model; S2, obtaining the volume fraction of rock composition according to the natural gamma ray logging, porosity logging and lithology density logging data of the target layer of the work area, calculating the bulk modulus and shear modulus of the solid phase according to the Hill average equation, obtaining the saturation of hydrate according to the resistivity logging data of the work area, then substituting the modulus of the solid phase and the saturation of the hydrate into the multi-scale rock physics model constructed in step S1, finally obtaining the longitudinal wave velocity and attenuation coefficient of the reservoir in the work area, comparing with the longitudinal wave velocity and attenuation coefficient obtained by acoustic logging, updating the saturation of hydrate until the predicted and measured longitudinal wave velocity and attenuation are within the set error range, and obtaining the predicted saturation of hydrate in the reservoir in the work area.

2. The method for calculating hydrate saturation in a well based on a multi-scale rock physics model of claim 1, wherein, In step S11, the jet flow model is expressed as: wherein is the bulk modulus of the inclusion, is the shear modulus of the inclusion, is the equivalent bulk modulus of the fluid in the hydrate inclusion and the pores, subscript , N is the different kind of inclusion; is the bulk modulus of the hydrate, subscript H is the hydrate, is the relaxation time, ω is the angular frequency, is the viscosity of the inclusion in the hydrate; The equivalent bulk modulus and equivalent shear modulus of hydrate are calculated by adding the jet flow model to the Kuster-Toksouz model, which is expressed as: wherein K is the bulk modulus of the hydrate, G is the shear modulus of the hydrate, G is the shear modulus of the hydrate, K is the bulk modulus of the hydrate, i V is the volume fraction of the inclusions; , K is the bulk modulus of the hydrate, , α R is the radius of the water and methane gas inclusions.

3. The method for calculating hydrate saturation in a well based on a multi-scale rock physics model of claim 1, wherein, In step S13, the cemented sandstone model is expressed as: wherein is the bulk modulus of dry rock, is the shear modulus of dry rock, is the bulk modulus of cement or hydrate, is the shear modulus of quartz grains or hydrate; parameter is proportional to the positive aspect of the two-grain assembly with cement, is proportional to the shear stiffness of the two-grain assembly with cement; parameter and depend on the content of cement and on the properties of the cement and the skeleton grains; n is the coordination number of the grains, taken as n = 8.5; is the initial porosity of the deposit; When the porosity is less than the critical porosity, the Hashin-Shtrikman model is expressed as: When the porosity is greater than or equal to the critical porosity, the Hashin-Shtrikman model is expressed as: wherein is the porosity; is the critical porosity, is the bulk modulus at the critical porosity, is the shear modulus at the critical porosity, is the Poisson's ratio of the mineral phase calculated from K 'and G 'is the Poisson's ratio of the mineral phase calculated from K 'is the bulk modulus of the solid phase or quartz grains, G 'is the shear modulus of the solid phase or quartz grains, P is the equivalent pressure, .

4. The method for calculating hydrate saturation in a well based on a multi-scale rock physics model of claim 2, wherein, In step S14: The bulk modulus and shear modulus of the second fluid phase obtained in step S12 and the bulk modulus and shear modulus of the dry rock framework of the contact cement hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and attenuation coefficient of the contact cement hydrate deposit; The bulk modulus and shear modulus of the second fluid phase obtained in step S12 and the bulk modulus and shear modulus of the dry rock framework of the contact cement hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and attenuation coefficient of the contact cement hydrate deposit; The bulk modulus and shear modulus of the second fluid phase obtained in step S12 and the bulk modulus and shear modulus of the dry rock framework of the contact cement hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and attenuation coefficient of the contact cement hydrate deposit; The bulk modulus and shear modulus of the second fluid phase obtained in step S12 and the bulk modulus and shear modulus of the dry rock framework of the contact cement hydrate obtained in step S13 are substituted into the Biot-Rayleigh double-pore theoretical model to obtain the longitudinal wave velocity and attenuation coefficient of the contact cement hydrate deposit.

5. The method for calculating hydrate saturation in a well based on a multi-scale rock physics model of claim 4, wherein, The wave equation of the Biot-Rayleigh double-pore theoretical model is: where N , A , , , i =1,2 is the rigidity coefficient, u is the solid average particle displacement, , are the first and second derivatives of u , respectively, , , denote the divergence of the displacement field of the solid, pore 1 fluid and pore 2 fluid, respectively; , i =1,2,3, j =1,2,3 are the density parameters, is the density of the host phase fluid, , are the porosities of the water-saturated and hydrate-saturated pores, respectively, , are the local porosities of the two regions, , i =1,2 are the dissipation parameters, is the fluid strain increment due to local flow processes, , are the first and second derivatives of , respectively, is the viscosity of the host phase fluid, is the permeability of the host phase fluid; if hydrates are the pore fluid, is the displacement when the first fluid phase is formation water, is the displacement when the second fluid phase is natural gas hydrate / free gas; if hydrates are the rock skeleton, is the displacement when the first fluid phase is fluid in the host skeleton, is the displacement when the second fluid phase is fluid in the hydrate inclusions. According to the plane wave analysis method, the analytical solution of the plane longitudinal wave is substituted into the above equation to obtain: wherein represents the wave number of the contact cement hydrate or the particle coating hydrate or the particle support hydrate or the pore filling hydrate; , represents the coefficient matrix of the equation, i = 1,2,3, j = 1,2,3; The calculation formula of the velocity and attenuation is expressed as: The calculation formula of the velocity and attenuation is expressed as: wherein is the attenuation coefficient; v is the P-wave velocity of the hydrate-bearing sediment; and the real and imaginary parts of k respectively.

Images