Time-dependent model for two-phase flow in ultra-high water-cut reservoirs: Time-varying permeability and relative permeability

Abstract

For the ultra-high water-cut reservoirs, after long-term water injection exploitation, the physical properties of the reservoir change and the heterogeneity of the reservoir becomes increasingly severe, which further aggravates the spatial difference of the flow field. In this study, the displacement experiments were employed to investigate the variations in core permeability, porosity, and relative permeability after a large amount of water injection. A relative permeability endpoint model was proposed by utilizing the alternating conditional expectation (ACE) transformation to describe the variation in relative permeability based on the experimental data. Based on the time dependent models for permeability and relative permeability, the traditional oil-water two-phase model was improved and discretized using the mimetic finite difference method (MFD). The two cases were launched to confirm the validation of the proposed model. The impact of time-varying physical features on reservoir production performance was studied in a real water flooding reservoir. The experimental results indicate that the overall relative permeability curve shifts to the right as water injection increases. This shift corresponds to a transition towards a more hydrophilic wettability and a decrease in residual oil saturation. The endpoint model demonstrates excellent accuracy and can be applied to time-varying simulations of reservoir physics. The impact of variations in permeability and relative permeability on the reservoir production performance yields two distinct outcomes. The time-varying permeability of the reservoir results in intensified water channeling and poor development effects. On the other hand, the time-varying relative permeability enhances the oil phase seepage capacity, facilitating oil displacement. The comprehensive time-varying behavior is the result of the combined influence of these two parameters, which closely resemble the actual conditions observed in oil field exploitation. The time-varying simulation technique of reservoir physical properties proposed in this paper can continuously and stably characterize the dynamic changes of reservoir physical properties during water drive development. This approach ensures the reliability of the simulation results regarding residual oil distribution.

Headnote

ABSTRACT

For the ultra-high water-cut reservoirs, after long-term water injection exploitation, the physical properties of the reservoir change and the heterogeneity of the reservoir becomes increasingly severe, which further aggravates the spatial difference of the flow field. In this study, the displacement experiments were employed to investigate the variations in core permeability, porosity, and relative permeability after a large amount of water injection. A relative permeability endpoint model was proposed by utilizing the alternating conditional expectation (ACE) transformation to describe the variation in relative permeability based on the experimental data. Based on the time dependent models for permeability and relative permeability, the traditional oil-water two-phase model was improved and discretized using the mimetic finite difference method (MFD). The two cases were launched to confirm the validation of the proposed model. The impact of time-varying physical features on reservoir production performance was studied in a real water flooding reservoir. The experimental results indicate that the overall relative permeability curve shifts to the right as water injection increases. This shift corresponds to a transition towards a more hydrophilic wettability and a decrease in residual oil saturation. The endpoint model demonstrates excellent accuracy and can be applied to time-varying simulations of reservoir physics. The impact of variations in permeability and relative permeability on the reservoir production performance yields two distinct outcomes. The time-varying permeability of the reservoir results in intensified water channeling and poor development effects. On the other hand, the time-varying relative permeability enhances the oil phase seepage capacity, facilitating oil displacement. The comprehensive time-varying behavior is the result of the combined influence of these two parameters, which closely resemble the actual conditions observed in oil field exploitation. The time-varying simulation technique of reservoir physical properties proposed in this paper can continuously and stably characterize the dynamic changes of reservoir physical properties during water drive development. This approach ensures the reliability of the simulation results regarding residual oil distribution.

Keywords:

Mimetic finite difference

Water flooding reservoir

Time-varying physical properties

Numerical simulation

(ProQuest: ... denotes formulae omitted.)

1. Introduction

After years of exploitation, a multitude of mature oil fields have now entered the ultra-high water-cut stage, resulting in the dispersal of the residual subsurface crude oil. Consequently, the task of exploration has become arduous. For reservoirs in the ultrahigh water-cut stage, the relentless injection of water over a prolonged period has caused a transformative shift in their physical attributes, leading to an escalation of their heterogeneity and further exacerbating the spatial difference of the flow field (Xu et al., 2010; Guo et al., 2006; Huang et al., 2013; Liu, 2011). Due to the long-term unchanged well pattern and fixed streamline in the injection-production network, a considerable amount of inefficient displacement and even ineffective water circulation have been caused (Jiang et al., 2016a, 2022; Sun et al., 2015; Zhang et al., 2019a, b). Therefore, studying the time-varying reservoir property parameters is crucial for accurately predicting the distribution pattern of remaining oil during high-water-cut stage.

Numerous studies have indicated that long-term water injection leads to changes in reservoir macroscopic, microscopic, and permeability parameters. In terms of macroscopic parameters, * Corresponding author. significant variations are observed in reservoir permeability and clay content, while median porosity and median grain size show relatively minor changes during the water injection process. In medium and low permeability oil reservoirs, favorable properties are enhanced after water flooding in water-wet reservoirs, while adverse properties deteriorate in oil-wet reservoirs. Additionally, the properties of high-porosity and high-permeability reservoirs continue to improve (Xu et al., 2010; Zhou et al., 2010; Zhu et al., 2004). Microscopic parameters include the rock particle framework, pore-throat network, clay minerals, etc. Among them, the pore-throat network is the most complex, exhibiting the greatest variation and having the most significant impact on reservoir fluid flow (Du, 2016; Shi et al., 2013). Furthermore, the content of clay minerals also changes with water injection. Many microscopic experimental studies indicate that, after long-term water flooding, the presence of clay minerals in the pores and on the surfaces of mineral particles significantly decreases in locations with coarser pore-throat diameters (Chen et al., 2012, 2016,; Song et al., 2018). Clay minerals, after hydration and expansion, migrate and eventually accumulate in smaller pore throats, even leading to porethroat plugging (Li et al., 2009; He et al., 2010), especially pronounced during the high-water-cut period in oil fields. Wettability and relative permeability are two of the most influential factors affecting reservoir permeability parameter (Huang et al., 2003; He et al., 2002). Regardless of whether the reservoir is oil-wet or water-wet, its initial water-wetting tendency is enhanced during long-term water injection development. Reservoirs exhibit strong water-wetting characteristics during high-water-cut conditions (Xu et al., 2010; Guo et al., 2006; Meng et al., 2010; Huang et al., 2001). Compared to pre-water flooding, the water phase overall permeability curve of the core shifts to the right after water flooding, with increased bound water saturation and reduced water phase relative permeability at the same water saturation. Oil-water equilibrium points shift to the right, indicating enhanced water-wetting tendencies (Wen et al., 2017; Jackson et al., 2003; Sun et al., 1996).

Although many studies have confirmed that reservoir physical properties may change in the process of water injection, current commercial simulators are unable to quantitatively describe changes in recovery efficiency caused by time-varying reservoir physical properties (Cui and Zhao, 2004; Xu et al., 2015, 2016), which cannot meet the actual demand for effective water injection in reservoir development in high water-cut period, and increases the difficulty of historical fitting (Eydinov et al., 2009; Zahoor and Derahman, 2013). Therefore, it is urgent to develop simulation methods that can characterize the influence of time-varying reservoir physical properties on reservoir development in high water-cut stage. As mathematical models and numerical simulations do not directly account for microscopic parameters, changes in reservoir properties over time are often reflected through variations in macroscopic parameters such as permeability, porosity, and relative permeability curves.

The key to construct time-varying numerical simulation of reservoir is to determine the variation law of reservoir physical parameters in the development process. Therefore, in the process of time-varying numerical simulation of reservoir, a method that can properly characterize the physical property changes is particularly important, and many scholars have carried out research on it. In order to correctly describe the changes of reservoir physical properties, some scholars adopted a staged numerical simulation method to accurately describe the distribution law of remaining oil in the ultra-high water-cut period by considering the changes of rock and fluid parameters (Gai et al., 2000; Gao et al., 2004; Huang et al., 2014). In order to solve the problem of discontinuous physical property changes caused by staged numerical simulation, Ju et al. (2006) and Jiang et al. (1996) proposed a time-varying simulation technology based on a function of water cut. However, this method has the following problems: when the reservoir enters the high water-cut stage, although the water-cut of the reservoir does not change much, injected water will have a strong scour effect on the reservoir, and the simulation parameters of time-varying reservoirs change dramatically, leading to a large simulation error. In addition, many scholars mainly characterize changes in reservoir properties by the cumulative injection of water within a unit volume or area. For example, the characterization methods of water passage multiple (Zhao, 2009; Li, 2007; Xu et al., 2021) and surface flux (Zhang et al., 2019a, b; Qiao, 2017; Jiang et al., 2016b, 2019) are defined as the sum of injected water or natural water passing through a unit pore volume or plane in porous media. Although this method is more reliable in its time-varying characterization of physical properties and overcomes the dependence on mesh size, it ignores the influence of porosity under the same water flux, resulting in inaccurate characterization of water flushing intensity (Lin et al., 2022).

The time-varying mathematical model used for studying the changing properties of reservoirs is primarily a set of equations coupling pressure and concentration equations, with the unknowns to be solved being primarily pressure and its saturation (concentration). Domestic and international research has extensively investigated numerical methods for such models. Over the past two decades, the most commonly used method for reservoir numerical simulation is the finite difference method (Li et al., 2021). In addition, the method of characteristics combined with the finite difference method, known as the characteristic finite difference method, has also been widely used for solving mathematical models of reservoirs (Zhou et al., 2021). While the finite difference method is convenient for handling such problems, it is not adept at addressing complex boundary issues. With the development of the finite element method, it has gradually been applied to numerical solutions in reservoirs (Zhang et al., 2013; Kim and Milind, 2000; Zhang et al., 2018; Noorishad and Mehran, 1982). However, finite element methods do not guarantee local mass conservation (Garipov et al., 2016; Ghorayeb and Firoozabadi, 2000). In contrast to finite element methods, the finite volume method and mixed finite element method ensure local mass conservation. However, for the former, constructing an accurate polygon mesh with a multipoint flux approximation is not easy (Terekhov et al., 2017; Bathe and Zhang, 2002). The latter, although improving algorithm accuracy compared to the finite volume method, reducing computational errors, and ensuring saturation accuracy (Douglas et al., 1983a, 1983b), is only applicable to certain types of elements due to the computation of mixed finite element basis functions. The mimetic finite difference (MFD) method is applicable to any polygon (Brezzi et al., 2005). Due to the robustness and flexibility of the MFD method, its application has significantly increased in recent years, being used in diffusion modeling (Singh, 2010; Lie et al., 2012; Campbell and Shashkov, 2001), twophase flow (Lipnikov et al., 2008), and more recently, state equation component simulations (Abushaikha and Terekhov, 2020). The comparison of various numerical methods is shown in Table 1.

As mentioned above, this paper presents a mathematical model of water flooding with time-varying physical properties based on the ACE method. Additionally, a numerical model utilizing the MFD method is established for numerical discretization. The ACE method utilized in this paper is entirely data-driven, deriving the relationship between dependent and independent variables without the need for a prior assumption of their relationship. Specifically, the prediction effectiveness is notably higher for segments with abrupt changes in parameters. The MFD method exhibits strong local conservation, supports polyhedral mesh, full tensor permeability, and offers the advantages of flexible division and reduced volume self-locking. Hence, it optimizes the modeling and simulation of reservoirs with time-varying properties. In comparison to the traditional time-varying model, the proposed model in this paper addresses issues such as the instability of calculation results and the incompatibility between traditional numerical methods and time-varying models observed in previous studies. The time-varying simulation technique for reservoir physical properties proposed in this paper can consistently and stably characterize the dynamic changes in reservoir physical properties during the water drive development in actual production. Therefore, the purpose of this study: (1) Establish water flooding simulation experiment to explore the permeability and the relation between relative permeability and injected water volume; (2) Based on this relationship, a multi-factor time-varying characterization model of relative permeability was established; (3) A time-varying numerical model is established based on the time-varying characterization model, and a fully implicit mimetic finite difference discretization scheme is adopted to ensure flux conservation and simulation efficiency. The aim of this study is to propose a new numerical simulation method to describe the timevarying mechanism of reservoir properties in water drive reservoirs. The rest of this article is divided into the following sections. In Section 2, the time-varying law of relative permeability is studied experimentally. In Section 3, the variation law of relative permeability is characterized based on alternating condition expectation transformation. In Section 4, the principle of MFD method and the corresponding numerical formula of discrete time-varying model are described in detail. In Section 5, several numerical cases of time-varying permeability and relative permeability are discussed to verify the correctness of the program, and the accuracy of this method in actual oilfield productivity prediction is further introduced. In Section 6, the summary and main conclusions were presented.

2. Experimental measurement and analysis

2.1. Relative permeability experimental measurement

The measurement of the relative permeability curve involves two fundamental laboratory experiments: the steady-state method and the unsteady-state method (Honarpour and Mahmood, 1988). In the steady-state method, a mixture of fluids is injected simultaneously into the core plug, and the relative permeability is determined using the one-dimensional Darcy flow theory. This method is characterized by a straightforward working principle and easy acquisition of parameters. However, it requires a stable pressure difference between the core inlet and outlet, as well as a constant fluid velocity, resulting in a long testing period. In contrast, the unsteady-state method involves injecting one fluid to displace another fluid. This method is based on the onedimensional displacement theory proposed by Buckley-Leverett, which states that the saturation distribution changes over time and distance during the displacement process. The relative permeability of each fluid in the core section varies with saturation, and the velocity of each fluid in the core section changes over time.

Consequently, the relative permeability curves can only be calculated by recording pressure and flow changes for each fluid during the displacement process. While this method significantly reduces the testing time, it necessitates more complex data processing.

This experiment aims to investigate the changes in relative permeability during the process of water flooding. The steady-state method, due to its significant pressure fluctuations and long stabilization time, is not suitable for effectively simulating the water flooding process. Therefore, in this paper, the unsteady-state method is utilized to address this limitation and provide more accurate results. By employing the unsteady-state method, the dynamic changes in relative permeability can be captured, allowing for a better understanding of the behavior of fluids during water flooding.

In this study, the unsteady-state method was employed to measure the oil-water relative permeability at various injected water volume. The experimental procedure adhered to the Chinese industry standard GB/T 28,912-2012. A sandstone core, obtained from a specific reservoir in Shengli Oilfield, was utilized, with a diameter of 2.5 cm and a length of 10 cm. The detailed core parameters are presented in Table 2.

The oil sample used in the experiment was a mixture of crude oil, degassed oil, and kerosene, with a viscosity of 6.25 mPa s and a density of 0.887 g/cm3 . The injection water was prepared based on the measured salinity of the oilfield, which was determined to be 32,000 mg/L, with a viscosity of 0.563 mPa s. The experimental procedure consisted of eight steps, and the experimental setup is illustrated in Fig. 1 (adapted from SY/T 5345-2007).

(a) Prior to testing, the core should be cleaned to remove all original fluid from the sample and dried.

(b) To measure permeability, weigh dry weight, measure length and diameter of rock samples, and measure the density and viscosity of simulated water and oil.

(c) Vacuum dried rock samples, saturate simulated formation water, then weigh the wet weight and calculate the effective pore volume.

(d) The permeability of the core after saturation of formation water should be measured for three consecutive times, and the relative error is less than 3%, and the water relative permeability should be calculated.

(e) Establish bound water saturation and measure oil relative permeability under bound water condition.

(f) Conduct water flooding experiment. When the displacement reaches 10 times the pore volume, determine the water relative permeability under the residual oil and end the experiment.

(g) Wash and dry the core, repeat step (c) to step (f); in step (f), control the water flooding multiple at 300, 600, 900, 1200, 2000 PV, and continue to step (b) to step (f).

(h) Calculate the relative permeability based on the JohnsoneBosslereNaumann method (Johnson et al., 1959).

2.2. Experimental result and analysis

In Fig. 2, the changes in the core after water flooding with different pore volume values are depicted. It is evident from the figure that the remaining oil within the core gradually decreases as the water flooding progresses. Fig. 3 presents the core permeability change multiple curve and porosity change multiple curve. As observed, both the core permeability and porosity show a gradual increase with an increasing volume of injected water. However, the growth rate gradually slows down over time. This can be attributed to the fact that as the injected water volume increases, the migration of particles becomes limited when the displacement velocity remains constant. Consequently, the rate of change in permeability and porosity weakens gradually and eventually reaches a stable state.

As shown in Fig. 4, it is evident that the relative permeability curve shifts to the right after water flooding. This shift indicates an increase in hydrophilicity, meaning that the affinity of the reservoir for water has increased. Furthermore, the residual oil saturation decreases after water flooding, implying that a larger portion of the original oil in the reservoir has been displaced by the injected water. Moreover, there is an increase in bound water saturation, indicating that water has been adsorbed onto the rock matrix and has become bound to it. These observations collectively highlight the significant impact of long-term water flooding on the physical properties and seepage characteristics of the reservoir.

3. Establishment of relative permeability characterization model

In this chapter, the ACE algorithm theory is employed to analyze the relative permeability curves. Five relative permeability curves ranging from 10 to 1200 PV are selected as the sample data. By applying the ACE algorithm, the optimal transformation model based on alternating conditional expectation is obtained. Subsequently, the model is used to predict the relative permeability curve at 2000 PV. The predicted values are then compared with the experimental values to assess the accuracy of the model.

3.1. Alternating condition expectation method

Alternating conditional expectation (ACE) was initially introduced by Breiman and Friedman (1985) and later refined by Xue et al. (1997). This method is a multiple regression approach used to estimate the optimal transformation. The primary objective of this transformation is to maximize the correlation coefficient between a dependent variable and multiple independent variables. What sets ACE apart from other conventional parametric regression methods is its data-driven nature. It relies solely on the data to establish the relationship between dependent and independent variables, without the need for any pre-assumptions regarding their relationship.

The expected regression model of alternating condition is

... (1)

where q is the function of the original dependent variable Y; 4i is a function of the original argument Xi. The n-dimension (X ¼ X1, X2, X3, ..., Xn) between the original dependent variable Y and the independent variable X is calculated by Eq. (1). The nonlinear regression problem is transformed into the problem of finding nþ1 one-dimensional optimal function.

3.2. Normalized mean relative permeability curve

The experiment result indicates that the injected water volume has a great effect on the shape of the relative permeability curves as shown in Fig. 2. In order to describe this influence, it is necessary to first know a standardized relative permeability curve, which can be obtained by dimensionless normalization of the original relative permeability curve. The following formula is used for normalization treatment (Zhang and Song, 2007).

... (2)

... (3)

... (4)

where Sw is the water saturation; Swc is the bound water saturation; Sor is the residual oil saturation; Krw is the water phase relative permeability; Kro is the oil phase relative permeability; KðSor Þ is the relative permeability of water phase under residual oil; KðSwcÞ is the relative permeability of oil phase in bound water; S* w is normalized water saturation; K* rwðS* wÞ is the normalized water phase relative permeability; K* roðS* wÞ is the normalized oil phase relative permeability

S* w is divided into m equal parts from 0 to 1, and K* roðS* wÞ and K* rwðS* wÞ of each relative permeability curve at the same S* w were obtained by interpolation. The average permeability curve of the reservoir was obtained according to Eqs. (5) and (6)

... (5)

... (6)

The normalized value of the data points of the original relative permeability curve was plotted into the normalized coordinate system with S* w as the abscissa and K* ro and K* rw as the ordinate.

Then the data points were returned to obtain the normalized dimensionless oil-water relative permeability curve shown in Fig. 5. The normalized relative permeability curve is regarded as the standard curve. According to the standard curve and the four basic characteristic parameters after changes, the changed relative permeability curve can be obtained in reverse.

3.3. Oil-water relative permeability curve endpoint value model

Since it is challenging to quantify how the relative permeability curve as a whole changes over time, many researchers often use changes in the endpoints of relative permeability curve to calculate the entire relative permeability curve (Verma et al., 1994; Wood et al., 1990; Xu et al., 2016). This research primarily investigates the modifications of the other three points because the relative permeability of the oil phase in bound water is all 1.

Porosity multiples, injection water volume, and permeability multiples were selected as independent variables to predict endpoints of the relative permeability curves, and the end point prediction model was obtained:

... (7)

In Eq. (7), Y is endpoints of the relative permeability curves; F is displacement multiple; Per is permeability multiple; Por is porosity multiple.

By inverting Eq. (7), the calculation relation of the dependent variable can be obtained, which is shown as follow:

... (8)

In the above formula, "100 represents the inverse transformation of 4i *ðXiÞ (i ¼ 1, 2, ..., n).

According to the relative permeability curve samples of water flooding simulation experiment, the optimal transformation relationship between the actual value and the transformed value of bound water saturation, residual oil saturation and the relative permeability of water phase under residual oil is obtained as shown in Fig. 6:

... (9)

... (10)

... (11)

Similarly, as shown in Fig. 7, the optimal transformations of the actual values and transform values of the permeability multiple, porosity multiple and displacement multiple respectively conform to the following relations.

The relationship between the optimal transformation function of the endpoint value of the relative permeability curve and the sum of the optimal transformation function of the parameters is shown in Fig. 8. In the figure, the horizontal axis represents the sum of the transformed values of each independent variable, while the vertical axis represents the corresponding transformed values of the independent variables.

Based on the analysis of Fig. 8, it is evident that a strong linear relationship exists between the optimal transformation function of the endpoint value of the relative permeability curve and the sum of the optimal transformation function of the parameters. Furthermore, the polynomial fitting slope between these two factors is 1, indicating a consistent relationship. Utilizing the polynomial fitting results, a prediction model for the relative permeability curve was established. This model incorporates the optimal transformation functions of the parameters to accurately predict the behavior of the relative permeability curve endpoints.

... (12)

... (13)

... (14)

In Eqs. (12)e(14), bn is the regression coefficient, X is the experimental value of the independent variable, and 4* i ðXiÞ is the optimal transformation expectation function of the independent variable. an is the optimal transformation polynomial fitting function of the independent variable. The regression coefficient is shown in Table 3.

The formulas for calculating the end points of the relative permeability curve are shown as follows:

... (15)

... (16)

... (17)

3.4. Validation of model

Based on the endpoint value model of the oil-water relative permeability curve and in combination with the normalized average relative permeability curve of the reservoir, a prediction model for the oil-water relative permeability curve is developed. The following are the basic steps involved in establishing this model:

(1) Firstly, according to several representative oil-water relative permeability curves, an average relative permeability curve is obtained by using the normalization method.

(2) According to the porosity multiple Por, air permeability multiple Per and displacement multiple F of a certain reservoir area, bound water saturation, residual oil saturation and water relative permeability under residual oil saturation of this area are obtained by using the end model of relative permeability curve.

(3) Eq. (18) was used to inversely calculate the water saturation Sw corresponding to the normalized water saturation S* w in the region.

... (18)

(4) By using the inverse transformation Eqs. (19) and (20), the oil phase relative permeability Kro and water phase relative permeability Krw under water saturation S* w were inversely calculated:

... (19)

... (20)

In order to validate the accuracy of the relative permeability prediction model established by the ACE method, the calculated results of the model were compared with the actual relative permeability results obtained from the water flooding simulation experiment. We calculate the root mean square error (RMSE) between the experimental relative permeability curve and the predicted relative permeability curve to quantify the overall performance of the model. The error between the experimental curve and the predicted curve is shown in Fig. 9. The calculated value of the relative permeability prediction model is in good agreement with the actual value, especially the prediction accuracy of the section with sharp changes in water content is significantly improved. The relative error between the calculated value and the actual value is small, which indicates that the prediction of the model has high accuracy.

4. Mimetic finite difference method

Considering the time-varying law of permeability and relative permeability, a time-varying reservoir model based on ACE characterization is established. The mathematical model is discretized by MFD method, and the discretization process is described in detail.

4.1. Governing equations

We consider a two-phase fluid flow system for a reservoir porous media. The governing equations are:

... (21)

... (22)

... (23)

... (24)

where So, Sw are saturation of phase oil, water respectively, %; KðtÞ is the absolute permeability vector, um2; Kro, Krw are the relative permeability of phase oil, water respectively; mo, mw are the viscosity of phase oil, water respectively, mPa s; po, pw are the pressure of phase oil, water respectively, MPa; ro, rw are the density of phase oil, water respectively, kg=m3; qvo, qvw are source of oil and water respectively, m3=s; vo, vw are velocity of phase oil and water respectively, m=s; g is gravity acceleration, m=s2; z is vertical depth, m; 4 is porosity, %;

For simplicity, we use the relative mobility as

... (25)

... (26)

Combining Eqs. (25) and (26) the velocity v:

... (27)

... (28)

The calculation process of the simulator is shown in Fig. 10. Compared with common reservoir numerical simulation software, the calculation process of each time step of the simulator is as follows: 1) According to the permeability and relative permeability variation law, calculate the conductivity of each grid; 2) Using the mimetic finite difference method to solve the pressure equation, calculate the water phase flow in each direction of each grid, and then obtain the directional plane flux and total surface flux of each grid; 3) Finite volume method was used to solve the saturation equation to obtain the water saturation of each grid. 4) Recalculate the permeability and relative permeability of each grid.

4.2. Mimetic finite difference method

Let U be a domain composed by non-overlapping polyhedrons UH ¼ fUig. For any Ui, Ui2UH, Ak ¼ Ui∩Uj is the interface of Ui and Uj. nk ¼ jAkjnbk is the area weighted normal vector, nbk is unit normal vector. First, we define the pi is the pressure at the centroid xi of element Ui and tk is the interface pressure at the centroid xk of interface Ak, as shown in Fig. 11.

By using MFD method to solve the pressure equation, the relation of Ui 2UH can be obtained as follows:

... (29)

... (30)

According to Darcy's law, normal flux u on the interface can be rewritten as:

... (31)

where ui ¼ ½u1; u2; ...; um T is the vector of fluxes on the interface; m is the number of interfaces of element Ui; L is the matrix of relative permeabilities; ei ¼ ð1; ...; 1Þ T 1m; Xi is the scalar quantity at the centroid of element Ui; ti is the vector of scalar quantities defined at the centroid of interfaces of element Ui; Ti is a positive definite matrix, which is a key part of MFD method. The quantities for MFD are shown in Fig. 11.

In order to get Ti, the scalar quantity Xi is supposed to be linear. Thus, we can get the expressions as follows:

... (32)

where x is the coordinate vector of centroid of element. The vector a and scalar b are constant vector and constant scalar quantity.

The paper examines the inclusion of pressure and gravity, along with their respective magnitudes at the centroid and interface. These quantities adhere to the following relationships:

... (33)

... (34)

where zi is the depth at the centroid of element Ui. tz;k;i is the interface depth at interface Ak. The coordinate vector of the centroid of element Ui is represented by xi , while the coordinate vector of the centroid of interface Ak is denoted as xk.

According to Darcy's law, the flux through the interface Ak of element Ui is

... (35)

where v is the velocity; vk is the Darcy velocity on the interface Ak; | Ak | is the area of interface Ak; nbk is the outward unit normal vector on interface Ak. Combining the Darcy velocity equation, the interface flux can be rewritten as:

... (36)

Based on the linear relationships (Eqs. (33) and (34)), substituting Eqs. (35) and (36) into Eq. (31) yields the following: Thus, we can reach the matrix form: TiX ¼ NK (38) where X ¼ ½x1 xi，...，xk xi，...，xm xi T; N ¼ ½jA1jnb1; ...; jAkjnbk; ...; jAmjnbm T Defining XðiÞ represent i-th coordinate of X, then NTX ij ¼ Xm k¼1jAkjnbk ðiÞ ðxk xiÞ ðjÞ (39) According to the divergence theorem, the above equation yields the following, NTX ij ¼ Xm k¼1jAkjnbk ðiÞ 1 jAkj ð Ak ðxk xiÞ ðjÞ ds ¼ Xm k¼1nbk ðiÞ ð Ak ðxk xiÞ ðjÞ ds ¼ Xm k¼1 bei ð Ak ðxk xiÞ ðjÞ nbk ds ¼ bei ð Ui Vðx xiÞ ðjÞ dU ¼ bei be T j jUij ¼ dijjUij (40) where jUij is the volume of element U i; bei ¼ 0，...， 1 ith ，...，0

; dij is Delta function. From Eq. (40), we can get NTX ¼ jU ijI (41) where I is d d identity matrix. Multiplying Eq. (38) with NT on both sides yields,

... (42)

Generally, in order to guarantee the existence of the inverse matrix of Ti, Ti is composed by two parts:

... (43)

...

Thus, we can reach the matrix form:

... (37)

... (38)

... (39)

According to the divergence theorem, the above equation yields the following,

... (40)

... (41)

where I is d d identity matrix. Multiplying Eq. (38) with NT on both sides yields,

... (42)

Generally, in order to guarantee the existence of the inverse matrix of Ti, Ti is composed by two parts:

... (43)

...

... (44)

where d denotes the spatial dimension; Ai is a diagonal matrix composed by interface areas. Qi ¼ orthðAiXiÞ is orthonormal basis, where Xi is a diagonal matrix composed by the distance from centroids to interface centroids. trðKÞ is the trace of K. Ii is an identity matrix.

After calculating the local matrix Ti, finally the interface flux ua;i of phase a can be rewritten as

... (45)

where ua represents the interface fluxes of phase a; T denotes the global matrix composed of fTig; Pa and z are element pressure vectors, capillary pressure vector and depth vector for elements, respectively; ta, x are interface pressure vector, capillary pressure vector and depth vector for interfaces, respectively; C is a block diagonal matrix, which i-th column is given by bei; D is related to the interface. The number of unit entries in the columns of D depends on the properties of interfaces. If element number i contain interface number j, then the corresponding entry Dij of D is 1, otherwise Dij is 0.

Considering the velocity continuity condition on the boundary surface of the element, the numerical calculation format of MFD can be obtained as follows:

... (46)

seepage velocity array; p is the element center pressure array; t is the central pressure array of the unit boundary surface; q ¼ ½qi is a unit source sink term of Ui; m is the total number of grid cells.

5. Applied computational tests

In this section, the accuracy of a time-varying simulator developed using the MFD method was demonstrated. Two injectionproduction models, namely the quarter five-point well pattern and the inverted five-spot well pattern, were established to assess the performance of the simulator under different scenarios. The models were considered with different variations: no time-varying properties, time-varying permeability, time-varying relative permeability, and combined time-varying permeability and relative permeability.

5.1. Validation tests results

Test 1: Quarter-five-point well pattern reservoir simulation.

In the first example, as shown in Fig. 12, we analyzed four different scenarios with varying reservoir properties over time using a quarter-five-point for production. The scheme details are provided in Table 4.

As shown in Fig. 13(a)e(d), the simulation results reveal significant differences among the four simulation scenarios. In Case 1, which does not consider time-varying effects, the remaining oil is distributed relatively uniformly throughout the reservoir. The flow lines are spread out, and the pressure distribution is relatively consistent. However, this scenario does not capture the impact of time-varying factors on reservoir behavior. In Case 2, where timevarying permeability is considered, the remaining oil tends to accumulate on both sides of the dominant channel in the reservoir. The flow lines are concentrated within this channel, and the pressure in the channel is affected by fluid pressure loss, resulting in a pressure drop. This effect highlights the importance of accounting for time-varying permeability in accurately modeling fluid flow behavior. Case 3 incorporates time-varying relative permeability, which improves the displacement efficiency of water flooding and reduces the residual oil saturation in the reservoir. Consequently, the remaining oil, pressure, and flow lines become more evenly distributed throughout the reservoir. In Case 4, both time-varying permeability and time-varying relative permeability are considered. This scenario results in the formation of a dominant channel, similar to Case 2, but the channel is relatively wider and exhibits lower remaining oil saturation compared to Case 2.

In Fig. 14(a) and (b), the recovery efficiencies of the four schemes are 59.54%, 52.29%, 61.02%, and 55.31%, respectively, while the water-free oil recovery periods for the four schemes are 1752, 1569, 2044, and 1715 d, respectively. It is evident that the increase in the water-free oil recovery period leads to an increase in the recovery rate. This can be explained from two perspectives. On one hand, the long-term water injection flushes increase the permeability of the reservoir, resulting in the formation of dominant channels within the reservoir, exacerbating the water cone phenomenon, and causing production wells to extract injected water earlier. On the other hand, changes in relative permeability reduce the water-oil mobility ratio, slowing down the flow of injected water, making oil easier to displace, thereby prolonging the water-free oil recovery period.

Test 2: Inverted five-spot reservoir simulation.

In the second example, we considered inverted five-point well pattern. Fig. 15 shows the numerical conceptual model of one-well injection, four-well production grid. The injection rate at the injection well qin is 0.004 PV/d, and the liquid production rate at the production well qout is 0.001 PV/d. The simulation time is ten years.

As shown in Figs. 16 and 17, the results of the inverted five-point well pattern are presented for the four different cases. Similar to the one injection one production pattern, the presence of dominant channels is observed in the case of time-varying permeability (Case 2). These dominant channels concentrate the remaining oil at the edges of the reservoir, making it difficult to displace effectively. In contrast, in the case of time-varying relative permeability (Case 3), the displacement of the reservoir leads to a more uniform distribution of the remaining oil content. The impact of time-varying physical properties in the reservoir is intensified due to the large injection and production rates in the inverted five-point well pattern. Fig. 17 shows the degree of recovery in each case, with Case 2 achieving a recovery efficiency of 49%, Case 3 reaching 70.6%, and Case 4 at 51.1%. Comparing these results to the one injection one production pattern, it is evident that the influence of permeability and relative permeability is enhanced in the inverted five-point well pattern. Specifically, the development effect is worse when considering time-varying permeability, while considering timevarying relative permeability leads to a better development effect.

5.2. Field study

K oilfield is a marine sandstone oilfield with abundant edge and bottom water energy. The porosity of the reservoir is 20.1%e30.5%, and the permeability is 52e157 mD. The distribution of permeability and well location are shown in Fig. 18. The oil has good physical property. The formation oil viscosity is 2.7 mPa s, and the oil volume coefficient is 1.05. The formation water viscosity is 0.36 mPa s, and the formation water volume coefficient is 1.028. The bound water saturation is 0.253, and the geological reserves used in the reservoir are 4130 104 m3 . The oil field has been steadily developed for 23 years and adopts the high-speed development strategy. The average daily fluid volume of a single well is up to 158 m3 . The current water-cut has reached 97%, and it is in the ultra-high water-cut period. At present, there are problems such as large water consumption, low water storage and difficult follow-up development.

As shown in Fig. 19, at the early stage of development, the high speed and low potential area is first formed near the water wells. As the water-cut increases, the low velocity and high potential area, where oil production potential is higher, tends to be limited to regions between the oil wells or the water wells. This suggests that the remaining oil is mainly concentrated in these areas, while the non-producing reserves areas become more dominant. When the water-cut reaches a high level, such as 95%, an invalid water circulation area forms near the water wells. This area represents regions where water circulation becomes less effective in displacing oil due to factors such as high-water saturation or increased permeability. In this stage, the production contribution from these areas may become negligible.

5.3. Discussion

The time-varying permeability and time-varying relative permeability have distinct effects on the reservoir. The changes in permeability caused by the erosion of injected water enhance the flow capacity of the dominant channels along the main stream lines, leading to the premature breakthrough of injected water into the production well. This negatively impacts reservoir development by reducing the efficiency of oil displacement. On the other hand, the flushing effect of water flooding alters the relative permeability curve. The changes in the endpoint of the relative permeability curve result in an increase in oil phase relative permeability and a decrease in water phase relative permeability at the same water saturation. This enhancement of oil flow capacity and weakening of water flow capacity, combined with a reduced water-oil mobility ratio, improve the sweep area of water flooding. Furthermore, the reduction in residual oil saturation leads to improved displacement efficiency and an increased degree of recovery. It is evident that time-varying permeability has a negative effect on reservoir development, while time-varying relative permeability has a positive effect. However, the impact of permeability is more significant compared to relative permeability. As a result, when considering the comprehensive time-varying effects, the production effect of the oilfield is still worse. Therefore, it is crucial to consider the timevarying effects comprehensively in reservoir simulation in order to better reflect the real reservoir conditions and optimize reservoir development strategies.

In summary, under two different productions well networks, the simulator can show good adaptability, and the simulation results meet the expected design requirements under different timevarying schemes, and the calculation results are accurate and effective. Meanwhile, as shown in Figs. 13 and 16, in the process of homogeneous simulation, the saturation field, pressure field and streamline field, the model calculation results of either the quarter well pattern or the reverse five-point well pattern have good symmetry, indicating that the simulator is reliable.

The new simulator's simulation data shows a high degree of matching with the actual production data of the field. This indicates that the new simulator is effective in accurately fitting the historical production behavior of the oilfield, particularly during the late stage of long-term water drive development. This accurate historical fitting is valuable for guiding future production measures and strategies in the oilfield. Compared to simulations that do not consider time variation, the new simulator demonstrates improved fitting effects for water rates and cumulative liquid production. Additionally, the predicted distribution of remaining oil aligns well with the actual production results in the early time steps. This indicates that simulation models can provide more accurate oilfield production prediction results, enhancing the understanding of reservoir behavior and supporting decision-making processes.

6. Conclusions

In this study, a time-varying prediction model is established based on core displacement experiments to investigate the effects of long-term waterflooding on reservoir physical properties and development performance. The following conclusions are drawn:

Long-term waterflooding leads to changes in core physical properties. Permeability and porosity increase logarithmically with the increase of injected water volume. The residual oil saturation decreases, and the hydrophilicity of the reservoir strengthens. The original relative permeability curve also undergoes changes.

An end value model of the oil-water relative permeability curve is established using the alternating condition expectation method. By combining this model with the normalized average relative permeability curve of the reservoir, a prediction model for the oilwater relative permeability curve is formed. The predicted values of the model closely match the experimental values, indicating its reliability.

Considering the time-varying effects of reservoir permeability and relative permeability, a new reservoir time-varying program was developed based on the mimetic finite difference method. A quarter-five-spot well pattern and a reverse five-spot well pattern were used to investigate and analyze the impact of different timevarying factors on development dynamics and residual oil distribution, thereby validating the correctness of the program. The newly developed simulator has helped improve the history matching of water flooding reservoirs in the K oilfield, thereby further enhancing the rationality and reliability of production forecasting. Overall, the history matching is relatively good.

CRediT authorship contribution statement

Shao-Chun Wang: Writing e original draft. Na Zhang: Conceptualization, Data curation, Funding acquisition, Methodology, Writing e original draft, Writing e review & editing. Zhi-Hao Tang: Data curation, Writing e review & editing. Xue-Fei Zou: Validation, Writing e review & editing. Qian Sun: Data curation, Writing e original draft, Writing e review & editing. Wei Liu: Conceptualization, Writing e original draft.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work is supported by Research project of Shengli Oifield Exploration and Development Research Institute (Grant No. 30200018-21-ZC0613-0125).

Sidebar

ARTICLE INFO

Article history:

Received 26 July 2023

Received in revised form

18 November 2023

Accepted 28 March 2024

Available online 12 April 2024

Edited by Meng-Jiao Zhou

Footnote

* Corresponding author E-mail address: [email protected] (N. Zhang).

References

References

Abushaikha, A.S., Terekhov, K.M., 2020. A fully implicit mimetic finite difference scheme for general purpose subsurface reservoir simulation with full tensor permeability. J. Comput. Phys. 406, 109194. https://doi.org/10.1016/ j.jcp.2019.109194.

Bathe, K.J., Zhang, H., 2002. A flow-condition-based interpolation finite element procedure for incompressible fluid flows. Comput. Struct. 80 (14e15), 1267e1277. https://doi.org/10.1016/S0045-7949(02)00077-9.

Breiman, L., Friedman, J.H., 1985. Estimating optimal transformations for multiple regression and correlation. J. Am. Stat. Assoc. 80 (391), 580e598. https:// doi.org/10.1080/01621459.1985.10478157.

Brezzi, F., Lipnikov, K., Shashkov, M., 2005. Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes. SIAM J. Numer. Anal. 43 (5), 1872e1896. https://doi.org/10.1137/04061395.

Campbell, J., Shashkov, M., 2001. A tensor artificial viscosity using a mimetic finite difference algorithm. J. Comput. Phys. 172 (2), 739e765. https://doi.org/10.1006/ jcph.2001.6856.

Chen, X., Liu, F., Wang, B., et al., 2012. Permeability changes under long-term water flooding and microscopic model study. Spec Oil Gas Reservoirs 19 (1), 66e69þ138. https://doi.org/10.3969/j.issn.1006-6535.2012.01.015 (in Chinese).

Chen, D., Li, J., Zhu, W., et al., 2016. Experimental research on reservoir parameters variation after water flooding for offshore unconsolidated sandstone heavy oil reservoirs. China Offshore Oil Gas 28 (5), 54e60. https://doi.org/10.11935/ j.issn.1673-1506.2016.05.009 (in Chinese).

Cui, C., Zhao, X., 2004. Numerical simulation method and its application when considering time-dependent effect of reservoir parameters. Hydrodynam Res Prog (Ser A) (S1), 912e915. https://doi.org/10.3969/j.issn.10004874.2004.z1.012 (in Chinese).

Douglas, J., Ewing, R.E., Wheeler, M.F., 1983a. The approximation of the pressure by a mixed method in the simulation of miscible displacement. RAIRO - Analyse Numerique 17 (1), 17 e33. http://www.numdam.org/item/M2AN_1983__17_1_ 17_0/.

Douglas, J., Ewing, R.E., Wheeler, M.F., 1983b. A time-discretization for a mixed finite element approximation of miscible displacement in porous media. RAIRO - Analyse Numerique 17 (3). http://www.numdam.org/item?id¼M2AN_1983__ 17_3_249_0.

Du, Q., 2016. Variation law and microscopic mechanism of permeability in sandstone reservoir during long-term water flooding development. Acta Pet. Sin. 37 (9), 1159e1164. https://doi.org/10.7623/syxb201609010 (in Chinese).

Eydinov, D., Gao, G., Li, G., et al., 2009. Simultaneous estimation of relative permeability and porosity/permeability fields by history matching production data. J. Can. Petrol. Technol. 48 (12), 13e25. https://doi.org/10.2118/132159-PA.

Gai, Y., Lv, D., Guo, Y., 2000. Numerical simulation by stages about the reservoir at high water cut period. Petroleum Geology and Recovery Efficiency 7 (1), 54e56. https://doi.org/10.13673/j.cnki.cn37-1359/te.2000.01.016 (in Chinese).

Gao, B., Peng, S., Huang, S., 2004. Numerical simulation of reservoir in the composition stage of sand layer 3 in Sha2 Member 2, Shengtuo Oilfield. Petrol. Explor. Dev. 31 (6), 82e84. https://doi.org/10.3321/j.issn:10000747.2004.06.021 (in Chinese).

Garipov, T., Karimi-Fard, M., Tchelepi, H., 2016. Discrete fracture model for coupled flow and geomechanics. Comput. Geosci. 20, 149e160. https://doi.org/10.1007/ s10596-015-9554-z.

Ghorayeb, K., Firoozabadi, A., 2000. Numerical study of natural convection and diffusion in fractured porous media. SPE J. 5, 12e20. https://doi.org/10.2118/ 51347-PA.

Guo, L., Wang, Y., Liu, W., et al., 2006. Variation law of reservoir parameters during waterflooding in Dagang oil field. Petroleum Geology & Experiment 28 (1), 85e90. https://doi.org/10.3969/j.issn.1001-6112.2006.01.017 (in Chinese).

He, F., Yu, T., Zhang, J., et al., 2002. Effects of water flooding on properties of reservoir and formation crude oil. J. Daqing Pet. Inst. 26 (2), 21e23þ133. https://doi.org/10.3969/j.issn.2095-4107.2002.02.007 (in Chinses).

He, W., Tan, W., Ma, C., 2010. Law and mechanism of characteristic changes of ultralow permeability reservoirs before and after waterflooding. J. Oil Gas Technol. (J. Jianghan Petroleum Inst.) 32 (5), 56e59þ400. https://doi.org/10.3969/ j.issn.1000-9752.2010.05.012 (in Chinese).

Honarpour, M., Mahmood, S.M., 1988. Relative-permeability measurements: an overview. J. Petrol. Technol. 40 (8), 963e966. https://doi.org/10.2118/18565-PA.

Huang, W., Wang, J., Fan, Z., 2001. A study on the relative permeability curves of the west area of Block 7 in Gudong oil field, 90-91þ11-1 Petrol. Explor. Dev. (3). https://doi.org/10.3321/j.issn:1000-0747.2001.03.027 (in Chinese).

Huang, W., Xu, H., Lu, X., 2003. Analysis of influencing factors of injected water on reservoir flow performance. Oil Gas Field Surface Eng 22 (8), 54. https://doi.org/ 10.3969/j.issn.1006-6896.2003.08.044 (in Chinese).

Huang, Y., Tang, Y., Li, L., et al., 2013. Changes in reservoirs with mid-low permeability after water flooding and analysis of influencing factors, 143-147þ0-1þ0 Journal of Oil and Gas Technology 35 (8). https://doi.org/10.3969/j.issn.10009752.2013.08.033 (in Chinese).

Huang, Y., Zhou, Z., Liu, Z., et al., 2014. Method for improving the numerical simulation accuracy of the ultra-high water cut reservoir in Puannan Three fault blocks. Petroleum Geology and Recovery Efficiency 21 (5), 65e68þ115. https:// doi.org/10.3969/j.issn.1009-9603.2014.05.015 (in Chinese).

Jackson, M.D., Valvatne, P.H., Blunt, M.J., 2003. Prediction of wettability variation and its impact on flow using pore to reservoir-scale simulations. J. Petrol. Sci. Eng. 39 (3), 231e246. https://doi.org/10.1016/S0920-4105(03)00065-2.

Jiang, R., Chen, Y., Deng, Y., et al., 1996. Numerical simulation of the reservoir physical properties in the second block of Shengtuo oil field. OGRT 3 (2), 50e56 (in Chinese).

Jiang, R., Qiao, X., Teng, W., et al., 2016a. Impact of physical properties time variation on waterflooding reservoir development. Fault-Block Oil Gas Field 23 (6), 768e771. https://doi.org/10.6056/dkyqt201606017 (in Chinese).

Jiang, R., Qiao, X., Teng, W., et al., 2016b. Numerical simulation of reservoir time variation based on surface flux. Special Oil Gas Reservoirs 23 (2), 69e72þ154. https://doi.org/10.3969/j.issn.1006-6535.2016.02.016 (in Chinese).

Jiang, R., Cui, Y., Hu, Y., et al., 2019. Numerical simulation of polymer flooding considering reservoir property time variation. Fault-Block Oil Gas Field 26 (6), 751e755. https://doi.org/10.6056/dkyqt201906015 (in Chinese).

Jiang, Q., Qu, Y., Wu, J., et al., 2022. Research progress of superior reservoir channels and interwell connectivity. Liaoning Chem. Ind. 51 (2), 261e265. https:// doi.org/10.14029/j.cnki.issn1004-0935.2022.02.011 (in Chinese).

Johnson, E.F., Bossler, D.P., Bossler, V.O.N., 1959. Calculation of relative permeability from displacement experiments. Transactions of the AIME 216 (1), 370e372. https://doi.org/10.2118/1023-G.

Ju, B., Fan, T., Zhang, J., et al., 2006. Variation of oil viscosity in water drive reservoirs and its influence on development results. Petrol. Explor. Dev. 33 (1), 99e102. https://doi.org/10.3321/j.issn:1000-0747.2006.01.023 (in Chinese).

Kim, J.G., Milind, M., 2000. Finite element, discrete-fracture model for multiphase flow in porous media. AIChE J. 46 (6), 1120e1130. https://doi.org/10.1002/ aic.690460604.

Li, X., 2007. Numerical Modeling Study of Reservoir Physical Characteristic Parameters Change after Waterflood Development Master thesis. China University of Petroleum. https://doi.org/10.7666/d.y1214411 (in Chinese).

Li, H., Wang, X., Liu, S., 2009. Study on variation of reservoir physical property parameters in old oil fields. J. Southwest Petrol. Univ. (Nat. Sci. Ed.) 31 (2), 85e89þ184e185. doi:10.3863/j.issn.1674-5086.2009.02.022 (in Chinses).

Li, R., Huang, Y., Li, Y., 2021. Research progress of reservoir seepage model and numerical simulation technology. Petrochemical Applications 40 (7), 11e15þ29. https://doi.org/10.3969/j.issn.1673-5285.2021.07.003 (in Chinese).

Lie, K.A., Krogstad, S., Ligaarden, I.S., et al., 2012. Open-source Matlab implementation of consistent discretisations on complex grids. Comput. Geosci. 16, 297e322. https://doi.org/10.1007/s10596-011-9244-4.

Lin, J., Jiang, R., Shen, Z., et al., 2022. Comprehensive characterization investigation of multiple time-varying rock-fluid properties in waterflooding development. ASME Journal of Energy Resources Technology 144 (7), 073001. https://doi.org/ 10.1115/1.4052166.

Lipnikov, K., Moulton, J.D., Svyatskiy, D., 2008. A multilevel multiscale mimetic (M3) method for two-phase flows in porous media. J. Comput. Phys. 227 (14), 6727e6753. https://doi.org/10.1016/j.jcp.2008.03.029.

Liu, X., 2011. Study on numerical simulation technology based on time varying physical properties. Petroleum Geology and Recovery Efficiency 18 (5), 58e62þ115. https://doi.org/10.3969/j.issn.1009-9603.2011.05.016 (in Chinese).

Meng, L., Ren, B., Ju, B., 2010. Research on water drive law of fault block reservoir in high water cut stage. J. Southwest Petrol. Univ. (Nat. Sci. Ed.) 32 (6), 139e142þ195. https://doi.org/10.3863/j.issn.1674-5086.2010.06.028 (in Chinese).

Noorishad, J., Mehran, M., 1982. An upstream finite element method for solution of transient transport equation in fractured porous media. Water Resour. Res. 18 (3), 588e596. https://doi.org/10.1029/WR018i003p00588.

Qiao, X., 2017. Research on Time-Varying Numerical Simulation and Flow Field Evaluation of Reservoir Physical Properties in Water Drive Reservoir. Master thesis (in Chinese), China University of Petroleum (East China).

Shi, C., Zhang, F., Chen, P., 2013. Affection of simulating water-flooding by water injection tests on reservoir properties. J Southwest Petrol Univ (Nat Sci Ed) 35 (5), 87e93. https://doi.org/10.3863/j.issn.1674-5086.2013.05.012 (in Chinese).

Singh, G., 2010. Mimetic finite difference method on gpu: application in reservoir simulation and well modeling. Master thesis, Institutt for Matematiske Fag.

Song, R., Wang, Y., Liu, J., 2018. Microscopic pore structure characterization and fluids transport visualization of reservoir rock. J. Southwest Petrol. Univ. (Sci. Technol. Ed.) 40 (6), 85e105. https://doi.org/10.11885/j.issn.16745086.2018.07.18.03 (in Chinese).

Sun, S., Han, J., Guo, Y., et al., 1996. Laboratory experiment on physical properties of flooding sandstone in Shengtuo Oilfield. J. Univ. Pet. (China) 20 (S1), 33e35 (in Chinese).

Sun, T., Xiong, Y., Wang, L., et al., 2015. Time-varying law of reservoir physical propertiesdpucheng oilfield Wen51 block lower shaer reservoir for example, pucheng. Sci. Technol. Eng. 15 (36), 52e55þ62. https://doi.org/10.3969/ j.issn.1671-1815.2015.36.009 (in Chinese).

Terekhov, K.M., Mallison, B.T., Tchelepi, H.A., 2017. Cell-centered nonlinear finitevolume methods for the heterogeneous anisotropic diffusion problem. J. Comput. Phys. 330, 245e267. https://doi.org/10.1016/j.jcp.2016.11.010.

Verma, M.K., Boucherit, M., Bouvier, L., 1994. Evaluation of residual oil saturation after waterflood in a carbonate reservoir. SPE Reservoir Eng. 9 (4), 247e253. https://doi.org/10.2118/21371-PA.

Wen, X., Dai, Z., Wang, H., et al., 2017. Physical properties of marine sandstone reservoir after long-term waterflooding. Spec Oil Gas Reservoirs 24 (1), 157e161. https://doi.org/10.3969/j.issn.1006-6535.2017.01.032 (in Chinses).

Wood, K.N., Tang, J.S., Luckasavitch, R.J., 1990. Interwell residual oil saturation at Leduc miscible pilot. In: SPE Annual Technical Conference and Exhibition. https://doi.org/10.2118/20543-MS.

Xu, C., Zhang, R., Ma, L., et al., 2010. Reservoir dynamic variation characteristics after water flooding and its influencing factors. Lithologic Reservoirs 22 (S1), 89e92. https://doi.org/10.3969/j.issn.1673-8926.2010.Z1.019 (in Chinese).

Xu, Q., Chen, Y., Hou, Y., et al., 2015. Research on time-varying numerical simulation of late-stage reservoirs with extremely high water cut. Drill. Prod. Technol. 38 (5), 41e43þ8. https://doi.org/10.3969/J.ISSN.1006-768X.2015.05.13 (in Chinese).

Xu, J., Guo, C., Jiang, R., et al., 2016. Study on relative permeability characteristics affected by displacement pressure gradient: experimental study and numerical simulation. Fuel 163, 314e323. https://doi.org/10.1016/j.fuel.2015.09.049.

Xu, W., Li, L., Chen, Z., 2021. Numerical simulation method based on pore network simulation and time-varying reservoir physical properties. Pet. Eng. Constr. 47 (S2), 173e178. https://doi.org/10.3969/j.issn.1001-2206.2021.S2.034 (in Chinese).

Xue, G., Datta-Gupta, A., Valko, P., et al., 1997. Optimal transformations for multiple regression: application to permeability estimation from well logs. SPE Form. Eval. 12 (2), 85e93. https://doi.org/10.2118/35412-PA.

Zahoor, M.K., Derahman, M.N., 2013. New approach for improved history matching while incorporating wettability variations in a sandstone reservoirdfield implementation. J. Petrol. Sci. Eng. 104, 27e37. https://doi.org/10.1016/ j.petrol.2013.03.025.

Zhang, J., Song, K., 2007. Relative permeability characteristic curve and its application. Acta Pet. Sin. 28 (4), 104. https://doi.org/10.3321/j.issn:02532697.2007.04.021 (in Chinese).

Zhang, N., Yao, J., Huang, Z., et al., 2013. Accurate multiscale finite element method for numerical simulation of two-phase flow in fractured media using discretefracture model. J. Comput. Phys. 242, 420e438. https://doi.org/10.1016/ j.jcp.2012.12.006.

Zhang, N., Wang, Y., Sun, Q., et al., 2018. Multiscale mass transfer coupling of triplecontinuum and discrete fractures for flow simulation in fractured vuggy porous media. Int. J. Heat Mass Tran. 116 (JAN.), 484e495. https://doi.org/10.1016/ j.ijheatmasstransfer.2017.09.046.

Zhang, F., Jiang, R., Cui, Y., et al., 2019. Numerical simulation method of fractured oil reservoir with time varying physical properties. Special Oil Gas Reservoirs 26 (4), 103e108. https://doi.org/10.3969/j.issn.1006-6535.2019.04.018 (in Chinese).

Zhang, W., Liu, B., Wang, X., et al., 2019. Research on the evolution of high conductive channels based on formation parameters variation. Xinjiang Oil and Gas 15 (4), 61e66þ4. https://doi.org/10.3969/j.issn.1673-2677.2019.04.013 (in Chinese).

Zhao, S., 2009. Numerical Simulation of Reservoir Physical Property Parameters in Waterflood Development Master thesis. China University of Petroleum. https:// doi.org/10.7666/d.y1542781 (in Chinese).

Zhou, Y., Hu, W., Zhang, Y., et al., 2010. Law of reservoir physical property changes both before and after waterfloodingdusing Es3 reservoir of south fault block Hu7 in huzhuangji oilfield for example. J. Oil Gas Technol. (J. Jianghan Petroleum Inst.) 32 (5), 25e30þ401. https://doi.org/10.3969/j.issn.10009752.2010.05.006 (in Chinese).

Zhou, Z., Hang, T., Jiang, T., et al., 2021. Mass conservative characteristic finite difference method for convectionediffusion equations. Int. J. Comput. Math. 98 (10), 2115e2136. https://doi.org/10.1080/00207160.2021.1876229.

Zhu, H., Du, Q., Li, Z., et al., 2004. Variation of physical characteristics and wettability of the reservoir in high water-cut stage. Petrol. Explor. Dev. (S1), 82e84. https:// doi.org/10.3321/j.issn:1000-0747.2004.z1.019 (in Chinese).