A Multi-Scale Numerical Simulation Method Considering Anisotropic Relative Permeability. (2025)

Author(s): Li Wu [1,2]; Junqiang Wang [1]; Deli Jia [2]; Ruichao Zhang [1,3]; Jiqun Zhang [2]; Yiqun Yan [2]; Shuoliang Wang (corresponding author) [1,*]

1. Introduction

Most of China’s oil reservoirs are fluvial deposits with firm reservoir heterogeneity. The differences in pore structure and fluid flow capacity in different directions clearly cannot be ignored. There have been many studies on the directional characterization of pore structure; however, the research on reservoir anisotropy, whether through physical experiments or mathematical modeling, mainly focuses on the anisotropy of absolute permeability. The idea of anisotropy in absolute permeability is widely recognized [1,2]. However, reservoir fluid is a multiphase flow, and the relative permeability is an essential primary parameter to characterize the fluid movement in numerical simulation. The shape of the relative permeability curve significantly affects the water cut, oil production, and residual oil distribution. Both theoretical and experimental studies have shown that relative permeability is an orientation-dependent tensor rather than a scalar, like absolute permeability. Since the discovery by Corey and Rathjens in 1956 that the relative permeability of the layered rocks was directional [3,4], a variety of scholars have investigated the anisotropy of relative permeability using core experiments, pore network modeling, and microscopic visualization of flow physics simulations [5,6,7]. The results show that relative permeability is closely related to the pore throat structure and microfractures, regardless of the representative elementary volume scale, laboratory scale, or field scale [8,9,10]. This anisotropy prevails and affects the value of relative permeability, making the phase flow ability in different directions inconsistent [11]. If the anisotropy is not considered in the numerical simulation, the multiphase flow process in the subsurface cannot be faithfully captured.

Nevertheless, although the anisotropic character of relative permeability has been recognized earlier, its application in numerical simulation still needs improvement. Commercial numerical simulation software (e.g., ECLIPSE 2022-12 and CMG-STARS 2022), widely used nowadays, typically offers the following methods for processing relative permeability. A common approach involves using a normalized relative permeability curve derived from representative samples to represent the entire zone. Additionally, assigning different relative permeability curves based on sedimentary microphase is a frequent practice. Calibration of the relative permeability curve is often performed using initial saturation as the reference endpoint. In all cases, relative permeability is consistently treated as a scalar, neglecting the directional variability. This simplification introduces significant bias in multiphase flow modeling, even though numerical simulation remains one of the most advanced and practical techniques in reservoir research. Existing commercial simulation software fails to account for the anisotropy of relative permeability, a critical factor in accurate reservoir modeling.

In actual reservoirs, the complexity of geological structures often hinders efficient reserve mobilization. Numerical simulation is a critical tool for reservoir management, with the accuracy of simulation models being essential for reliable forecasting. Theoretically, these models are validated against analytical solutions and laboratory experiments; however, in practical oilfield applications, achieving model accuracy requires a detailed history-matching process. This process adjusts the simulated reservoir state to closely reflect actual conditions by refining well logging data and geological parameters and compensating for incomplete historical data. Physics-based methods are continually applied to update the understanding of geological structures and reservoir properties. The quality of history matching is crucial for generating accurate future predictions and identifying zones with residual oil regions [12]. In recent years, Li and Wang et al. have established reservoir simulators considering anisotropic relative permeability based on experiment data [13]. The study by Li et al. shows that numerical simulation methods incorporating anisotropic absolute permeability and anisotropic relative permeability produce more accurate history matching for natural reservoirs. These findings highlight that considering these specific mechanisms significantly improves the accuracy of reservoir models [14,15,16,17]. Nevertheless, the study by Li and Wang only integrated the experimentally derived anisotropic relative permeability to the control equation, employing fully implicit discretization and the Gauss–Seidel iteration without improving the solution method. Including anisotropic relative permeability will substantially increase the complexity of numerical simulation [18]. As a result, although the accuracy of the simulator developed by Li and Wang has been improved, maintaining computational efficiency remains challenging [19]. Global unconventional oil resources are estimated to be comparable in volume to conventional resources. As conventional oil and gas fields progress into their late stages of development, unconventional resources are assuming a more critical role in production. Interest in unconventional oil and gas has surged, and the rapid growth of unconventional oil and gas development in recent years has led to significant advancements in technologies like horizontal drilling and multi-stage fracturing. Numerical simulation of unconventional hydrocarbons has become a prominent focus in the industry, driven by the need for more accurate modeling and efficient development [20,21,22]. The development of unconventional hydrocarbons faces many challenges, including the macroscopic and microscopic heterogeneity of reservoirs, multi-scale and multi-medium characteristics, complex fluid properties, state characteristics as well as complex flow characteristics and coupled exploitation mechanisms under unconventional development modes [19,23,24]. For example, the numerical simulation of hydraulic fracture propagation in unconventional oil and gas fields is a challenging computational problem because it involves multiple physical processes [25]. Typically, the CPU time is the limiting factor when the simulation model is established. Hence, there is an urgent need to develop multi-scale numerical simulation technologies to improve simulation efficiency while considering complex mechanisms.

To build simulators that consider relative permeability anisotropy while improving simulation efficiency, this paper proposed an improved multi-scale finite volume (IMsFV) method. It developed a multi-scale numerical simulator that enables the consideration of reservoir anisotropy. Compared with the traditional procedure, the Egg model verifies that the newly developed simulator can substantially speed up the simulation while guaranteeing computational accuracy. Moreover, the built-up multi-scale simulator demonstrated an influence of anisotropy relative permeability on waterflooding. Therefore, the method can be used as an effective tool for fine simulation and accurate residual oil prediction in anisotropic reservoirs, which is of great significance for development in the waterflooding field.

2. Methods

This study introduces a multi-scale reservoir simulation method incorporating anisotropic relative permeability. Initially, improvements were made to both the internal reservoir flow model and the well model based on the oil–water two-phase black-oil model, specifically considering anisotropic relative permeability. Then, the model was decoupled using a sequential solution into pressure and transport equations. An IMsFV method was employed for the pressure equation, with prolongation and restriction operators used to enable iterative solving. The transport equation was solved through flux reconstruction. The main technical roadmap is illustrated in Figure 1.

2.1. Numerical Simulation with Anisotropic Relative Permeability

2.1.1. Anisotropic Relative Permeability

Relative permeability, defined as the ratio of effective permeability to absolute permeability (Equation (1)), is a fundamental parameter crucial for characterizing the behavior of multiphase fluid flow. (1)K[sub.ra]=Ka/K where K[sub.ra] is the relative permeability of phase a (a represents oil or water); K[sub.a] is the permeability of phase a; and K is the absolute permeability. In addition, in this paper, K[sub.ra(x)],K[sub.ra(y)],K[sub.ra(z)] represent the relative permeability of phase a and S[sub.or(x)],S[sub.or(y)],S[sub.or(z)] represent the residual oil saturation in the x-, y-, and z-directions, respectively.

Relative permeability directly affects the accuracy and reliability of reservoir numerical simulation results. In complex reservoir conditions, pore structure, fluid dynamics, and rock–fluid interactions introduce differences in relative permeability, making it a direction-dependent tensor. Figure 2 shows the test results for relative permeability from a core sample, while Table 1 specifies the absolute permeability values for the core in each direction and provides the endpoints of the phase percolation; the specific testing method refers to the research paper by Li et al. from our team [13].

Table 1 and Figure 2 clearly show a positive correlation between relative and absolute permeability within the same core. The relative permeability of both oil and water is directional, with higher residual oil saturation (S[sub.or]) occurring in the direction of greater permeability. Residual oil saturation increases sequentially in the x, y, and z order, while water relative permeability of residual oil saturation (k[sub.rw](S[sub.or])) decreases sequentially in the x, z, and y order. When 1-S[sub.or(x)]<S[sub.w]<1-S[sub.or(y)], oil can no longer flow along the x-direction but it can still flow in the y- and z-directions, leading to a continued decrease in oil saturation. Concurrently, as water relative permeability in the x-direction increases, this results in significant water outflow along the x-direction, while oil is restricted to seepage in the y–z plane. When 1-S[sub.or(y)]<S[sub.w]<1-S[sub.or(z)], oil loses the ability to flow in the x- and y-directions and the water relative permeability continues to rise. When S[sub.w]=1-S[sub.or(z)], oil completely loses the ability to flow in all directions, and the system reaches the “residual oil state” where oil cannot flow.

2.1.2. Fundamental Mathematical Model

In this paper, a two-phase black-oil model for oil and water is developed. In the model, both oil and water phases can be represented as a single component, the properties of which only change with variations in temperature and pressure. However, the internal composition of the component remains unchanged, and the basic assumptions for modeling are as follows.

(1) Only oil and water phases are preset in the reservoir;

(2) No mass exchange between oil and water;

(3) Fluid flow conforms to Darcy’s law;

(4) Reservoir is isothermal;

(5) Phase equilibrium completed instantaneously.

The control equations for oil and water are, respectively, given by Equations (2):(2)?(?basa)/?t=-?·(b[sub.a]??[sub.a])+b[sub.a]q[sub.a] where ? is the porosity, fraction; b[sub.a] is the inverse formation value factor of phase a; s[sub.a] is the saturation of phase a, volume fraction; ??[sub.a] are the Darcy velocities of phase a, cm/s; q[sub.w] represents the volumetric production or injection rate, respectively, s[sup.–1].

The phase velocities can be obtained based on Darcy’s law, (3)??[sub.a]=-kraK/µa(?p[sub.a]-?[sub.a]g?z) where µ[sub.a] is the viscosity of phase a, kg/ms; p[sub.a] is the pressure of phase a, bar; ?[sub.a] is the density of phase a, kg/m[sup.3]; g is the gravity acceleration, kg·m/s[sup.2]; z is the vertical depth, m.

The number of unknown variables in Equations (2)–(3) exceeds the number of equations, therefore, subsidiary Equations (4) and (5)—for saturation and capillary force—are provided to satisfy the closure condition for the unknown variable. (4)?as[sub.a]=1 (5)p[sub.o]-p[sub.w]=p[sub.cow] where the capillary pressure between oil–water p[sub.cow] is a function of saturation.

2.1.3. The Black-Oil Model Considering Anisotropic Relative Permeability

The MATLAB Reservoir Simulation Toolbox (MRST 2022b), provided by MathWorks in Natick, US, was used for conducting the studies. The basic principle of MRST is to discretize the mass conservation equation and then characterize it in divergence form, using the upstream weights method and gradient operator to compute the flow and gradient, respectively. After determining the flow velocity at the grid surface, the divergence operator is used to compute the flow velocity for each grid block. The final step involves accumulating the inflow–outflow within each block to balance the mass conservation [26,27]. Conventional relative permeability in the equations is replaced to assess the impact of anisotropic relative permeability on development in the black-oil model.

Flow within the Reservoir

If anisotropic relative permeability is ignored, the grid mobility can be obtained directly by the upstream weighting method. The relative permeability of each grid cell is expanded from one to three (K[sub.ra(x)],K[sub.ra(y)],K[sub.ra(z)]). Since mobility is defined as the ratio of permeability to viscosity, the mobility of each grid cell changes from a single value to three distinct values. This article employs the finite volume method to calculate the fluxes. However, the grid block has three velocity components, and only the component aligned with the face’s direction is utilized to compute the flux for each grid cell face. Therefore, when considering anisotropy and using the upstream weighting method to calculate the mobility of each grid cell, it is necessary to identify the upstream grid block and the orientation of the grid surface. Ultimately, only the mobility in one specific direction is retained for each grid cell.

For example, if the upstream grid face lies in the x-direction, the mobility in the x-direction of that grid is retained, and the relative permeability in the x-direction of the upstream grid is used for Equation (6). (6)K[sub.ra,ij]={K[sub.ra(x),i]p[sub.i]=p[sub.j]K[sub.ra(x),j]p[sub.i]<p[sub.j]

The same applies to the y- and z-directions; see Equation (7) and Figure 3. (7)K[sub.ra,ij]={K[sub.ra(y),i]p[sub.i]=p[sub.j]K[sub.ra(y),j]p[sub.i]<p[sub.j]orK[sub.ra,ij]={K[sub.ra(z),i]p[sub.i]=p[sub.j]K[sub.ra(z),j]p[sub.i]<p[sub.j] where K[sub.ra,ij] represents the relative permeability of phase a in cell (i, j).

Well Models

The well modeling follows the Peaceman model [28], the most commonly used in commercial numerical simulators, see Equation (8). (8)q[sub.ai]=?[sub.ai]WI[sub.i][(p[sub.ai]-p[sub.bh])-?[sub.ai]g(z[sub.i]-z[sub.ref])] where q[sub.ai] is the volumetric production or injection rate of phase a in cell i; s[sup.–1] ?[sub.ai] is the mobility of phase a in cell i, mD·kg/ms; WI[sub.i] is the well index in cell i; p[sub.ai] is the pressure of phase a in cell i, bar; p[sub.bh] is the bottom hole pressure, bar; ?[sub.ai] is the density of phase a in cell i, kg/m[sup.3]; z[sub.i] is the vertical depth in cell i, m; z[sub.ref] is the vertical reference depth, m.

The mobility and well index WI are calculated by Equations (9) and (10). (9)?[sub.ai]=krai/µai (10)WI[sub.i]=2pkh/ln(reqrw)

For each well grid, the effective relative permeability needs to be calculated. In this paper, well perforation in different directions is calculated using the harmonic mean method, as shown in Equations (11)–(13), respectively, which is consistent with the Peaceman model. (11)k[sub.ra]=[square root of k[sub.ra(y)]k[sub.ra(z)]] (12)k[sub.ra]=[square root of k[sub.ra(x)]k[sub.ra(z)]] (13)k[sub.ra]=[square root of k[sub.ra(x)]k[sub.ra(y)]] where K[sub.ra(x)],K[sub.ra(y)],K[sub.ra(z)] represent the relative permeability of phase a in the x-, y-, and z-directions, respectively.

2.2. Sequential Splitting

An additional challenge is achieving an efficient solution for the model once special mechanisms are considered. Commercial numerical simulators primarily use fully implicit discretization to transform a nonlinear system into a linearized form, which is then solved using Jacobi or Gauss–Seidel iterative methods. However, as simulation considerations become more complex and grid resolution increases, the complexity of the numerical simulations grows significantly, making efficiency requirements increasingly demanding. Wallis and Gries et al. [29,30] substituted an iterative solver with a Constrained Pressure Residual (CPR) preconditioner for a direct solver to obtain higher computational performance. However, the CPR preconditioner is ineffective, according to Møyner and Cusini et al. [31,32]. A sequential solution method is applied to decouple the black-oil equations (Equations (2)–(3)) into distinct pressure and transport equations, enabling a multi-scale solution for the pressure equation.

2.2.1. The Pressure Equation

The pressure equation is first obtained by applying fully implicit discretization to Equations (2). (14)R[sub.a]=1/?t[(?basa)[sup.n+1]-(?basa)[sup.n]]+?·(ba??a)[sup.n+1]-(baqa)[sup.n+1]=0 where R[sub.a] represents the pressure equations in the residual form of phase a.

The saturation-independent pressure equation is obtained from Equation (4), (15)R[sub.p]=Rw/bwn+1+Ro/bon+1=0

Then, considering the Peaceman well model as shown in Equation (8), the following equations are formulated:(16)(Rq,a)[sub.w]=(qas)[sub.w]-?c?C[sub.w]b[sub.a][c]q[sub.a][c]=0 (17)(RC)[sub.w]=p[sub.bh,w]-C[sub.bh,w]or(RC)[sub.w]=q[sub.w][sup.s]-C[sub.q,w] where R[sub.q,a] represents well production equations in the residual form of phase a; C[sub.bh,w] and C[sub.q,w] represent the well control conditions for each well.

Based on Equations (14)–(17), the pressure equation in residual form O[sub.e] could be written as follows:(18)O[sub.e]=[Rp,Rq,w,Rq,o,RC][sup.T]=0

Using the oil-phase pressure p[sub.o] as the primary variable, the Newton–Raphson method is employed to solve for the unknown state quantities e=[po,qws,qos,pbh][sup.T] in the pressure equation.

2.2.2. The Transport Equation

To derive the transport equations for updating saturation, it is assumed that the state variables from the pressure equation, which satisfy the convergence conditions, have already been obtained. Additionally, the relationship between the total flux ?? and the phase flux ??[sub.a] is governed by ??=?a??[sub.a].

Based on these state variables, the total flux is calculated by summing Darcy’s law for each phase using the pressure equation in Equation (18) and substituting it into Equation (19) to determine the phase velocity in the conservation equation. (19)??[sub.a]=f[sub.a](??+K???[sub.?](??g?z+?P[sub.c])) where f[sub.a] is partial flow coefficient, fraction.

The phase velocity from Equation (19) is substituted into Equations (2), thereby solving the saturation of the different phases. In the case of the oil–water two-phase black-oil model, the saturation of the oil and water phases must satisfy the condition that their sum equals one. Therefore, it is only necessary to solve for oil saturation. Let d=[so][sup.T], where d represents the column vector composed of the primary variables of the transport equation.

2.2.3. Solving the Nonlinear Problem

The Newton–Raphson method is used for the iterative solution. It is assumed that the state variables (e[sup.n], d[sup.n]) for the current time step have been determined. The computation processes for the subsequent time steps are as follows.

Step 1. Set the iteration counters ? and ? for the pressure and transport equations, both initialized to zero. The state of the current timestep is used as the initial guess for the solution, that is, e[sub.0]=e[sup.n] and d[sub.0]=d[sup.n], for the pressure equation and the transport equation, respectively;

Step 2. Solve the linearized pressure equation and update the pressure variables accordingly. (20)e[sub.?+1]=e[sub.?]-J[sub.e](e?,d?,en,dn)[sup.-1]O[sub.e](e[sub.?],d[sub.?],e[sup.n],d[sup.n]) where J[sub.e] refers to the Jacobian of the residual equation O[sub.e];

Step 3. Verify whether the pressure iteration converges or not by using the following criteria. (21)||p?+1-p?||[sub.8]=e[sub.p](maxp[sub.?]-minp[sub.?])||R[sub.q]||[sub.8]=e[sub.q]?R[sub.C]?[sub.8]=e[sub.C]

Set ?=?+1. If convergence has not been reached, repeat Step 2. Otherwise, proceed to the next step after calculating the total flux ??[sub.?];

Step 4. Solve the linearized transport system based on the updated pressure. (22)d[sub.?+1]=d[sub.?]-J[sub.d](e?,d?,???,en,dn)[sup.-1]O[sub.d](e[sub.?],d[sub.?],??[sub.?],e[sup.n],d[sup.n])

Check whether the saturation solution has converged using the maximum and the total error for the phases under consideration;

Step 5. Evaluate whether the saturation iteration has converged using the following criteria. (23)?tRa?[sub.8]=e[sub.V]b[sub.a][sup.avg]||?tRa||1/||?||1=e[sub.M]b[sub.a][sup.avg]

Set ?=?+1. Repeat Step 4 if convergence has not been achieved; otherwise, proceed to the next step;

Step 6. Optionally, check the pressure residual against a specified tolerance for the outer iteration, as follows:(24)||O[sub.e](e[sub.?],d[sub.?],e[sup.n],d[sup.n])||[sub.8]=e[sub.it] where O[sub.e] represents the pressure equations in residual form.

If not converged, return to Step 2. Otherwise, complete the timestep.

The pressure and saturation at each time step are determined using the iterative process described above. The pressure is calculated by solving the system of discrete equations using the Newton–Raphson method, as outlined in Equation (20). The computational cost primarily involves constructing the Jacobi matrix and performing iterative updates of the variable increments.

2.3. Multi-Scale Formulation for Pressure Equations

In reservoir numerical modeling, the choice of time discretization significantly affects pressure solutions. Commonly used approaches include fully implicit methods, implicit pressure explicit saturation (IMPES) schemes, and sequential splitting techniques designed to linearize the nonlinear system. To solve the pressure equation system, iterative algorithms like Jacobi iteration, Gauss–Seidel iteration, Newton–Raphson, conjugate gradient, and generalized minimal residual methods are employed. The fundamental goal is to iteratively reduce the residual until convergence to a sufficiently small value is achieved. Two primary methods are used to improve the convergence rate of the iterative solution of the pressure equation system. Based on the predict–correct concept, the first method involves the CPR operator. During the initial phase of the iterative process, a decoupled linear equation system is solved to determine the coefficients, variables, and source terms that depend solely on pressure. This forms a reduced-order equation system, which is solved as the prediction phase and combined with the correction phase. The second method, inspired by the divide and conquer principle, utilizes the Algebraic Multigrid Method (AMG) [33]. AMG enhances the convergence rate of the iterative algorithm by transferring information across grids at different levels. Scholars such as Jenny, Lunati, and Møyner have proposed multi-scale methods based on this AMG framework. The paper adopts an IMsFV method for solving the nonlinear pressure equation to accelerate the iteration loop. This multi-scale method requires two sets of grids: a fine grid based on the static properties of the reservoir {Oif}[sub.i=1][sup.nf]—on which the models in this paper are built—and a coarse grid {Ojc}[sub.j=1][sup.nc] derived from the fine grid, as shown in Figure 4. The coarse mesh consists of multiple fine meshes, each associated with only one coarse mesh. The construction of the multi-scale grids follows the approach outlined by Olav et al. [34,35].

This paper focuses on developing a multi-scale method for the pressure equation. For a rapid solution of the transport equation, the domain decomposition method proposed by Kozlova et al. [36] may be consulted.

2.3.1. Multi-Scale Solver for Pressure

The primary solution object for the pressure equation is the pressure increment. Equation (20) is first re-expressed to develop a multi-scale solver for pressure. (25)J?p=r where J is the Jacobian matrix of discrete equations; r is the vector of residuals.

Since Equation (25) may exhibit significant coefficient variations and considerable coupling between adjacent cells, a direct solution is avoided. Instead, an approximation of the pressure is pursued using a multi-scale expansion, enabling an efficient solution on a coarser grid system {Ojc}[sub.j=1][sup.nc]. The process is as follows:

Assume the existence of a prolongation operator P and a restriction operator R, based on the concept of algebraic multigrid. First, the prolongation operator establishes the relationship between the pressure increment at the fine and coarse scales. (26)?p˜?p[sub.f]=P?p where P is the prolongation operator (matrix); ?p[sub.f] are fine-scale pressure increments.

Solving fine-scale pressure increments is then reformulated as solving for coarse-scale pressure increments. (27)(RJP)?p[sub.c]=J[sub.c]?p[sub.c]=Rr=r[sub.c] where R is restriction operator (matrix); ?p[sub.c] is a coarse-scale increment.

2.3.2. The Restriction and Iterative Prolongation Operators

The critical aspect of the procedure is the construction of the prolongation operator P and restriction operator R.

Lunati and Lee et al. [37] demonstrated the relationship between the basis functions and the prolongation operator, showing that the regular arrangement of full-domain basis functions represents the prolongation operator. A classic method for constructing a multi-scale prolongation operator P involves concatenating solutions to a series of locally computed flows [38,39]. For example, in the classical MsFV method, a series of local flows are solved to construct a prolongation operator by creating a dual mesh of coarse-scale grids [40]. In this paper, Restricted-Smoothed Basis Functions are employed to iteratively construct the prolongation operator P, following the approaches of Møyner and Lie, and the prolongation operator and the basis functions are considered equivalent. This method avoids complex setups for flow problems constrained by boundary conditions and addresses the MsFV method’s limitation in constructing high-quality basis functions on unstructured grids.

Assumed that the matrix used for iteration has the property that the row sum is zero and the initial prolongation operator has the property that the row sum is one. Therefore, (28)J=1/2{J+J[sup.T]-diag[(J+J[sup.T])1]} where 1 denotes a column vector with elements equal to 1, and diag is an operator that expands the vector into a diagonal matrix of the same size.

Given the following initial conditions, the initial basis functions of the coarse grid are set to a constant value within the grid and zero elsewhere, i.e., P[sup.0]=R[sup.T]. The following procedure is applied to iteratively update P. (29)P[sup.n+1]=P[sup.n]-M(?Ds[sup.-1]JP[sup.n]) where ? is a relaxation factor; D[sub.s] is the smoother constructed from the diagonal elements of the matrix J. The restrictor M limits the smoothing of the basis functions to a support region bounded by a core of coarse-scale grids and eight peripheral grid centroids, ensuring the accuracy of the prolongation.

For the restriction operator R, a straightforward construction involves using the control volume restriction:(30)(Rc?)[sub.ji]={1,ifx[sub.i]?O[sub.j][sup.b]0,otherwise

Jenny applied this restriction operator in an MsFV method and demonstrated that it ensures mass conservation [40]. In addition, Wang et al. proposed the Galerkin restriction operator R[sub.G]=P[sup.T] to improve convergence speed [41]. Some researchers combine the characteristics of both operators by applying the Galerkin restriction operator in the early stage and the standard restriction operator in the later stage to ensure mass conservation. This combination strategy also applies to the multi-scale method proposed in this paper. However, as Møyner and Lie [26,42] demonstrated, both operators exhibit the same convergence performance when prolongation operators are constructed using Restricted-Smoothed Basis Functions. Therefore, the restriction operators in this paper are built by Equation (30).

2.3.3. Iterative Multi-Scale and Flux Reconstruction

Pressures computed by multi-scale methods typically approximate fine-scale pressures as the prolongation operator does not account for the changes in density and mobility caused by pressure variations. Direct use of the prolongation operator from Equation (29) to calculate fine-scale pressure increments introduces a certain degree of relative error. Therefore, this paper adopts an iterative approach to update the pressure increments. The multi-scale pressure approximation is iteratively corrected using ILU(0) as a smoother, following the method proposed by Wang et al. [41].

The defect is expressed as follows:(31)r[sup.n]=J?p[sup.n]

Let s[sup.n]=U[sup.-1](L[sup.-1]r[sup.n]), the result of using the ILU(0) smoother to the defect with an initial guess of zero can be expressed as outlined below:(32)?p[sup.n+1]=?p[sup.n]+P[J[sub.c][sup.-1]R(r[sup.n]-Js[sup.n])]+s[sup.n]

The fine-scale pressure can be updated by applying a small number of multi-scale iterations at each nonlinear iteration step using Equation (32), thereby eliminating the need for the Newton–Raphson method to solve the linearized system. Once the pressure changes satisfy the conditions in Equation (24) and the multi-scale pressure solver is complete, the fluxes on the fine-scale grid can be reconstructed from the approximate pressures.

According to Jenny et al. [37,40], the multi-scale method ensures conservation on coarse-scale grids but has difficulty guaranteeing mass conservation under fine-scale approximate pressure fields. To solve this problem, local flow with a Newman boundary is solved on the coarse-scale grid to reconstruct a pressure field that ensures mass conservation in the fine-scale system.

Face fluxes are corrected based on whether the fine-scale grids connected to the grid faces belong to the boundary of the coarse grid or the interior grid [26,31]. (33)?^[f]={?(e[sub.?-1],?p^,d[sub.?]),ifB[C[sub.1](f)]=B[C[sub.2](f)]?(e[sub.?-1],?p[sub.?],d[sub.?]),ifB[C[sub.1](f)]?B[C[sub.2](f)] where f represents a specific grid face; C[sub.1] and C[sub.2] represent the grids connected to that face; and B determines the coarse-scale grid to which the grid belongs; p[sub.?] is the approximate pressure obtained through the multi-scale solution method described in Section 2.3.1; p^ is the reconstruction pressure, which is needed to solve the local flow problem with the Neumann boundary conditions based on p[sub.?]. The form of the operator for p^ can be found in studies by Møyner et al.

Once the flux reconstruction is finished, return to step 3 in Section 2.2.3 to continue the iteration for computing the transport equations.

3. Results and Discussion

This section builds on the methodology outlined in Section 2, utilizing MRST (2022b) to develop a multi-scale numerical simulator that accounts for relative permeability anisotropy.

3.1. Model Validation

The results of the traditional fully implicit method (FI) and the multi-scale method (MS) are compared using the Egg model. The Egg model—an open geological model widely adopted by researchers—is a benchmark for reservoir simulation, optimization, of reservoir development plans, and inverse modeling. It is primarily utilized to simulate two-phase oil–water flow, with waterflooding being the main production mechanism due to the absence of an aquifer and gas cap [43,44,45,46]. It is a heterogeneous waterflooding reservoir model with eight injection and four production wells, and the model’s permeability and well distribution are shown in Figure 5. A detailed description of the model is provided by J. D. Jansen [47]. The model contains 18,553 active grids. In the multi-scale solver, the 60 × 60 × 7 fine-scale grid is uniformly coarsened to a 12 × 12 × 1 coarse-scale grid, resulting in a scale-up factor of 5 × 7.

3.1.1. Accuracy

The method’s computational accuracy is demonstrated by comparing it with fully implicit fine-scale solutions using the Egg model.

Figure 6 and Figure 7 display the pressure fields solved at fine scales using the traditional fully implicit method (FI) compared to those solved using the multi-scale method (MS) established in this paper at simulation times of 6 months, 1 year, and 3 years, respectively. The field maps reveal excellent agreement between the traditional fully implicit method and the multi-scale method. This confirms the computational reliability of the multi-scale approach for the nonlinear pressure equations.

Figure 8 and Figure 9 present the oil saturation solutions obtained from the traditional fully implicit and multi-scale methods, respectively. The minimal difference between the two results demonstrates the robustness and accuracy of the multi-scale solver.

In addition to the field maps, the differences in the production curves of each well between the two methods are compared. Taking the PROD2 as an example, Figure 10 and Figure 11 present the daily oil production and water cut curves obtained by the two methods. The curves for PROD2 from both methods are nearly identical, demonstrating a good fit between the results. Although there is some deviation in the oil production and water cut curves of PROD2 in the early stage, the overall shape and trend are similar. Statistical analysis reveals a mean relative error of 9.79% in daily oil production and 7.12% in water cut, providing precise measures of variability in the respective processes. This comparison further verifies the reliability of the proposed multi-scale method under isotropic relative permeability conditions.

3.1.2. Efficiency

This section uses the Egg model to investigate the efficiency of the multi-scale approach to the nonlinear pressure equation. The efficiency is demonstrated by comparing the simulation times.

As shown in Figure 12, for each iteration time step, the linear solution time and the simulation time for the multi-scale method (MS) exhibit a substantial advantage over the traditional fully implicit method (FI). The multi-scale method substantially improves the iterative solution speed of the model.

Table 2 compares the total computation time for the traditional fully implicit (FI) and the multiscale (MS) method. The multiscale method demonstrates a significant improvement in total computation efficiency. Specifically, the simulation time is reduced to less than 30%, and the linear solution time is reduced to less than 1/20 of that required by the fully implicit method under the same simulation conditions.

3.2. Numerical Results with Anisotropic Relative Permeability

In this sub-section, the developed simulator detailed the effect of anisotropic relative permeability on waterflooding.

In the base case, a heterogeneous one-quarter Five-Spot Problem in a 3D reservoir model was considered, where a producer and an injector were positioned at the diagonal corners. The model size was 250 × 250 m. The 50 × 50 × 10 fine-scale grids were uniformly coarsened to 10 × 10 × 2 coarse-scale grids, resulting in an upscaling factor of 5 ×5. Other reservoir parameters are presented in Table 3.

To investigate oil–water displacement and residual oil distribution under conditions of absolute permeability anisotropy, relative permeability anisotropy, and their combined effects, a comparative analysis of five cases was conducted using x- and y-directional anisotropy as an example. Case 1 represents an isotropic homogeneous model as the baseline for comparison. Case 2 considers only absolute permeability anisotropy (Kx > Ky), while Case 3 focuses solely on relative permeability anisotropy. Case 4 addresses absolute and relative permeability anisotropy, with Kx > Ky, and Case 5 examines dual anisotropy but with Kx < Ky. The specific configurations are presented in Table 4.

In all cases, the reservoir is initially saturated with oil, and water is injected to displace the oil. Water is injected from the lower left corner, and oil is produced from the upper right corner. The water injection rate is maintained at a constant 30 m[sup.3]/day, and the reservoir fluid is produced at the same rate under reservoir conditions.

The anisotropic relative permeability curves in the simulation are shown in Figure 13. In all the cases of isotropic permeability, the x-direction curves in Figure 13 are used.

The relative permeability curve indicates that S[sub.wc(x)] ˜ S[sub.wc(y)] = S[sub.wc], and S[sub.or(x)] > S[sub.or(y)], as shown in Table 5.

As previously discussed in Section 2.3, scholars like Jenny, Lunati, and Møyner have proposed various multi-scale solution methods, and comparisons of pressure and saturation field maps are techniques frequently employed by researchers. For instance, Zhou and colleagues utilized the top layer of the SPE 10 model to analyze the Line Drive scenario; concurrently, one-quarter of a five-spot injection problem within a 2D reservoir model was used to investigate the distribution of pressure and saturation at different dimensionless characteristic times [48]. Møyner et al. validated the effectiveness of the proposed multi-scale method using the SPE10 model and the Gullfaks field. Subsequently, saturation profiles from a two-dimensional waterflooding example were employed to analyze changes in the characterization of the water drive front over different time intervals [42]. Chaabi and Cusini et al., taking SPE10 as an example, compared the saturation and pressure field distributions at each time step of three multi-scale methods with those of a fine-scale numerical simulation method. The results demonstrate that these multi-scale methods provide highly accurate approximations of the saturation calculation outcomes from the fine-scale simulations [49,50]. The comparison of saturation field maps from different methods or the same method applied to various models at different times has consistently demonstrated that multi-scale solution methods closely approximate fine-scale solutions. This highlights the superiority of the multi-scale approaches.

Based on the simulation results across different cases, the distributions of pressure (Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18) and oil saturation (Figure 19, Figure 20, Figure 21, Figure 22 and Figure 23) were compared at 3 months, 1 year, 3 years, and at the time of water breakthrough (t[sub.wb]).

The pressure and oil saturation distribution results across the five cases, as shown in Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, Figure 19, Figure 20, Figure 21, Figure 22 and Figure 23, exhibit significant differences. Both absolute and relative permeability anisotropy substantially impact the distribution of pressure and oil saturation. These findings align with those reported by Li but contrast with the conclusions of Gomez-Hernandez, who stated that pressure is not significantly influenced by relative permeability anisotropy [13].

In Case 1, depicted in Figure 14, the pressure between the injection and production wells is uniformly distributed, with the pressure field exhibiting strict symmetry along the diagonal connecting the injection and production wells. A uniform water flood occurs in this case, resulting in a saturation field, and the residual oil is symmetrically distributed along the diagonal, as shown in Figure 19.

In Case 2, a pronounced counterclockwise distortion is observed at each time point compared to Case 1, as shown in Figure 15, due to the directionality of absolute permeability. The absolute permeability in the x-direction is higher than in the y-direction, leading to more rapid pressure transmission in the x-direction from the start. Figure 20 demonstrates that the injected water preferentially flows in the x-direction, resulting in faster and more complete flooding. At the same time, more residual oil accumulates in the y-direction, where the permeability is lower. The anisotropy of the absolute permeability significantly impacts the pressure and oil saturation fields throughout the production cycle.

In Case 3, the pressure distribution, as shown in Figure 16, deviates from the symmetry observed in Case 1 but exhibits distinct behavior compared to Case 2. The counterclockwise distortion of the pressure field becomes apparent primarily in the middle and late stages (after t = 3 years), while the early-stage pressure field remains relatively similar to that of Case 1. Additionally, the degree of pressure field distortion in this case is significantly less pronounced than in Case 2. The saturation field, as illustrated in Figure 21 follows a similar pattern. In the early stage, the saturation field closely resembles that of Case 1. Still, in the middle and late stages, the displacement in the x-direction accelerates, causing more residual oil to accumulate in the y-direction. Unlike the pressure field, the variability in the saturation field and residual oil distribution in the later stages is greater than in Case 2. The relative permeability curve indicates that S[sub.wc(x)] ˜ S[sub.wc(y)] = S[sub.wc], and S[sub.or(x)] > S[sub.or(y)], as shown in Figure 13. In the early stage, when water saturation S[sub.w] is less than 1- S[sub.or(x)], oil and water flow in the x- and y-directions. Water shows slightly higher mobility in the x-direction, resulting in marginally faster pressure transmission, though this difference remains subtle. In the middle and later stages, when S[sub.w] exceeds 1- S[sub.or(x)], a substantial volume of water flows out in the x-direction, leading to much faster pressure transmission than the y-direction. This effect appears as a counterclockwise distortion in the pressure distribution. The saturation profile reveals that waterflooding progresses more rapidly along the x-direction while residual oil accumulates in the y-direction. Anisotropic relative permeability shifts pressure propagation during water injection toward regions of higher water-phase relative permeability, leaving more residual oil in areas with lower water-phase relative permeability.

In summary, the analyses suggest that anisotropic absolute permeability and anisotropic relative permeability significantly affect pressure distribution and oil–water flooding through different mechanisms. While both factors influence the entire production cycle, relative permeability anisotropy has a more pronounced impact during the middle and high water cut periods. Absolute permeability anisotropy predominantly influences the pressure field, whereas relative permeability anisotropy exerts a greater effect on oil–water flow dynamics and residual oil distribution.

The previous analyses examined the individual effects of absolute and relative permeability anisotropy. In actual reservoirs, both anisotropies typically coexist. Figure 17 and Figure 22 present the pressure and oil saturation fields under the combined influence of absolute and relative permeability anisotropies when both the dominant permeability and the higher water-phase relative permeability are aligned in the x-direction. In contrast, Figure 18 and Figure 23 depict these fields when the dominant permeability is oriented in the y-direction while the higher water-phase relative permeability remains in the x-direction.

Figure 17 shows a pressure field similar to that in Figure 15, indicating that incorporating relative permeability anisotropy does not substantially alter the pressure field compared to considering absolute permeability anisotropy alone. Figure 18, which represents a 90° rotation of Figure 15, confirms that the effect of relative permeability anisotropy on the pressure field is minimal compared to the impact of absolute permeability anisotropy. These observations support the conclusion that absolute permeability anisotropy primarily influences the pressure field, with relative permeability anisotropy having a lesser effect.

The trends in the saturation fields shown in Figure 20 and Figure 22 are similar, with residual oil concentrating predominantly in the y-direction. However, Figure 22 demonstrates that when the effects of anisotropic absolute permeability and anisotropic relative permeability are combined in the same direction, the saturation field and residual oil enrichment are more pronounced in that dominant direction. In contrast, Figure 23 shows that the saturation field and residual oil distribution become asymmetric when anisotropies with different dominant directions are combined. During the early stages of water injection, the saturation field shifts notably in the y-direction. Over time, this shift lessens as anisotropic relative permeability enhances flooding in the x-direction, thereby reducing the uneven distribution of residual oil saturation.

In summary, anisotropic relative and absolute permeability have distinct effects on multi-phase flow and cannot be considered equivalent or interchangeable. In practical reservoir applications, it is crucial to comprehensively evaluate these anisotropies’ combined impact on water injection development.

Water cut curves (Figure 24) and cumulative oil production curves (Figure 25) were compared to quantitatively analyze the impact of anisotropy on the time of water breakthrough, final water cut, and yields. Figure 24 shows that water breakthrough occurs in Case 1, then Case 2, Case 3, and finally Case 4 and Case 5. Figure 25 indicates that Case 1 yields the lowest production, followed by Case 2, while Cases 3, 4, and 5 have similar, higher yields. This result contrasts with the expectation that homogeneous reservoirs should exhibit the most complete displacement, the slowest water flooding, and the most favorable production outcomes.

The observed results arise from the simplicity of the injection and production model used in the study. The wells are positioned diagonally in this model. In Case 1, the homogeneous reservoir allows for uniform water drive, leading to the fastest water flooding along the diagonal and the lowest production. In Case 2, where Kx > Ky, the injected water floods rapidly in the x-direction. This results in reduced effectiveness of flooding along the diagonal, causing a delay in water breakthrough and an increase in the pre-breakthrough flooding area, which extends the production period. Despite this, once water breakthrough occurs, the water cut increases rapidly, leading to a higher final water cut than the homogeneous reservoir, with less favorable production outcomes. This suggests that permeability anisotropy forms dominant flow channels, hindering efficient waterflooding. In Case 3 to Case 5, the directionality of absolute or relative permeability results in extended flooding times along the diagonal. These cases exhibit distinct variations in water breakthrough timing and the rate of water cut increase. Although the final water cut is lower and production higher in these cases compared to Case 1, these outcomes are likely due to the model’s simplicity and well network configuration. In real reservoirs, the trends may vary. Nevertheless, the observed qualitative trends argue that anisotropy relative permeability significantly impacts final yield.

4. Conclusions

This study introduces a multi-scale numerical simulation method for oil–water two-phase flow accounting for relative permeability anisotropy. The method’s accuracy and performance are validated using the Egg model. Furthermore, a typical model is built to analyze the influence of absolute and relative permeability anisotropy on waterflooding performance.

(1) The anisotropy of relative permeability has been proved, yet most existing reservoir simulators neglect its impact. This study presents a numerical simulation method, developed using MRST and based on the black-oil model, which explicitly incorporates the effects of anisotropic relative permeability;

(2) Following the development of the multi-scale model, this paper adopts a sequential solution method to enhance computational efficiency. The method leverages the parallelism of the MsFV method to solve the pressure equation. The Egg model’s numerical validation confirms that the proposed IMsFV method significantly improves the iterative speed while preserving result accuracy;

(3) The impact of anisotropy on waterflooding is investigated using the developed multi-scale numerical simulator. Findings indicate that anisotropic absolute and relative permeability play critical roles in waterflooding, though they operate via distinct mechanisms and cannot be equated, substituted, or entirely counterbalanced;

(4) Compared to the traditional numerical simulator, the IMsFV method proposed accurately captures the anisotropy of relative permeability in reservoirs. This advancement is crucial for high-resolution numerical simulation of fluvial-phase sedimentary reservoirs, improved characterization of flow regimes, and more precise prediction of residual oil distribution.

Despite its significant applications, this study identifies specific limitations that will be addressed in future research:

(1)This paper’s developed multi-scale reservoir numerical simulation method is limited to two-phase (oil and water) black-oil models. Future work aims to extend the model to encompass three-phase black-oil and component models, facilitating faster multi-scale solutions;

(2)The analysis of relative permeability anisotropy using the established multi-scale numerical simulator employs a simplified model setup. Future research will enhance the application of this simulation method in real reservoirs.

Author Contributions

Conceptualization, L.W. and S.W.; methodology, J.W.; software, J.Z.; validation, L.W., J.W., and S.W.; formal analysis, R.Z.; investigation, Y.Y.; resources, R.Z.; data curation, Y.Y.; writing—original draft preparation, L.W.; writing—review and editing, J.W.; visualization, J.Z.; supervision, S.W.; project administration, D.J.; funding acquisition, D.J. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

Deli Jia, Jiqun Zhang and Yiqun Yan were employed by PetroChina. Ruichao Zhang was employed by Shandong Institute of Petroleum and Chemical Technology. All authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflicts of interest.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Acknowledgments

The authors would like to show sincere thanks to those who have contributed to this research.

Glossary

Nomenclature

Physical Quantities ? Porosity, fraction s a Saturation of phase a, volume fraction b a Inverse formation value factor of phase, dimensionless ? ? a Darcy velocities of phase a, cm/s p a Pressure of phase a, bar k r a Relative permeability of phase a, dimensionless k r a ( x ) Relative permeability of phase a in the x-direction K Absolute permeability, mD µ a Viscosity of phase a, kg/ms g Gravity acceleration, kg.m/s2 z Vertical depth, m q o, q w Volumetric production/injection rate, s–1 ? a Density of phase a, kg/m3 p c Capillary pressure, bar ? a i Mobility of phase a, mD kg/ms f w Water cut, fraction f a Partial flow coefficient, fraction Sor Residual oil saturation, fraction Sor(x) Residual oil saturation in the x-direction, fraction W I i Well indices, dimensionless twb Time of water breakthrough, year Domain and Grid B [ C ( f ) ] Coarse block containing cell c containing face f c Cell number { O i f } i = 1 n f The fine-scale grid domain v { O j c } j = 1 n c The coarse-scale grid domain Discretization and Discrete Quantities R w Phase pressure equations in residual form O e Pressure equations in residual form R q, a Well production equations in residual form R C Well control equations in residual form Linear Algebra A, b, x (system) matrix, left-hand side, and vector of unknowns D Diagonal matrix J e / J d J Jacobian matrix of discrete equations P prolongation operator (matrix) R Restriction operator (matrix) r Vector of residuals e Vector with all unknowns in pressure equation d Vector of unknown quantities for the transport equation e p / e q / e C Iterative convergence constant for the pressure equation e M / e V Iterative convergence constant for the transport equation e i t External iterative convergence constant Subscripts a Fluid phase (a=w,o) b Basis function o Oil w Water i, j, k Cell numbers bh Bottomhole pressure x, y, z x-, y-, z-directions f Fine grid c Coarse grid ref Reference value Superscripts n Timestep b Basis function 0 Initial moment s Standard conditions

References

1. J. Farquharson; F.B. Wadsworth Upscaling permeability in anisotropic volcanic systems., 2018, 364,pp. 35-47. DOI: https://doi.org/10.1016/j.jvolgeores.2018.09.002.

2. T. Li; M. Li; X. Jing; W. Xiao; Q. Cui Influence mechanism of pore-scale anisotropy and pore distribution heterogeneity on permeability of porous media., 2019, 46,pp. 594-604. DOI: https://doi.org/10.1016/S1876-3804(19)60039-X.

3. A.T. Corey; C.H. Rathjens Effect of Stratification on Relative Permeability., 1956, 8,pp. 69-71. DOI: https://doi.org/10.2118/744-G.

4. B. Busahmin; R.R. Karri; S. Tyson; M. Jami Modeling of a long sand-pack for heavy crude oil through depletion tests utilizing methane gas., 2021, 7,pp. 188-198. DOI: https://doi.org/10.1016/j.petlm.2020.07.002.

5. Crotti; J.A. Rosbaco Relative Permeability Curves: The Influence of Flow Direction and Heterogeneities.,

6. C.C. Ezeuko; S.R. McDougall; I. Bondino; G. Hamon Anisotropic Relative Permeabilities for Characterising Heavy-Oil Depletion Experiment.,

7. L. Paterson; S. Painter; X. Zhang; W.V. Pinczewski Simulating Residual Saturation and Relative Permeability in Heterogeneous Formations., 1998, 3,pp. 211-218. DOI: https://doi.org/10.2118/50938-PA.

8. S. Bakhshian; S.A. Hosseini; L. Lake CO2-brine relative permeability and capillary pressure of Tuscaloosa sandstone: Effect of anisotropy., 2020, 135,p. 103464. DOI: https://doi.org/10.1016/j.advwatres.2019.103464.

9. Z. Lei; T. Liu; C. Xie; M. Wang; Z.-K. Zhang Predictions of Relative Permeability for Low Permeability Reservoirs and its Scale Effect.,

10. A.K. Pergament; P.Y. Tomin The study of relative phase-permeability functions for anisotropic media., 2012, 4,pp. 1-9. DOI: https://doi.org/10.1134/S2070048212010097.

11. X. Pei; Y. Liu; W. Gu; C. Jian Characterization and Analysis of Anisotropic Relative Permeability., 2021, 27,pp. 579-596. DOI: https://doi.org/10.2118/206724-PA.

12. N. Raoult; S. Beylat; J.M. Salter; F. Hourdin; V. Bastrikov; C. Ottlé; P. Peylin Exploring the Potential of History Matching for Land Surface Model Calibration., 2024, 2024,pp. 1-33. DOI: https://doi.org/10.5194/gmd-17-5779-2024.

13. C. Li; S. Wang; Q. You; C. Yu A New Measurement of Anisotropic Relative Permeability and Its Application in Numerical Simulation., 2021, 14, 4731. DOI: https://doi.org/10.3390/en14164731.

14. M.H. Sedaghat; S. Azizmohammadi; S.K. Matthäi Tensor Analysis of the Relative Permeability in Naturally Fractured Reservoirs., 2020, 25,pp. 162-184. DOI: https://doi.org/10.2118/197064-PA.

15. A.V. Blonsky; D.A. Mitrushkin; I.Y. Kudryashov; V.V. Plynin Computation of Absolute and Relative Permeability Full Tensors for Fractured Reservoirs.,

16. X. Pei; Y. Liu; Z. Lin; Y. Mao; L. Xue Anisotropic Characteristics of Relative Permeability in Sedimentary Reservoirs.,

17. R. Li; Z. Han; L. Zhang; J. Zhou; S. Wang; F. Huang Numerical Determination of Anisotropic Permeability for Unconsolidated Hydrate Reservoir: A DEM–CFD Coupling Method., 2024, 12, 1447. DOI: https://doi.org/10.3390/jmse12081447.

18. E. Keilegavlen; J.M. Nordbotten; A.F. Stephansen Simulating Two-phase Flow in Porous Media with Anisotropic Relative Permeabilities.,

19. N.C. Kokkinos; D.C. Nkagbu; D.I. Marmanis; K. Dermentzis; G. Maliaris Evolution of Unconventional Hydrocarbons: Past, Present, Future and Environmental FootPrint., 2022, 15,pp. 15-24. DOI: https://doi.org/10.25103/jestr.154.03.

20. H. Wang; F. Ma; X. Tong; Z. Liu; X. Zhang; Z. Wu; D. Li; B. Wang; Y. Xie; L. Yang Assessment of global unconventional oil and gas resources., 2016, 43,pp. 925-940. DOI: https://doi.org/10.1016/S1876-3804(16)30111-2.

21. C. Zou; G. Zhai; G. Zhang; H. Wang; G. Zhang; J. Li; Z. Wang; Z. Wen; F. Ma; Y. Liang et al. Formation, distribution, potential and prediction of global conventional and unconventional hydrocarbon resources., 2015, 42,pp. 14-28. DOI: https://doi.org/10.1016/S1876-3804(15)60002-7.

22. O.A. Hisseine; A. Tagnit-Hamou Development of ecological strain-hardening cementitious composites incorporating high-volume ground-glass pozzolans., 2020, 238,p. 117740. DOI: https://doi.org/10.1016/j.conbuildmat.2019.117740.

23. C. Zou; G. Zhang; S. Tao; S. Hu; X. Li; J. Li; D. Dong; R. Zhu; X. Yuan; L. Hou et al. Geological features, major discoveries and unconventional petroleum geology in the global petroleum exploration., 2010, 37,pp. 129-145. DOI: https://doi.org/10.1016/S1876-3804(10)60021-3.

24. C. Jia; M. Zheng; Y. Zhang Unconventional hydrocarbon resources in China and the prospect of exploration and development., 2012, 39,pp. 139-146. DOI: https://doi.org/10.1016/S1876-3804(12)60026-3.

25. S. Mao; K. Wu; G. Moridis Integrated simulation of three-dimensional hydraulic fracture propagation and Lagrangian proppant transport in multilayered reservoirs., 2023, 410,p. 116037. DOI: https://doi.org/10.1016/j.cma.2023.116037.

26. S. Krogstad; K.A. Lie; O. Møyner; H.M. Nilsen; X. Raynaud; B. Skaflestad MRST-AD—An Open-Source Framework for Rapid Prototyping and Evaluation of Reservoir Simulation Problems.,

27. K.A. Lie; S. Krogstad; I.S. Ligaarden; J.R. Natvig; H.M. Nilsen; B. Skaflestad Open-source MATLAB implementation of consistent discretisations on complex grids., 2011, 16,pp. 297-322. DOI: https://doi.org/10.1007/s10596-011-9244-4.

28. D.W. Peaceman Interpretation of Well-Block Pressures in Numerical Reservoir Simulation with Nonsquare Grid Blocks and Anisotropic Permeability., 1983, 23,pp. 531-543. DOI: https://doi.org/10.2118/10528-PA.

29. J.R. Wallis Incomplete Gaussian Elimination as a Preconditioning for Generalized Conjugate Gradient Acceleration.,

30. S. Gries; K. Stüben; G.L. Brown; D. Chen; D. Collins Preconditioning for Efficiently Applying Algebraic Multigrid in Fully Implicit Reservoir Simulations., 2014, 19,pp. 726-736. DOI: https://doi.org/10.2118/163608-PA.

31. O. Møyner; K.-A. Lie A Multiscale Restriction-Smoothed Basis Method for Compressible Black-Oil Models., 2016, 21,pp. 2079-2096. DOI: https://doi.org/10.2118/173265-PA.

32. M. Cusini; A.A. Lukyanov; J.R. Natvig; H. Hajibeygi Constrained pressure residual multiscale (CPR-MS) method for fully implicit simulation of multiphase flow in porous media., 2015, 299,pp. 472-486. DOI: https://doi.org/10.1016/j.jcp.2015.07.019.

33. A. Brandt Algebraic multigrid theory: The symmetric case., 1986, 19,pp. 23-56. DOI: https://doi.org/10.1016/0096-3003(86)90095-0.

34. O. Møyner; K.A. Lie The Multiscale Finite-Volume Method on Stratigraphic Grids., 2014, 19,pp. 816-831. DOI: https://doi.org/10.2118/163649-PA.

35. O. Møyner; K.A. Lie The Multiscale Finite Volume Method on Unstructured Grids.,

36. A. Kozlova; Z. Li; J.R. Natvig; S. Watanabe; Y. Zhou; K. Bratvedt; S.H. Lee A Real-Field Multiscale Black-Oil Reservoir Simulator., 2016, 21,pp. 2049-2061. DOI: https://doi.org/10.2118/173226-PA.

37. I. Lunati; S. Lee An Operator Formulation of the Multiscale Finite-Volume Method with Correction Function., 2009, 8,pp. 96-109. DOI: https://doi.org/10.1137/080742117.

38. M. Mosharaf-Dehkordi; M.T. Manzari Effects of using altered coarse grids on the implementation and computational cost of the multiscale finite volume method., 2013, 59,pp. 221-237. DOI: https://doi.org/10.1016/j.advwatres.2013.07.003.

39. O. Møyner; K.A. Lie A multiscale two-point flux-approximation method., 2014, 275,pp. 273-293. DOI: https://doi.org/10.1016/j.jcp.2014.07.003.

40. P. Jenny; S.H. Lee; H.A. Tchelepi Multi-scale finite-volume method for elliptic problems in subsurface flow simulation., 2003, 187,pp. 47-67. DOI: https://doi.org/10.1016/S0021-9991(03)00075-5.

41. Y. Wang; H. Hajibeygi; H.A. Tchelepi Algebraic multiscale solver for flow in heterogeneous porous media., 2014, 259,pp. 284-303. DOI: https://doi.org/10.1016/j.jcp.2013.11.024.

42. O. Møyner; H. Tchelepi A multiscale restriction-smoothed basis method for high contrast porous media represented on unstructured grids., 2016, 304,pp. 46-71. DOI: https://doi.org/10.1016/j.jcp.2015.10.010.

43. O. Nwanwe; N.C. Izuwa; N.P. Ohia; A. Kerunwa; U. Nnaemeka Determining optimal controls placed on injection/production wells during waterflooding in heterogeneous oil reservoirs using artificial neural network models and multi-objective genetic algorithm., 2024, DOI: https://doi.org/10.1007/s10596-024-10300-2.

44. X. Tian; O. Volkov; D. Voskov An advanced inverse modeling framework for efficient and flexible adjoint-based history matching of geothermal fields., 2024, 116,p. 102849. DOI: https://doi.org/10.1016/j.geothermics.2023.102849.

45. P. Kor; A. Hong; R.B. Bratvold Reservoir Production Management with Bayesian Optimization: Achieving Robust Results in a Fraction of the Time., 2023,pp. 1-21. DOI: https://doi.org/10.2118/217985-PA.

46. C.S. Lee; F.P. Hamon; N. Castelletto; P.S. Vassilevski; J.A. White Multilevel well modeling in aggregation-based nonlinear multigrid for multiphase flow in porous media., 2024, 513,p. 113163. DOI: https://doi.org/10.1016/j.jcp.2024.113163.

47. J.D. Jansen; R.M. Fonseca; S. Kahrobaei; M.M. Siraj; v.G. Essen; P.M. Van den Hof The egg model—A geological ensemble for reservoir simulation., 2014, 1,pp. 192-195. DOI: https://doi.org/10.1002/gdj3.21.

48. H. Zhou; S.H. Lee; H.A. Tchelepi Multiscale Finite-Volume Formulation for the Saturation Equations., 2011, 17,pp. 198-211. DOI: https://doi.org/10.2118/119183-PA.

49. O. Chaabi; M.A. Kobaisi Algorithmic Monotone Multiscale Finite Volume Methods for Porous Media Flow., 2024, 499,p. 112739. DOI: https://doi.org/10.1016/j.jcp.2023.112739.

50. M. Cusini; C. van Kruijsdijk; H. Hajibeygi Algebraic dynamic multilevel (ADM) method for fully implicit simulations of multiphase flow in porous media., 2016, 314,pp. 60-79. DOI: https://doi.org/10.1016/j.jcp.2016.03.007.

Figures and Tables

Figure 1: Flow chart of the technology roadmap. [Please download the PDF to view the image]

Figure 2: Measured relative permeability in x-, y-, and z-directions from a core sample in a lateral bar sedimentary reservoir. [Please download the PDF to view the image]

Figure 3: The schematic diagram used to determine the direction of the grid surface. [Please download the PDF to view the image]

Figure 4: Diagram describing the multi-scale grid (small rectangles of different colors represent fine grids, while black rectangles represent coarse grids. Different colors indicate varying permeability values). [Please download the PDF to view the image]

Figure 5: Reservoir displaying the permeability distribution (shown in red-blue gradient colors), as well as the position of the injectors (blue cylinders) and producers (red cylinders) of the Egg Model (Refer to J. D. Jansen, 2014 [47]). [Please download the PDF to view the image]

Figure 6: Traditional fully implicit solution of pressure at different times for the Egg model. [Please download the PDF to view the image]

Figure 7: Multi-scale solution of pressure at different times for the Egg model. [Please download the PDF to view the image]

Figure 8: Traditional fully implicit solution of saturation at different times for the Egg model. [Please download the PDF to view the image]

Figure 9: Multi-scale solution of saturation at different times for the Egg model. [Please download the PDF to view the image]

Figure 10: Comparison between oil production versus time. [Please download the PDF to view the image]

Figure 11: Comparison between water-cut versus time. [Please download the PDF to view the image]

Figure 12: The comparison of the computation time of each step. [Please download the PDF to view the image]

Figure 13: Anisotropic relative permeability curves in the model. [Please download the PDF to view the image]

Figure 14: Multi-scale solution of pressure at different times for Case 1 (P1 represents the oil production well, and I1 represents the water injection well, the same applies figures below). [Please download the PDF to view the image]

Figure 15: Multi-scale solution of pressure at different times for Case 2. [Please download the PDF to view the image]

Figure 16: Multi-scale solution of pressure at different times for Case 3. [Please download the PDF to view the image]

Figure 17: Multi-scale solution of pressure at different times for Case 4. [Please download the PDF to view the image]

Figure 18: Multi-scale solution of pressure at different times for Case 5. [Please download the PDF to view the image]

Figure 19: Multi-scale solution of oil saturation at different times for Case 1. [Please download the PDF to view the image]

Figure 20: Multi-scale solution of oil saturation at different times for Case 2. [Please download the PDF to view the image]

Figure 21: Multi-scale solution of oil saturation at different times for Case 3. [Please download the PDF to view the image]

Figure 22: Multi-scale solution of oil saturation at different times for Case 4. [Please download the PDF to view the image]

Figure 23: Multi-scale solution of oil saturation at different times for Case 5. [Please download the PDF to view the image]

Figure 24: Water cut curves for different cases. [Please download the PDF to view the image]

Figure 25: Cumulative oil production curves for different cases. [Please download the PDF to view the image]

Table 1: Basic parameters of a cubic core.

Table 2: The comparison of the total computation time by traditional fully implicit and multi-scale methods.

Table 3: Simulation model parameters of the base case.

Table 4: Simulation cases.

Table 5: Endpoints of anisotropic relative permeability.

Author Affiliation(s):

[1] School of Energy, China University of Geosciences (Beijing), Beijing 100083, China; [emailprotected] (L.W.); [emailprotected] (J.W.); [emailprotected] (R.Z.)

[2] Research Institute of Petroleum Exploration & Development, PetroChina, Beijing 100083, China; [emailprotected] (D.J.); [emailprotected] (J.Z.); [emailprotected] (Y.Y.)

[3] Shandong Institute of Petroleum and Chemical Technology, Dongying 257061, China

Author Note(s):

[*] Correspondence: [emailprotected]; Tel.: +86-133-7013-8719

DOI: 10.3390/pr12092058