A method for predicting the imbibition recovery rate of shale oil reservoirs

Abstract

The present invention discloses a method for predicting the imbibition recovery rate of shale oil reservoirs, comprising: (1) respectively regressing the oil-water relative permeability curves to obtain the oil-water relative permeability; calculating the average water saturation, initial water saturation, and residual oil saturation; obtaining the average oil-water relative permeability of the shale oil reservoir; (2) obtaining the fractional flow function and solving it; (3) obtaining the reverse imbibition water saturation profile; (4) establishing the mass conservation equations of shale fractures and matrix at the field scale for predicting the imbibition recovery rate. Through the prediction method established in the present invention, the deficiencies in the traditional method of only considering capillary force and ignoring boundary slip can be overcome, and the reverse imbibition length and imbibition capacity can be accurately determined, thereby providing a basis for predicting the field recovery rate and optimizing fracturing fluids, etc.

Description

Technical Field

[0001] The present invention belongs to the field of oil exploration and development, and more particularly, to a method for predicting the imbibition recovery factor of shale oil reservoirs. Background Art

[0002] During the development of shale oil, the spontaneous imbibition of water into the matrix blocks is an important development mechanism for enhancing oil recovery. In the existing technologies, the oil recovery factor is usually studied based on laboratory experiments, but this method is unstable and cannot predict the macroscopic oil recovery factor of the entire target area. At the same time, the existing research models cannot systematically predict the reverse imbibition recovery factor, and there are problems such as insufficient consideration of factors and inability to scale up to the field scale.

[0003] Therefore, there is an urgent need for a method for predicting the imbibition recovery factor of shale reservoirs that can consider the exchange rate of oil and water. Summary of the Invention

[0004] To solve the above problems, the present invention provides a method for predicting the imbibition recovery factor of shale oil reservoirs, including:

[0005] (1) Regress the oil-water relative permeability curves respectively to obtain the oil-water relative permeability;

[0006]

[0007] Taking the logarithm of Equation (37) gives:

[0008]

[0009] Obtaining the regression coefficients β w 、β o ;

[0010] ③ Using Equation (39) to calculate the average water saturation, initial water saturation, and residual oil saturation:

[0011]

[0012] Using the calculation results of Equation (39) and calculating the average oil-water relative permeability of the shale oil reservoir according to Equation (37);

[0013] In the formula:

[0014] k rw —— Water relative permeability, dimensionless;

[0015] k rw,,max —— Maximum water relative permeability, dimensionless;

[0016] k ro —— Oil relative permeability, dimensionless;

[0017] kro，max —— The maximum oil relative permeability, dimensionless;

[0018] S w —— Water saturation, %;

[0019] S wi —— Initial water saturation, %;

[0020] S or —— Residual oil saturation, %;

[0021] k ro —— Oil relative permeability, dimensionless

[0022] β w 、β o —— Regression coefficients;

[0023] (2) Obtain the fractional flow function F(S w ) and solve it:

[0024]

[0025] where F(S w ) is the fractional flow function under capillary-driven flow, and F(S w ) is given by Equation (49):

[0026]

[0027] where q w —— Water imbibition velocity of the aqueous phase, m / s;

[0028] S w,max —— Maximum water saturation, %;

[0029] (3) Obtain the reverse imbibition water saturation profile through Equation (68);

[0030]

[0031] In the formula: x —— The length corresponding to each water saturation in the core;

[0032] i —— The micro-element segment number in the core;

[0033] C —— Imbibition capacity characterization coefficient,

[0034] t —— Imbibition time, s;

[0035] φ —— Porosity;

[0036] (4) Establish the mass conservation equations for shale fractures and matrix at the field scale and predict the imbibition recovery factor;

[0037]

[0038]

[0039] Wherein: the subscripts f and m represent fracture and matrix, respectively;

[0040] q t ——Total flow velocity, m / s;

[0041] t—s;

[0042] T—Transfer function, 1 / s;

[0043] In order to obtain the transfer function expression of the imbibition dimensionless time t D The inverse imbibition oil recovery is defined as:

[0044]

[0045] Wherein: the subscript m represents matrix;

[0046] R—Predicted recovery;

[0047] R ∞ —Ultimate recovery;

[0048]

[0049] t D —Imbibition dimensionless time;

[0050] α is an empirical value obtained through experiments and is independent of wettability and rock properties;

[0051] Further, the imbibition recovery is predicted by using Equation (79):

[0052]

[0053] Advantages of the present invention:

[0054] (1) The present invention establishes a model for predicting the inverse imbibition recovery, in which the exchange rate of oil and water, i.e., the recovery, is described by a transfer function. Therefore, a prediction method for describing the relationship between the recovery of matrix blocks and time is provided, and the relationship between the recovery and time can be predicted only by basic fluid and rock physical property parameters.

[0055] (2) The present invention comprehensively considers the multi-scale pore size distribution of shale imbibition, the boundary slip of micro-nano pores, and the particularity of clay infiltration imbibition flow, and provides an effective boundary slip length calculation model that comprehensively considers the true slip between the aqueous phase and the pore wall surface and the apparent slip between the near-wall water and the bulk water.

[0056] (3) The present invention introduces the fractal theory and the oil-water two-phase seepage theory. Based on the B-L material balance equation, a reverse imbibition model for shale multi-scale pores is established, which can characterize the reverse imbibition ability and imbibition length of shale.

[0057] (4) The method for dividing the reservoir imbibition mode based on the relative size of the imbibition length and the fracture spacing in the present invention uses the reverse imbibition ability coefficient to develop a corresponding dimensionless time scale model, realizing the upgrade from the imbibition law of shale multi-scale pores to the prediction of the imbibition displacement recovery factor of shale reservoirs. Using the method of the present invention, the salinity of the fracturing fluid and the soaking time can be optimized according to the mineral composition and transformation degree of the reservoir to improve the recovery factor of shale oil wells. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0059] Figure 1 It is a schematic diagram of correcting the boundary condition by the slip length according to the present invention;

[0060] Figure 2 It is a schematic diagram of the post-fracture imbibition of a shale reservoir according to the present invention;

[0061] Figure 3 It is a schematic diagram of four boundary conditions for the reverse imbibition of a core according to the present invention;

[0062] Figure 4 It is a schematic diagram of the calculation result of the water saturation according to the present invention;

[0063] Figure 5 It is a schematic diagram of the calculation of the dimensionless capillary force function according to the present invention;

[0064] Figure 6 It is a schematic diagram of the calculation of the water saturation profile according to the present invention;

[0065] Figure 7 It is a schematic diagram of the change of the recovery factor with time according to the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0066] Combined with the description of the accompanying drawings and the specific embodiments of the present invention, the details of the present invention can be more clearly understood. However, the specific embodiments of the present invention described herein are only for the purpose of explaining the present invention and cannot be construed in any way as a limitation of the present invention. Under the teaching of the present invention, those skilled in the art can conceive any possible variations based on the present invention, and these should all be regarded as belonging to the scope of the present invention.

[0067] The present invention provides a method for predicting the spontaneous imbibition recovery rate of shale oil reservoirs, which can overcome the deficiencies of traditional methods that only consider capillary force and ignore boundary slip. Using the method provided in the present invention, the reverse imbibition length and imbibition capacity can be accurately determined, thereby providing a basis for predicting the field recovery rate and optimizing fracturing fluids. A method for predicting the spontaneous imbibition recovery rate of shale oil reservoirs includes the following steps:

[0068] First, obtain the physical property parameters of the shale reservoir, pore size distribution parameters, basic characteristic parameters of clay minerals, and fracturing fluid performance parameters.

[0069] The physical property parameters of the shale reservoir include mineral composition, wetting angle, porosity, and formation temperature; the pore size distribution parameters include the maximum pore size, minimum pore size, fractal dimension of pore size distribution, and fractal dimension of tortuosity; the basic characteristic parameters of clay minerals include density, cation exchange capacity, etc.; the fracturing fluid performance parameters include salinity, viscosity, and interfacial tension. Establish a calculation model for the microscopic imbibition force of shale, including: a calculation model for capillary force and gravity, and a calculation model for the osmotic pressure of clay minerals.

[0070] (1) Calculation model for capillary force and gravity

[0071] According to the direction of capillary imbibition, there is no gravity in the horizontal direction, and imbibition is only affected by capillary force, while in the vertical direction, it is affected by both capillary force and gravity. According to the reciprocal of the Bond number in Table 1 (i.e., ) the contribution of capillary force or gravity to dominate the shale imbibition can be predicted. The formula for the number is shown in Equation (1), where the numerator is the capillary force and the denominator is the gravity.

[0072]

[0073] In the formula:

[0074] σ—Interfacial tension between two phases, mN / m;

[0075] φ—Porosity, %;

[0076] K—Permeability, mD;

[0077] Δρ—Density difference between two phases, g / cm 3 ;

[0078] H—Height of porous medium, cm;

[0079] C—Geometric size coefficient (≈0.4 for circular pipe).

[0080] Table 1 Classification intervals of the reciprocal of Bond number

[0081]

[0082] The value is affected by pore size: The main pore size distribution interval of conventional sandstone is >2μm; that of tight sandstone is 30nm - 2μm; and that of shale is 0.7nm - 0.1μm. Conventional sandstone The values are all greater than 1, and those of shale are significantly greater than 5, that is, the dominant self-imbibition force of shale is capillary force and gravity can be ignored.

[0083] The capillary force can be calculated by the L - Y equation, as shown in Equation (2).

[0084]

[0085] In the formula:

[0086] P c —Capillary force, MPa;

[0087] σ—Oil - water interfacial tension, N / m;

[0088] θ—Wetting angle;

[0089] d—Pore diameter, m.

[0090] (2) Calculation model of osmotic pressure of clay minerals

[0091] Based on the mathematical relationship of thermodynamics, the calculation formula of osmotic pressure is as shown in Equation (3).

[0092]

[0093] In the formula:

[0094] P π —Osmotic pressure, KPa;

[0095] a f 、a i —Activity of external fracturing fluid and original formation water respectively, dimensionless, and the activity of fresh water is 1;

[0096] R—Gas constant, 0.08314 (L·KPa) / (mol·K);

[0097] T—Formation temperature, K;

[0098] V w— The molar volume of water, with a value of 0.018 L / mol.

[0099] E π — The semipermeable membrane efficiency, and the calculation formula is shown in Equation (4).

[0100]

[0101] Among them, K s is the distribution coefficient of the solute in the semipermeable membrane:

[0102] K s = C a / C s (5)

[0103] In the formula:

[0104] C s — The arithmetic mean of the solute in the solutions on both sides of the semipermeable membrane, mol / L;

[0105] C a — The anion concentration in the pores of the semipermeable membrane, mol / L.

[0106]

[0107] In the formula:

[0108] C sh — The molar concentration of the solute in the original formation water of the shale formation, mol / L;

[0109] C f — The molar concentration of the solute in the external working fluid (such as fracturing fluid) in the fracture, mol / L.

[0110] Among them, C a can be calculated by the following formula:

[0111]

[0112] In the formula:

[0113] E CEC — The cation exchange capacity of clay minerals, mmol / 100 g;

[0114] ρ clay — The density of clay minerals, g / cm 3 , and the clay mineral composition can be obtained through whole-rock analysis and calculated by weighted average;

[0115] φ C — The porosity of the semipermeable membrane in shale clay minerals, that is, the porosity of clay minerals, %.

[0116] φ c = Vc × ρ clay (8)

[0117] In the formula:

[0118] V c — Pore volume of clay minerals in shale, cm 3 / g, obtained by calculating the multiple pore splitting in Section 2.6.

[0119] C in formula (4) c is the cation concentration in the semipermeable membrane and is calculated by the following formula:

[0120] C c = C a + E CEC ρ clay (1 - φ c ) (9)

[0121] R in formula (4) is the friction coefficient ratio, and the specific definition is:

[0122]

[0123] In the formula:

[0124] f — represents the friction coefficient, and c, a, w, m represent cations, anions, water, and semipermeable membrane. f ij represents the frictional resistance between 1 mole of component i and an infinite amount of component j.

[0125] The friction coefficient f of cation - water can be obtained by referring to the literature cw , and the friction coefficient f of the anion - water system aw ; R ca-m can be obtained from the ratio of the hydrodynamic radii of cations and anions (for NaCl, R ca-m = 1.8); f am can be obtained by diffusion coefficient testing. For an abnormally loose semipermeable membrane, f am << f aw , and for an abnormally dense semipermeable membrane, the friction coefficient between the anion and the semipermeable membrane is very large, f am ≥ f aw .

[0126]

[0127] It can be seen from formula (11) that when a f < a i , the osmotic pressure is the driving force, and the external working fluid flows into the shale pores; when a f = a i , the osmotic pressure is 0, and the fluid neither flows in nor out; when a f < a iWhen the osmotic pressure is the resistance, the original water in the shale flows out.

[0128] In practical applications, the salinity difference is often directly calculated. As shown in Equation (12), when ΔC is less than 1 mol / L, the calculation error is less than 5%:

[0129] P π = εE π RTΔC = εE π RT(C sh -C f )(12)

[0130] In the formula:

[0131] ε—the number of ions after the solute ionization (e.g., the value of NaCl is 2), dimensionless.

[0132] Establish a reverse imbibition flow model for shale multi-scale pores, including: establishing a calculation model for the matrix apparent permeability considering the boundary slip effect of micro-nano pores; based on the Buckley-Leverct material balance equation and using the two-phase oil-water seepage theory, establishing a reverse imbibition flow model for multi-scale pores dominated by capillary force and osmotic pressure; characterizing the reverse imbibition capacity and imbibition length of shale multi-scale pores.

[0133] (1) Characterization of shale matrix permeability

[0134] First, based on the single nano-pore fluid flow model with effective slip effect and the fractal basic theory, establish a shale reservoir matrix fluid model. The total flow rate Q of all nano-capillaries within the cross-sectional area of the shale unit can be obtained by summing the flow rates of all capillaries, that is, Equation (13):

[0135]

[0136] On the macroscopic scale, the boundary slip effect is not likely to occur, and even if it occurs, it can be ignored. At this time, it is treated as a model without boundary and without slip; when entering the nano-microscopic scale, a slip velocity is generated at the boundary, and the boundary slip length of the pore wall surface needs to be corrected, as shown in the attached instruction Figure 1 As shown. For the steady laminar incompressible viscous fluid in a circular tube, by correcting the flow boundary conditions, the continuous fluid mechanics is still effectively applicable to describe the fluid flow in micro-nano pores.

[0137] Assuming the fluid is incompressible and laminar, according to the Navier-Stokes equation, the mathematical model of the fluid flow in a single circular tube can be obtained, as shown in Equation (14)

[0138]

[0139] The velocity at the center of the circular tube remains unchanged, and the pore wall surface has a slip flow velocity. The boundary conditions are as shown in Equation (15).

[0140]

[0141] In the formula:

[0142] L s ——Solid-liquid boundary fluid slip length, m.

[0143] Substitute Equation (15) into Equation (14), and the velocity distribution equation is obtained as shown in Equation (16).

[0144]

[0145] Integrate Equation (16) to obtain the flow rate considering the boundary slip effect.

[0146]

[0147] The viscosity of boundary water is quite different from that of bulk water, and obvious slip will occur at the interface between boundary water and bulk water. Therefore, in practical applications, the effective slip length should be used to replace the restricted fluid slip length considering the true slip effect and apparent slip effect.

[0148]

[0149] In the formula:

[0150] L se ——Boundary effective slip length, m;

[0151] L sa ——Apparent slip length, m;

[0152] μ ∞ ——Bulk water viscosity, Pa·s;

[0153] d——Micro-nano capillary diameter, m;

[0154] μ d ——Boundary water effective viscosity, Pa·s.

[0155] It can be seen from Equation (18) that the effective boundary slip length depends on the pore surface wettability, fluid effective viscosity and pore size.

[0156] When the water phase flow in micro-nano scale capillaries increases the order of its nano-scale pores, the fluid viscosity can no longer accurately describe the fluid viscosity near the pore wall. The effective viscosity is determined by the viscosity of the fluid flowing through the core and the interface area, and can be written as Equation (19).

[0157]

[0158]

[0159]

[0160] In the formula:

[0161] μ d —— Effective viscosity of fluid in micro-nano scale capillary, Pa·s;

[0162] μ i —— Viscosity of boundary fluid, Pa·s;

[0163] A id —— Area occupied by boundary water, m 2 ;

[0164] d c —— Critical thickness of boundary water, m, and the fitting value of experimental test is 0.7 nm;

[0165] μ ∞ —— Viscosity of bulk water, Pa·s;

[0166] A td —— Cross-sectional area of capillary, m 2 .

[0167] The fractal dimension D f is used to characterize the roughness of the porous medium. The larger D f , the more complex the pore structure is.

[0168]

[0169] In the formula:

[0170] d—— Diameter of micro-nano pore, nm;

[0171] d min ,d max —— Minimum and maximum pore diameters of shale, nm;

[0172] D f —— Pore size fractal dimension, 0 < D f < 2.

[0173] The negative sign in formula (22) indicates that the number of capillaries decreases with the increase of pore diameter, and -dN > 0.

[0174] Substituting equations (17), (18) to (22) into formula (13), the total flow rate of fluid in the unit cross-sectional area can be obtained:

[0175]

[0176] Formula (23) is the fractal calculation model for flow transmission of the fluid in the shale reservoir to be analyzed, and the effective viscosity is μ dis a variable related to the pipe diameter of the shale reservoir. Therefore, it is difficult to further simplify Equation (23) by integration. For the convenience of solution, the shale micro-nano pore size distribution is discretized into J tiny units, and the flow rate Q within each unit (d min,i ≤d i ≤d max,i ) can be written as Equation (24): i It can be written as Equation (24):

[0177]

[0178] In the formula:

[0179]

[0180]

[0181]

[0182]

[0183]

[0184] Then, by algebraically summing the flow rates of each tiny unit in Equation (24), the total volume flow rate expression can be obtained:

[0185]

[0186] In the formula:

[0187] Q f ——The total volume flow rate, nm 3 / s;

[0188] J——The number of tiny units discretized from the nano-pore diameter distribution, unit.

[0189] The total flow rate expression (fluid flow model) is obtained:

[0190]

[0191]

[0192] In the formula:

[0193] d max,i ——The maximum diameter of the nano-pore in the i-th section, nm;

[0194] d min,i ——The minimum diameter of the nano-pore in the i-th section, nm;

[0195] μ d,i ——The effective viscosity of the fluid, mPa·s.

[0196] Calculate the permeability of the shale reservoir to be analyzed based on the hydrodynamic model and Darcy's law. According to the generalized Darcy's law, the present invention can obtain the flow equation of the porous medium of the shale reservoir:

[0197]

[0198] The apparent permeability of shale is:

[0199]

[0200] In the formula:

[0201] k——Apparent permeability of shale, μD.

[0202] Combining Equation (30) and Equation (32), the shale permeability can be obtained as:

[0203]

[0204] Counter-current imbibition flow model

[0205] First, consider the one-dimensional counter-current flow of two immiscible and incompressible phases through a homogeneous porous medium. In the counter-current flow, water invades a closed system and oil leaves in the opposite direction. The following assumptions are made:

[0206] ① The system is homogeneous;

[0207] ② The fluid is incompressible;

[0208] ③ The inlet capillary pressure is zero and there is no capillary back pressure;

[0209] ④ It conforms to the traditional multi-phase Darcy's law;

[0210] ⑤ The compressibility of the rock is not considered;

[0211] ⑥ Before the imbibition front reaches the far boundary of the sample, the solution is a function of the parameter ;

[0212] ⑦ The oil phase and the water phase are immiscible during the imbibition process.

[0213] The one-dimensional mass conservation equation of the incompressible fluid can be expressed as

[0214]

[0215] In the formula:

[0216] S w ——Water saturation, dimensionless;

[0217] q w ——Imbibition velocity of the water phase, m / s;

[0218] x—the self - absorption length, m;

[0219] t—the self - absorption time, s.

[0220] φ—the porosity

[0221] The Darcy velocity equation of the wetting phase is shown in Equation (35):

[0222]

[0223] In the formula:

[0224] λ w —the mobility of the wetting phase, 1 / Pa·s, λ w =k rw / μ w ;

[0225] λ nw —the mobility of the non - wetting phase, 1 / Pa·s, λ o =k ro / μ o ;

[0226] λ t —the total mobility, λ t =λ o +λ w , 1 / Pa·s;

[0227] q t —the total flow velocity, m / s. q t =q o +q w For reverse self - absorption, water invades a closed system, so the direction in which the non - wetting phase leaves is opposite to that of water. In this case, q nw =-q w , that is, q t =0:

[0228] P c —the capillary pressure, defined as the difference between the non - wetting phase pressure and the wetting phase pressure, MPa;

[0229] P π —the self - absorption osmotic pressure of clay pores hydration, MPa;

[0230] ρ w 、ρ o —the densities of the wetting phase and non - wetting phase fluids, g / cm 3 ;

[0231] g x —the acceleration of gravity, m / s.

[0232] Neglecting gravity, we get

[0233]

[0234] Fitting relative permeability using formulas:

[0235] ① Select representative oil-water relative permeability curves;

[0236] ② Regress the oil-water relative permeability curves respectively using the following formulas:

[0237]

[0238] Taking the logarithm of Equation (37) gives:

[0239]

[0240] The regression coefficients β w and β o can be obtained using Equation (38).

[0241] ③ Calculate the average S w , S or , and S wi using the following formula:

[0242]

[0243] Using the calculation results of Equation (39), the average oil-water relative permeability of the shale reservoir can be obtained according to Equation (37).

[0244] In the formula:

[0245] k rw —— Water relative permeability, dimensionless;

[0246] S w —— Water saturation, %;

[0247] S wi —— Initial water saturation, %;

[0248] S or —— Residual oil saturation, %;

[0249] k ro —— Oil relative permeability, dimensionless.

[0250] The Buckley-Leverett fractional flow is defined as: f w = λ w / λ t , and Equation (36) becomes:

[0251]

[0252] Therefore, for reverse imbibition, the wetting phase velocity equation is as shown in Equation (41).

[0253]

[0254] Equation (41) becomes:

[0255]

[0256] Substituting (21) into Equation (42), Equation (42) is rearranged as:

[0257]

[0258] In the traditional Buckley-Leverett analysis, if capillary forces are ignored, the solution to this flow equation is a function of v = x / t, where v is the dimensionless wave speed. However, in capillary-controlled flow, it has been proposed by previous researchers that the distance x representing the imbibition front is expressed in terms of √t rather than t. Based on this, here the imbibition distance x is used to represent the imbibition front, √t is used in place of t, and x is proportional to √t. Therefore, we define the scaling factor ω as:

[0259]

[0260] S w = S w (ω)(45)

[0261] Since the imbibition distance is proportional to √t, the imbibition volume must also be proportional to √t. Therefore, there exists a proportionality constant C (with units of m / √t) that expresses the aqueous-phase imbibition velocity q w as a function of 1 / √t.

[0262]

[0263]

[0264] Where:

[0265] C — Imbibition capacity characterization coefficient,

[0266] L(t) — Imbibition length of the wetting phase, m;

[0267] C characterizes the rate of shale imbibition velocity and is a parameter that describes the intrinsic properties of the rock / fluid. The value of C depends on shale and fluid properties such as wettability, fluid viscosity, relative permeability, and absolute permeability. If one knows how to calculate C, there is a simple method to quantify how shale and fluid properties affect the imbibition capacity of shale in the reverse imbibition mode.

[0268] The key new step in finding the general solution is to express the solution in terms of the derivative of the fractional flow function F:

[0269]

[0270] F is the fractional flow function in the case of capillary-driven flow, and we assume that F is given by the ratio of the water saturation to the maximum value of the function:

[0271]

[0272] Differentiate Equation (48):

[0273]

[0274] To analyze and calculate the cumulative water absorbed into the porous medium for each water saturation (S wi to S S range), the diffusion equation in linear form will be used to redefine the derivative of the water saturation: w0 )

[0275]

[0276]

[0277] Substitute Equations (51) and (52) into Equation (43) and rearrange to obtain:

[0278]

[0279] To eliminate the second derivative, integrate Equation (53).

[0280]

[0281] Then, according to the Buckley-Leverett analysis, for some capillary fractional flow F (1 ≥ F ≥ 0) and constant C, ω can also be defined as:

[0282]

[0283] Here, introduce the factor 2C / φ, where C is a constant to make F dimensionless. Differentiate ω with respect to S w to give:

[0284]

[0285] Substitute Equations (55) and (56) into Equation (54) to obtain:

[0286]

[0287] Rearrange to obtain the second-order ordinary differential equation for F as shown in Equation (58).

[0288]

[0289] Equation (58) is an implicit non-linear second-order ordinary differential equation. Different from the complete numerical solution of the partial differential equation that controls in space and time, equation (58) can be simply solved. The total amount of fluid entering the system is no longer the boundary condition, but is now embedded in the constant C. Equation (58) is discretized in the saturation space, which is a very fast and accurate method of solution.

[0290] From this, F and C can be determined. Once F′ and C are obtained, the water saturation profile is obtained. F, F′ and C can be determined based on the simple backward difference numerical solution of the equation.

[0291] Equation (58) is only valid before the saturation front reaches the boundary. The failure time t * can be obtained by setting the imbibition distance x(S w , t*) = L in a rock sample with length L.

[0292]

[0293] In the formula:

[0294] t * ——Effective imbibition time, s. When t ≤ t * , the backward difference solution method is effective.

[0295] The solution method for obtaining F and C based on equation (58) is introduced below. The backward difference gives:

[0296]

[0297]

[0298] Let F(S w ) = X. Substituting equations (60) and (61) into equation (58) gives:

[0299]

[0300] Equation (62) is a quadratic equation of ax 2 + bx + c = 0. Using the solution of the quadratic equation, the solution of the counter-phase imbibition equation is obtained:

[0301]

[0302] When obtaining the starting values F(S w -ΔS w ) and F(S w -2ΔS w ), the approximate solution of F(S w,strat ) can be found. Take F(S w )w,max ) = 1 and F(S w,max -ΔS w ) are used as two initial values. Among them, it is assumed in the text that the maximum saturation is not transported into the medium, so F′(S w,max ) = 0.

[0303] F(S w,max -ΔS w ) = F(S w,max ) - ΔS w F′(S w,max ) = 1 (64)

[0304] Through an iterative process, the concept of backward difference approximation is used to find the unknown constant C. By choosing an initial value for C and then continuously changing it until certain convergence conditions of F are satisfied. By reviewing some physical properties of the solution, two equivalent convergence criteria can be obtained.

[0305] The first criterion: F is a fractional flow function, from which Equation (65) is obtained.

[0306] F(S wir ) = 0 (65)

[0307] The second criterion: Due to material balance, the saturation curve must be equal to the total pore absorption volume.

[0308]

[0309]

[0310] All in all, the solution process of the reverse imbibition analytical solution is as follows:

[0311] ① Determine F” through backward difference approximation;

[0312] ② Iteratively determine a finite number of n water saturation point calculation nodes;

[0313] ③ Iterate the constant C until F(S w ) converges to the correct solution.

[0314] During the reverse imbibition process in a one-dimensional horizontal medium, each saturation moves over time, and the imbibition distance is shown in Equation (68):

[0315]

[0316] In the formula: x—the length corresponding to each water saturation in the core, m.

[0317] A prediction model for the imbibition oil recovery rate of shale reservoirs is established, including: dividing the imbibition mode of the reservoir based on the relative size of the fracture spacing and the imbibition length; developing a corresponding dimensionless time scale model using the reverse imbibition capacity coefficient; establishing a mathematical relationship between the recovery rate and the dimensionless imbibition time based on the fracture-matrix transfer function, so as to realize the upgrade from the multi-scale pore imbibition data of shale cores to the prediction of reservoir recovery rate.

[0318] (1) Division of imbibition flow modes

[0319] After large-scale hydraulic fracturing of shale reservoirs, the matrix is surrounded by fractures and saturated with fracturing fluid. It is generally believed that shale reservoirs are hydrophilic, and the fracturing fluid imbibes into the matrix pores through the fracture walls. The imbibition occurs in two ways: forward and reverse, and the imbibition characteristics are as follows:

[0320] ① During the forward imbibition process, the fracturing fluid is inhaled from the core end face under the drive of capillary force and osmotic pressure. The shale oil in the pore throat is only affected by the imbibition force at one end and is displaced to the other end face of the core for discharge.

[0321] ② The reverse imbibition process is more complex and includes two oil displacement forms: ① The fracturing fluid first imbibes along the pore throat surface and tries to occupy the pore corners, and the shale oil is displaced to the center of the pore throat and then displaced out, with the oil and water flowing in opposite directions; ② Due to the heterogeneity of the pore size distribution, the fracturing fluid preferentially imbibes into the fine pore throats with large capillary force, and the shale oil is discharged from the thick pore throats on the same side of the core.

[0322] The imbibition mode is different for different reservoir transformation degrees. The water content, porosity and permeability of shale reservoirs are extremely low. Relative to fractures, they are impermeable boundaries. The Figure 2 description in the specification appendix

[0323] As Figure 2 In the No. ① area in (a), the distance between two fractures is far, and the matrix length is much higher than the reverse imbibition distance, which is equivalent to Figure 2 the one-end open (OEO) boundary condition shown in (b①). At this time, only the reverse imbibition occurs.

[0324] Figure 2 In the No. ② area in (a), the distance between two fractures is close, and the reverse imbibition occurs simultaneously from both end faces of the core matrix, which is Figure 2 the two-end open (TEO) boundary condition shown in (b②). At this time, there is also the forward imbibition of the fracturing fluid to displace the shale oil, that is, the forward imbibition, but the reverse imbibition still dominates.

[0325] Figure 2 In the No. ③ area in (a), the matrix is surrounded by fractures, and the distance between the surrounding fractures is less than the reverse imbibition distance, which is Figure 2(Boundary condition of the open - around (AFO) as shown in (b③), where reverse imbibition occurs on all four side - end faces of the core.)

[0326] Figure 2 (As shown in (b④), the boundary condition is that the sides of the core are open and the two end - faces are of matrix (TWC). At this time, reverse imbibition occurs on both sides.)

[0327] Generally speaking, reverse imbibition always dominates in fractured shale reservoirs. The micro - structure of shale reservoirs is highly heterogeneous. During the imbibition - displacement process of oil, it is very easy to form bypass flow. Shale oil can only be effectively exploited by relying on reverse imbibition. In terms of the degree of macro - reservoir transformation: when the fracture spacing is less than 2 times the reverse - imbibition action distance, reverse imbibition will occur simultaneously on both end - faces of the core; otherwise, reverse imbibition will only occur on one side; when the fracture spacing is less than the forward - imbibition action distance, forward imbibition will occur concomitantly.)

[0328] (2) Reverse dimensionless time - scale model

[0329] When a water - saturated core invades water, the volume of the displaced oil changes with time. In order to describe the influence of key parameters such as interfacial tension, fluid viscosity, porosity, permeability, and core size on the experimental results, many scholars have developed the dimensionless time - scale theory of water - oil reverse imbibition.)

[0330] Mattax et al. believed that when the volume of the liquid - phase imbibition - displacement of the oil - phase is constant, the imbibition time is proportional to the square - root of the imbibition distance, and the ratio is a constant value. They first introduced the concept of dimensionless imbibition time and used micro - pore imbibition to predict the recovery rate of shale matrix.)

[0331]

[0332] For cores and all reservoir matrices with the same physical properties and geometric shapes, the same recovery curve will be obtained. This means that the imbibition test data of the core can measure the imbibition of all reservoir matrix blocks with the same shape and rock type, and this process will be through the dimensionless time - scale of imbibition.)

[0333]

[0334] In the formula:

[0335] α——Unit conversion factor, 3.16×10 -4 .

[0336] In order to relate the dimensionless time - scale to the physical properties of the oil - water two - phase and the rock during the reverse imbibition process, the imbibition length of the wetting phase is used for imbibition scaling.)

[0337]

[0338] In the formula:

[0339] t D —— Dimensionless time, s.

[0340] Further explanation of the dimensionless time:

[0341]

[0342] Substituting Equation (72) into Equation (71), the dimensionless time expression can be obtained.

[0343]

[0344] In the formula:

[0345] t c —— Characteristic time, dimensionless.

[0346] Next, a more in-depth analysis of the characteristic length will be carried out. First, the characteristic length is normalized, where

[0347]

[0348] In the formula:

[0349] L c —— Self-aspiration characteristic length, m;

[0350] V b —— Core volume, m 3 ;

[0351] A i —— Self-aspiration area in the i-th direction, m 2 ;

[0352] l Ai —— Distance from the self-aspiration front to the non-flow boundary from the open surface, m.

[0353] The self-aspiration mechanism under different core boundary conditions is involved here. Reflected in the two-dimensional plane, there are mainly 4 types of self-aspiration boundaries in the real reservoir, as Figure 3 shown.

[0354] ① All faces open (AFO): The fracturing fluid flows into and out of the core from all directions.

[0355] ② One end open (OEO): The fracturing fluid flows horizontally into the core from one side, and the other side of the core is sealed.

[0356] ③ Two ends open (TEO): The fracturing fluid flows from one end to the other in the central axis direction of the core, while the top and bottom are sealed.

[0357] ④ Two ends closed (TEC): The fracturing fluid flows vertically from top to bottom, while the other side is completely isolated.

[0358] Figure 3 Four boundary conditions for the reverse imbibition of cores are described as follows: (a) is one - end open (OEO), (b) is two - ends open (TEO), (c) is two - ends closed (TEC), and (d) is all - faces open (AFO).

[0359] Among these four boundary conditions, the TEO boundary will undergo forward imbibition; the AFO, OEO, and TEC boundaries are mainly characterized by reverse imbibition. By introducing the imbibition characteristic length parameter, the influence of core shape and size on the imbibition law can be fully considered, and the imbibition direction at the core scale can be made consistent with that of the fracture - matrix system.

[0360] (3) Field - scale recovery model

[0361] The migration of oil in the shale matrix pores to the fractures depends on imbibition. The fracturing fluid enters the pores under the action of capillary force and osmotic pressure to displace the oil. The fluid transport between the matrix and the fractures is modeled by the transfer function T, upgrading the imbibition oil - displacement experimental data at the core scale to the field scale.

[0362] In the fracture - pore dual model, the mass conservation equations for the fractures and the matrix are as follows:

[0363]

[0364]

[0365] In the formula:

[0366] T - transfer function, 1 / s.

[0367] To obtain the expression of the transfer function with respect to the dimensionless imbibition time t D The reverse imbibition oil - displacement recovery factor is defined as:

[0368]

[0369] In the formula:

[0370] R ∞ — Ultimate recovery factor, taken as 0.9 in this paper.

[0371]

[0372] Furthermore, a prediction model for the imbibition recovery factor of the fracture - matrix system is established:

[0373]

[0374] In the formula:

[0375] α ≈ 70. This empirical value is obtained through experiments and is independent of wettability and rock properties.

[0376] Predicting recovery factor based on imbibition capacity

[0377] According to the VG parameters and substituting them into the relative permeability calculation model in Equation (37), the relative permeability curve can be obtained, as Figure 4 shown; the fractional function F(S w ) can be iteratively solved by Equation (48). Subsequently, based on Equation (68), the reverse imbibition water saturation profile can be obtained, as Figure 6 shown. The above process is the solution steps of the reverse imbibition flow equation. Since the imbibition distance - water saturation profiles at different times are different, in order to obtain a saturation curve with a unified profile for the same core, the experimental data is scaled by the square root of time, with the abscissa ω = x / √t.

[0378] Table 2 Basic data

[0379]

[0380] From Figures 4 to 7 the calculation results shown, it can be seen that as the imbibition time increases, the water absorption profile advances, the capillary force decreases accordingly, the imbibition driving force gradually decreases, and the imbibition oil displacement rate first increases and then decreases. This is because the initial water saturation in the core is relatively low, and the fracturing fluid quickly occupies the small pore volume under the action of capillary force; as the water saturation in the core matrix increases, the seepage resistance increases, and the imbibition oil displacement efficiency gradually decreases.

[0381] Although the disclosed embodiments of the present invention are as above, the content described is only an embodiment adopted for the convenience of understanding the present invention and is not intended to limit the present invention. Any person skilled in the art within the technical field to which the present invention pertains can make any modifications and changes in the form of implementation and details without departing from the spirit and scope disclosed by the present invention. However, the scope of patent protection of the present invention shall still be subject to the scope defined by the appended claims.

Claims

1. A method for predicting the imbibition recovery rate of shale oil reservoirs, comprising: (1) Respectively regress the oil-water relative permeability curves to obtain the oil-water relative permeability; Taking the logarithm of Equation (37) gives: Obtain the regression coefficient β using Equation (38) w , β o ; Use Equation (39) to calculate the average water saturation, initial water saturation, and residual oil saturation: Using the calculation results of Equation (39), calculate the average oil-water relative permeability of the shale oil reservoir according to Equation (37); In the formula: k rw —— water relative permeability, dimensionless; k rw,max —— maximum relative water permeability, dimensionless; k ro —— Oil relative permeability, dimensionless; k ro，max —— The maximum relative oil permeability, dimensionless; S w —— water saturation, %; S wi —— Initial water saturation, %; S or —— Residual oil saturation, %; k ro —— Relative permeability of oil, dimensionless β w and β o —— regression coefficients; (2) Obtain the fractional flow function F(S w ) and solve it: where F(S w ) is the fractional flow function under capillary-driven flow, and F(S w ) is given by Equation (49): where qw is the imbibition velocity of the water phase, m / s; S w,max —— maximum water saturation, %; (3) Obtain the reverse imbibition water saturation profile through Equation (68); In the formula: x is the length corresponding to each water saturation in the core; i is the micro-element section number in the core; C——Coefficient for characterizing self-priming ability, t is the imbibition time, s; φ is the porosity; (4) Establish the mass conservation equations for shale fractures and matrix at the field scale and predict the imbibition recovery rate; In the formula: the subscripts f and m respectively represent fracture and matrix; qt is the total flow rate, m / s; t—s; T is the transfer function, 1 / s; In order to obtain the transfer function expression for the dimensionless time t of self-aspiration D the reverse self-aspiration oil displacement recovery factor is defined as: In the formula: the subscript m represents matrix; R is the predicted recovery rate; R ∞ — Ultimate recovery factor; t D — Self-priming dimensionless time; α is an empirical value obtained through experiments and is independent of wettability and rock properties; Further predict the imbibition recovery rate using Equation (79):

Images