As shown in Figure 1, there is a water injection well in a rectangular closed formation, which is injecting water into the reservoir with high pressure and large volume. The entire physical model is divided into two areas. The inner area is a region with higher seepage capacity. It is a reformation area, formed due to water injection that exceeds the formation fracture pressure during high-pressure water injection. Concurrently, microseismic monitoring data reveal the reformation area to be elliptical. The outer region refers to the initial formation. In order to describe the two-phase fluid flow during high-pressure water injection in low permeability reservoirs, some idealizations and assumptions are made as follows.

(1) the fluid flow within the reservoir consists of oil-water two-phase flow under isothermal conditions;

This study used two independent numerical modules to calculate the saturation distribution and pressure distribution during high-pressure water injection of low permeability reservoirs (Zhang et al., 2022). Based on the finite element method, the saturation changes of the oil phase and the water phase were calculated through the phase transport in porous media module. By improving Darcy’s law module (Liu et al., 2019), the pressure field of the oil-water mixture was calculated taking into account the threshold pressure gradient. Then the calculated pressure field and saturation field were coupled to obtain a complete mathematical model.

(1) Phase Transport in Porous Media

Based on the law of conservation of mass, the continuity equation of the fluid in the reservoir can be expressed as: ∂ ρ l ϕ s l ∂ t + ∇ ⋅ ρ l v l = ρ l q l l = o , w (1) where ρ is the density, t is the time, v is the velocity vector, s is the saturation, ϕ is the porosity, q refers to the sink or source term in the reservoir, the subscripts o and w represent oil phase and water phase.

The kinematic equations of the oil and water phases can be obtained under the condition of considering the capillary forces and they can be expressed as Eqs 2, 3: v o = − k k r o μ o ∇ p o (2) v w = − k k r w μ w ∇ p w (3) where v o is the oil phase velocity, v w is the water phase velocity, k is the permeability of porous media, k r o is the oil relative permeability, k r w is the water relative permeability, μ o is the viscosity of the oil phase, μ w is the viscosity of the water phase, p w is the pressure of the water phase, p o is the pressure of the oil phase, ∇ represents the Hamiltonian operator. The oil phase pressure and water phase pressure satisfy Eq. 4: p o = p w + p c (4) where p c is the capillary pressure, it can be calculated from Eq. 5 (Blunt, 2017): p c = p e n s w − s w i 1 − s w i − s o r λ (5) where p e n is the entry capillary pressure, s w is the water saturation, s w i refers to the irreducible water saturation, s o r refers to the residual oil saturation, λ represents the capillary pressure exponent.

The control equations can be obtained, by substituting Eqs. 2, 3 into Eq. 1, − ρ l k k r l μ l ∇ 2 p l = s l ρ l ϕ C ∂ p l ∂ t + ρ l ϕ ∂ s l ∂ t   l = o , w (6) where p is the pressure, C is the compressibility.

Besides, the saturation of the oil and water phases satisfies Eq. 7, s o + s w = 1 (7)

(2) Oil and water phase pressure calculation method

To calculate the pressure distribution of the oil-water mixture, the properties of the mixture must first be obtained. Chen et al. (Chen et al., 2006) proposed the calculation method of mixture density and viscosity, which can be expressed as: ρ ¯ = ∑ l ρ l s l (8) μ ¯ = ρ ¯ ∑ l k r l ρ l μ l (9)

Darcy’s law for the oil-water mixture can be determined as: v ¯ = − k μ ¯ ∇ p ¯ (10) where v ¯ refers to the velocity of the oil-water mixture, p ¯ represents the pressure for the oil-water mixture. Meanwhile, the injection phase is different from the development phase. In the physical simulation process of reservoir development, the reservoir energy is gradually depleted and the pore pressure becomes smaller as crude oil is extracted. The reservoir seepage capacity is reduced and its macro performance is reduced permeability. However, in the injection stage, the continuous increase of injected fluid effectively replenishes the reservoir energy and improves the seepage capacity of the reservoir. The macroscopic manifestation is the decrease of permeability. Therefore, we define a macroscopic stress sensitivity coefficient to characterize permeability changes. The greater the stress sensitivity coefficient, the greater the variation range of permeability. According to the method proposed by Pedrosa et al. (Pedrosa, 1986; Zhu et al., 2021), the relationship between permeability, stress sensitivity coefficient and reservoir pressure can be expressed as: k = k i e − α p i − p (11) where α is the stress sensitivity factor. The compressibility of rock and fluid is considered by Eqs. 12, 13: ρ l = ρ i 1 + C l p − p i (12) ϕ = ϕ i 1 + C t p − p i (13) where C l is the fluid compressibility, C t is the rock compressibility. The continuity equation can be expressed as, ∇ ⋅ ρ ¯ v ¯ + ∂ ρ ¯ ϕ ∂ t = 0 (14)

The discrete fracture model is used to describe the flow of fluid in fractures. The kinematic equation of fluid flow in fractures can be expressed as: v ¯ f = − k f μ ¯ ∇ F p ¯ f (15)

Where v ¯ f refers to the velocity for the oil-water mixture of the fracture, ∇ F p ¯ f refers to the pressure gradient tangent to the fracture surface, and k f is the permeability of the fracture. Therefore, the governing equation of the fracture can be expressed as ∇ F ⋅ d f ρ ¯ v ¯ f + d f ∂ ρ ¯ ϕ f ∂ t = 0 (16)

Where ϕ f refers to the porosity of the fracture, and d f refers to the width of the fracture. The average pressure of the oil-water mixture can be calculated by Eqs. 14, 16, and then the pressure of the oil phase and water phase can be calculated by Eqs. 4, 5, 17 (Chen et al., 2006): p ¯ = s o p o + s w p w (17)

At the initial moment, the reservoir pressure and water saturation are equal everywhere. The initial condition can be expressed as: p = p i , t = 0 (18) s w = s w i , t = 0 (19)

The outer boundary of the model ∂ Ω 1 is assumed to be closed, it can be written as follows: − n ⋅ ρ v ¯ ∂ Ω 1 = 0 (20)

Constant pressure water injection was used in the simulation. The inner boundary conditions can be expressed as: p = p i n j (21)

Solution Procedure for Evaluation of High Pressure Injection Capacity in Low Permeability Reservoirs Shown in Figure 2.

This study simulates high-pressure water injection using COMSOL Multiphysics software. In the simulation, the phase transport in porous media module and the Darcy’s law module were selected. The two modules achieve coupled solutions through the multiphysics interface multiphase flow in porous media. As shown in Figure 3, the geometric model is a square reservoir of 800 m × 800 m, and the inner area is an elliptical reformation area with a long semi-axis of 150 m and a short semi-axis of 60 m. A water injection well with constant pressure is located at the center of the reservoir. Due to the symmetry of the physical model, in order to reduce the calculation amount of the model and increase the calculation speed, only half of the reservoir area is calculated. The meshing method selects unstructured mesh, as shown in Figure 3. In order to accurately describe the fluid flow near the wellbore and fracture, mesh refinement was performed around the wellbore and fracture. The numerical model contains a total of 1,375 domain cells and 175 boundary elements. Other parameter settings in the numerical model are shown in Table 1.

For low permeability reservoirs, the effective treatment of threshold pressure gradient is one of the key issues in high-pressure water injection simulation. The kinematic equation for considering the influence of gravity in the Darcy’s Law module of COMSOL Multiphysics can be described as: v l , β = − k l , β μ l ∂ p l ∂ β − ρ ⋅ g β   β = x , y (22) where g β is the gravitational acceleration in the reservoir, the subscripts x and y represent the direction. Following the method proposed by Liu et al. (Liu et al., 2019), we utilize piecewise functions to replace the gravitational acceleration in Darcy’s equation of motion to simulate the threshold pressure gradient. This replacement can be expressed as follows: g β = λ β ρ ∂ p ∂ β > λ β − λ β ρ ∂ p ∂ β < − λ β ∂ p ∂ β 1 ρ − λ β ≤ ∂ p ∂ β ≤ λ β (23)

The piecework function of Eq. 23 can be written into COMSOL by programming to realize the numerical processing of threshold pressure gradient.

This study focuses on the development of low-permeability oil fields in the Bohai Bay Basin in East China. The Bohai Bay Basin is a complex near-coast Mesozoic and Cenozoic continental rift basin that contains substantial petroleum reserves (Guo et al., 2010; Yu et al., 2015). The Bohai Bay Basin has experienced complex structural evolution. The target area Binnan Oilfield is located in the Lijin Depression in the northwest of the Dongying Depression, adjacent to the Binxian Uplift to the west and the Chenjiazhuang Uplift to the north. Figure 5 shows the structural location map of the study area. The target formation is mainly composed of sandstone, siltstone and thin limestone interlayers, and the porosity ranges from 1% to 17.9% with an average porosity of 9%. The permeability range of the reservoir is 0.1–11.457 mD. According to the mercury injection curve in this area, the pore throat radius is mainly distributed between 0.0025μm and 1.5 μm (Song et al., 2012; Wang et al., 2016). Other simulation parameters are shown in Table 3. It is worth mentioning that the target block is a low-permeability oil reservoir that has not yet been put into production.