Implicit Pressure Explicit Saturation Research Articles

Physics-informed neural network simulation of two-phase flow in heterogeneous and fractured porous media

Physics-informed neural networks (PINNs) have received great attention as a promising paradigm for forward, inverse, and surrogate modeling of various physical processes with limited or no labeled data. However, PINNs are rarely used to predict two-phase flow in heterogeneous and fractured porous media, which is critical to lots of subsurface applications, due to the significant challenges in their training. In this work, we present an Enriched Physics-Informed Neural Network (E-PINN) to overcome these barriers and realize the simulation of such flow. Specifically, the Embedded Discrete Fracture Model (EDFM) is adopted to explicitly represent fractures, and then the finite volume method (FVM) instead of the Automatic Differentiation (AD) is used to evaluate spatial derivatives and construct the physics-informed loss function, so that the flux continuity between neighboring elements with different properties (e.g. matrix and fracture) can be defined rigorously. Besides, we develop a novel physics-informed neural network (NN) architecture adopting the adjacency-location anchoring, adaptive activation function, skip connection and gated updating to enrich the pressure information and enhance the learning ability of NN. Additionally, the initial and boundary conditions are constrained through a hard approach, which encodes them into network design, to improve the accuracy and efficiency of network training. In order to further reduce the difficulty of training, the Implicit-Pressure Explicit-Saturation (IMPES) scheme is used to calculate pressure and saturation, in which only the pressure needs to be solved by training NN. Finally, the superiority and applicability of E-PINN to complex practical problems is demonstrated through the simulations of immiscible displacement in 2D/3D heterogeneous and fractured reservoirs.

A novel streamline simulation method for fractured reservoirs with full-tensor permeability

In this work, we develop a novel streamline (SL) simulation method that integrates seamlessly within the embedded discrete fracture model (EDFM). The novel SL-based method is developed based on a hybrid of two-point flux approximation (TPFA) and mimetic finite difference (MFD) methods, which is applicable to a two-phase anisotropic flow in fractured reservoirs. We refer to this novel approach as EDFM-TPFA-MFD-SL. The approach is operated in an IMplicit Pressure Explicit Saturation (IMPES) manner. First, this work establishes a novel EDFM utilizing a hybrid TPFA-MFD scheme to solve the pressure equation for phase flux approximation. Subsequently, we introduce a practical streamline tracing workflow designed to derive the distribution of streamlines within the reservoir domain and the time-of-flight distribution along each individual streamline. This feature allows for the parallel computation of water saturation along the streamlines. Two numerical examples are presented to validate the superiority of the proposed EDFM-TPFA-MFD-SL method compared to the existing streamline-based EDFM on cases with full-tensor permeability. The results show that the proposed method could significantly mitigate the numerical dissipation and reduce the computation costs. Another numerical example demonstrates the effectiveness of the proposed method in dealing with complex fracture networks and providing rapid flow diagnostics, indicating its significant potential for real-world field applications.

Numerical investigation of gas and sand production from hydrate-bearing sediments by incorporating sand migration based on IMPES method

Accurately predicting gas and sand production from oceanic gas hydrate reservoirs considering multi-field coupling effects is crucial. To address this, we first systematically summarized the multi-field coupling effect of gas hydrate-bearing sediments (GHBS) during sand production and gas production. After that, we developed a thermo-hydro-mechanical coupled model that incorporates both sand detachment and sand migration (THMS) to describe the complex behaviors of fluid and solid flow. We then solved this model using the implicit pressure-explicit saturation (IMPES) method and incremental scheme for a cylindrical core and presented a detailed algorithm derivation process for reproducibility. The newly developed model was validated through comparison with available experimental data. Subsequently, several parameter effects, including outlet pressure, initial hydrate saturation, load stress, and sand detachable ratio, on the behaviors of fluid and solid flow were studied. The computation results indicate that high sand production can be induced by low outlet pressure, high initial hydrate saturation, high load stress, and high sand detachable ratio. Meanwhile, high cumulative gas production can be caused by low outlet pressure, low initial hydrate saturation, high load stress, and high sand detachable ratio.

Discontinuous Galerkin method for hybrid-dimensional fracture models of two-phase flow

In this paper, we present a discontinuous Galerkin (DG) method for the simulation of incompressible and immiscible two-phase flow in fractured porous media with capillary pressure. The hybrid-dimensional fracture models are considered in which fractures are treated as (d−1)-dimensional interfaces in the d-dimensional domain and fluid exchange between fractures and surrounding matrix is taken into account. The approximation schemes combining the DG method with the fully implicit backward Euler time discretization are formulated for models with a single fracture and intersecting fractures network, and the stability and unique solvability are established. Optimal error estimates in the discrete H1-norm for saturation and pressure are derived. Numerical experiments in two-dimensional geometry are conducted to verify the accuracy of our theoretical analysis. In addition to the backward Euler formulation, we also implement the algorithm with a kind of conventional IMplicit Pressure Explicit Saturation (IMPES) method.

A locally conservative multiscale method for stochastic highly heterogeneous flow

In this paper, we propose a local model reduction approach for subsurface flow problems in stochastic and highly heterogeneous media. Generally, we apply the mixed generalized multiscale finite element method (MGMsFEM) to solve the single-phase flow equation. As for the two-phase flow, we use mixed finite element method (MFEM) with piecewise constant bases on a fine mesh to deal with the transport equation and a standard Implicit Pressure Explicit Saturation (IMPES) scheme is utilized in the time-marching process. To guarantee the mass conservation, we consider the mixed formulation of the flow problem and aim to solve the problem in a coarse grid to reduce the complexity of a large-scale system. We decompose the entire problem into a training and a testing stage, namely the offline coarse-grid multiscale basis generation stage and online simulation stage with different parameters. In the training stage, a parameter-independent and small-dimensional multiscale basis function space is constructed, which includes the media, source and boundary information. The key part of the basis generation stage is to solve some special local problems. With the parameter-independent basis space, one can efficiently solve the concerned problems corresponding to different samples of permeability field in a coarse grid without repeatedly constructing a multiscale space for each new sample. A rigorous analysis on convergence of the proposed method is presented. In particular, we consider a generalization error, where bases constructed with one source will be used to a different source. In the numerical experiments, we apply the proposed method for both single-phase and two-phase flow problems. Simulation results for both 2D and 3D representative models demonstrate the high accuracy and efficiency of the proposed model reduction techniques.

Numerical solution of extended black-oil model incorporating capillary effects based on a high-resolution central scheme

The purpose of this paper is to present an efficient numerical method to solve the extended black-oil model’s equations with considering the capillary pressure effects. The governing equations are formulated based on a global pressure concept to account for the capillary effects. The developed equations are conservative and suitable for the control volume approach. To solve the equations, an implicit pressure-explicit saturation scheme is used. The component transport equations are decomposed to a pure convection term driven by the global pressure and a pure diffusion part affected by the capillary pressures. A high-resolution central scheme is employed to discretize the convection term of the component transport equations. The results of the proposed model are compared with several semi-analytical and numerical solutions that were previously published. It is shown that the suggested model is computationally efficient and easy to be implemented for solving convection-diffusion problems. The results illustrate reasonable accuracy in comparison with the reference solutions. As a result of simplicity, low computational cost, and decent accuracy, this method is appropriate for solving complicated problems of reservoir engineering.

A physics-informed convolutional neural network for the simulation and prediction of two-phase Darcy flows in heterogeneous porous media

The physics-informed neural network (PINN) is a general deep learning framework for simulating flows with limited or no labeled data. In the current study, we develop a physics-informed convolutional neural network (PICNN) for simulating transient two-phase Darcy flows in heterogeneous reservoir models with source/sink terms in the absence of labeled data, where the finite volume method (FVM) is adopted to approximate the PDE residual in the loss function such that flux continuity between neighboring cells of different properties is defined rigorously, and a well model is adopted to approximate the high pressure gradient near sources or sinks. The implicit-pressure explicit-saturation (IMPES) scheme is employed such that only a single CNN needs to be trained per time step. Dirichlet boundary conditions are not a mandatory requirement for PICNN-based implicit solver but act as labeled data that can help enhance accuracy. The proposed approach is validated in homogeneous and heterogeneous reservoirs and aspects including efficiency and accuracy are discussed. In addition, we demonstrate that the CNN structure can be trained as a data-driven surrogate for two-phase Darcy flows given sufficient labeled samples.

Numerical simulation of oil-water displacements in naturally fractured reservoirs using a hybrid grid MPFA-streamline (SLS-HyG) formulation

The numerical simulation of fluid displacements in naturally fractured reservoirs is rather expensive in terms of CPU requirements and challenging, due to the random spatial distribution of fractures of different sizes and shapes. In this paper, we present a numerical formulation to simulate the 2-D two-phase flow of oil and water in naturally fractured oil reservoirs using unstructured quadrilateral meshes. In our Hybrid-Grid MPFA-Streamline (SLS-HyG) formulation, the elliptic pressure equation is discretised by a finite volume method with a Multipoint Flux Approximation using a Diamond stencil. This formulation is able to handle a full permeability tensor representation in an arbitrary polygonal grid. For the advective saturation problem, we used the streamline-based method, which decouples the transport equations into several 1-D problems solved along each streamline. These formulations were implemented in the framework of the Hybrid-Grid (HyG) representation of fractures, where the fractures in the physical domain are represented by 1-D objects in 2-D domains and extended to the same dimension of the matrix problem. The coupling between the flow and the transport problem is done by a modified implicit pressure explicit saturation procedure. We have solved some reference problems from the literature, obtaining promising results when compared to other classical formulations.

Boltzmann transformation of radial two-phase black oil model for tight oil reservoirs

Tight oil accumulates in impermeable reservoir rocks, often shale or tight sandstones. The flow behaviour of tight oil in unconventional reservoirs is described by peculiar complexities such as the typical low permeability to viscosity ratio and the dissolution of some gases in the oil phase. Reservoir simulations that consider these complexities negligible stand the potential of poorly characterizing the reservoir flow dynamics. The adoption of similarity transformation effectively reduces the complexities associated with the flow equations through spatial variable (r) and temporal variable (t). The Boltzmann variable left(xi =dfrac{r}{sqrt{t}}right) is introduced to facilitate the reformulation of transient two-phase flow phenomenon in a radial geometry. The technique converts the two-phase Black oil model (thus highly nonlinear partial differential equations (PDEs)) to ordinary differential equations (ODEs). The resulting ODEs present a reduced form on the flow model which is solved by finite difference approximations (the Implicit-Pressure-Explicit-Saturation (IMPES)) scheme. No loss of vital flow characteristics was observed between the Black oil model and the similarity transform flow model. Furthermore, the similarity approach facilitated the determination of pressure and saturation equations as unique functions of the Boltzmann variable. This derivation is applied to an infinitely acting reservoir where the Boltzmann variable tends to infinity (xi rightarrow infty). Finally, this case study’s analytical solution formulated critical relations for fluid flow rate and cumulative production, which are useful for single-phase flow analysis.

A very high-order flux reconstruction approach coupled to the MPFA-QL finite volume method for the numerical simulation of oil-water flows in 2D petroleum reservoirs

In recent years, high-order methods for the discretization of the advection dominated saturation equation, and that results from the modelling of the 2D Oil-Water displacement through porous formations, have received renewed attention from the CFD community, as a suitable alternative to the traditional first-order upwind (FOU) finite volume (FV) scheme, which is the standard in industry. A key aspect of high-order methods applied to numerical simulation of petroleum reservoirs is the lower computational cost for a comparable degree of accuracy against their lower-order counterparts. In this context, in the present paper, for the first time in literature, we propose the use of the very high order Flux Reconstruction (FR) formulation for the discretization of the advective terms in the saturation equation using unstructured quadrilateral meshes. To solve the elliptic pressure equation, we use the recently developed Multipoint Flux Approximation method with a Quasi-Local Stencil (MPFA-QL). This method is very robust and able to handle highly heterogeneous and anisotropic permeability tensors. The system of equations is solved in a segregated manner by using the classical Implicit Pressure Explicit Saturation (IMPES) procedure. To properly couple the cell-centered MPFA-QL method with the FR formulation, it is necessary to obtain an adequate velocity reconstruction throughout the control volumes of the mesh. As the cell-centered finite volume method naturally delivers fluxes across control surfaces of the primal grid, then, we use a reconstruction operator based on the first-order Brezzi–Douglas–Marini (BDM) space, and the Piola transformation, to get the complete knowledge of the velocity field throughout the domain. To eliminate spurious oscillations around the discontinuities, which are typical in high-order methods, we have used a hierarchical multidimensional limiting process (hMLP). To compute fluxes over the control surfaces we have used the Roe-E approximate Riemann solver. Finally, through some representative examples, we show the accuracy, efficiency and robustness of our new high order numerical methodology to model oil-water displacements in highly heterogeneous and anisotropic petroleum reservoirs.

A multiphase multicomponent flow and transport model for liquid aerosol filtration in coalescing fibrous filters

Pore-scale CFD simulations of liquid aerosol filtration in coalescing fibrous filters are currently limited by severe computational constraints. In this paper, we present a new multiphase multicomponent framework for modeling such dynamic filtration process at the Darcy-scale. Moving from the pore- to the Darcy-scale allows for significant reduction in computational efforts with limited repercussion on the representation of the filter geometry. In our model, two phases coexist within the filter, namely the injected oil-mist and the coalesced oil phase. The oil-mist phase consists of a mixture of a gas phase and oil droplets of different sizes. The governing equations are the standard ones of Darcy flow and transport in porous media, where the filtration processes, such as droplet capture and drainage, are accounted for through ad hoc mass exchange terms, which are coherent with the Jump & Channel model by Kampa et al. (2014). The equations are strongly non-linear and are solved through the IMplicit Pressure Explicit Saturation (IMPES) algorithm. Cell-centered finite volumes are used for discretization in space. The model is validated by replicating experiments of oil-mist filtration available in the literature for both wettable and non-wettable media, as well as for combinations of media with different wettability. Simulations are then performed for a variety of operational conditions, demonstrating both the accuracy and robustness of our implementation.

An adaptive time stepping algorithm for IMPES with high order polynomial extrapolation

Two-phase flows in oil reservoirs can be modeled by a coupled system of elliptic and hyperbolic partial differential equations. The transport velocity of the multiphase fluid system is related to the pressure through Darcy’s law and it is coupled to a conservation law for the saturation variable of one of the phases. A time step of the classical IMPES (IMplicit Pressure Explicit Saturation) method consists of first solving the elliptic problem for pressure and Darcy velocity, and then updating the saturation with an explicit numerical scheme for conservation laws. This method is very computationally costly, since the time-consuming elliptic solver must be invoked at time intervals defined by the stability limit of the hyperbolic solver. A popular variant is not to update the velocity at all hyperbolic time steps, but to skip a fixed number C of them, with C determined by the user. In this work we propose a more accurate and systematic procedure for time stepping in IMPES codes. The velocity is updated at all transport time steps, though the elliptic solver is only invoked every C steps. In the time steps at which the elliptic problem is not solved, the velocity is extrapolated from previously computed values with polynomials of high degree. Further, we introduce an error estimator that allows for the number C to be adaptively determined without user intervention. The algorithm was tested in several relevant benchmark problems. This allowed for the optimization of its parameters and comparisons with previous variants. The results show that the proposed algorithm is very stable, reliable and time-cost effective. It is also easily implemented in pre-existent IMPES codes.

Fully Implicit Dynamic Pore‐Network Modeling of Two‐Phase Flow and Phase Change in Porous Media

AbstractDynamic pore‐network model (PNM) has been widely used to model pore‐scale two‐phase flow. Numerical algorithms commonly used for dynamic PNM including IMPES (implicit pressure explicit saturation) and IMP‐SIMS (implicit pressure semi‐implicit saturation) can be numerically unstable or inaccurate for challenging flow regimes such as low capillary number (Ca) flow and unfavorable displacements. We perform comprehensive analyses of IMPES and IMP‐SIMS for a wide range of flow regimes under drainage conditions and develop a novel fully implicit (FI) algorithm to address their limitations. Our simulations show the following: (1) While IMPES was reported to be numerically unstable for low Ca flow, using a smoothed local pore‐body capillary pressure curve appears to produce stable simulations. (2) Due to an approximation for the capillary driving force, IMP‐SIMS can deviate from quasi‐static solutions at equilibrium states especially in heterogeneous networks. (3) Both IMPES and IMP‐SIMS introduce mass conservation errors. The errors are small for networks with cubic pore bodies (less than 1.4% for IMPES and 1.2% for IMP‐SIMS). They become much greater for networks with square‐tube pore bodies (up to 45% for IMPES and 46% for IMP‐SIMS). Conversely, the new FI algorithm is numerically stable and mass conservative regardless of the flow regimes and pore geometries. It also precisely recovers the quasi‐static solutions at equilibrium states. The FI framework has been extended to include compressible two‐phase flow, multicomponent transport, and phase change dynamics. Example simulations of two‐phase displacements accounting for phase change show that evaporation and condensation can suppress fingering patterns generated during invasion.

Development of coupled two phase flow and geomechanics model to predict stress evolution in unconventional reservoirs with complex fracture geometry

Understanding stress changes associated with reservoir depletion is very important for developments of unconventional reservoirs such as completion of infill well. Fluid flow in reservoir often contains multiple phases mostly a mix of oil and water phase. Multiphase flow in porous media also involves in many applications including water injection. In this paper, we have upgraded our coupled fluid flow and geomechanics with Embedding Discrete Fracture Model (EDFM) to two phase flow using OpenFOAM, an open-source platform for CFD, in order to expand the model to wider range of applications including characterization of stress evolution associated with depletion of two-phase flow in unconventional reservoirs with complex fracture geometries. Coupled geomechanics and fluid-flow was developed based on well-known fixed-stress split, which is unconditionally stable and computationally efficient to simulate how stress changes during reservoir depletion. EDFM was coupled to the model to gain capability of simulating complex fracture geometries using structured grids. The two phase model is based on Implicit Pressure Explicit Saturation (IMPES) method. Validations of our single phase model and EDFM were conducted in our previous paper. Validation of two phase model was done against analytical solution of Buckley-Leverett and Capillary-gravity equilibrium problems. In addition, validation of EDFM part was conducted against local grid refinement. The model was used to run performance test on water injection-oil production problem with geomechanics and capillary effect. Non-planar fractures under production was also used to test for EDFM part. Finally, we applied our model in simulation of 3D multiple fractures under production in a fully 3D reservoir. The results indicate that created fracture geometries have a significant effect on stress distribution and reorientation induced by depletion. Characterization of stress evolution due to reservoir depletion in two phase flow will provide critical guidelines for optimization of completion designs especially in infill well drilling and further development in unconventional reservoirs.

Implicit linearization scheme for nonstandard two‐phase flow in porous media

SummaryIn this article, we consider a nonlocal (in time) two‐phase flow model. The nonlocality is introduced through the wettability alteration induced dynamic capillary pressure function. We present a monotone fixed‐point iterative linearization scheme for the resulting nonstandard model. The scheme treats the dynamic capillary pressure functions semiimplicitly and introduces an L‐scheme type stabilization term in the pressure as well as the transport equations. We prove the convergence of the proposed scheme theoretically under physically acceptable assumptions, and verify the theoretical analysis with numerical simulations. The scheme is implemented and tested for a variety of reservoir heterogeneities in addition to the dynamic change of the capillary pressure function. The proposed scheme satisfies the predefined stopping criterion within a few number of iterations. We also compared the performance of the proposed scheme against the iterative implicit pressure explicit saturation scheme.

A 6M digital twin for modeling and simulation in subsurface reservoirs

A 6M digital twin for modeling and simulation in subsurface reservoirs

Application of Implicit Pressure-Explicit Saturation Method to Predict Filtrated Mud Saturation Impact on the Hydrocarbon Reservoirs Formation Damage

Hydrocarbon reservoirs’ formation damage is one of the essential issues in petroleum industries that is caused by drilling and production operations and completion procedures. Ineffective implementation of drilling fluid during the drilling operations led to large volumes of filtrated mud penetrating into the reservoir formation. Therefore, pore throats and spaces would be filled, and hydrocarbon mobilization reduced due to the porosity and permeability reduction. In this paper, a developed model was proposed to predict the filtrated mud saturation impact on the formation damage. First, the physics of the fluids were examined, and the governing equations were defined by the combination of general mass transfer equations. The drilling mud penetration in the core on the one direction and the removal of oil from the core, in the other direction, requires the simultaneous dissolution of water and oil flow. As both fluids enter and exit from the same core, it is necessary to derive the equations of drilling mud and oil flow in a one-dimensional process. Finally, due to the complexity of mass balance and fluid flow equations in porous media, the implicit pressure-explicit saturation method was used to solve the equations simultaneously. Four crucial parameters of oil viscosity, water saturation, permeability, and porosity were sensitivity-analyzed in this model to predict the filtrated mud saturation. According to the results of the sensitivity analysis for the crucial parameters, at a lower porosity (porosity = 0.2), permeability (permeability = 2 mD), and water saturation (saturation = 0.1), the filtrated mud saturation had decreased. This resulted in the lower capillary forces, which were induced to penetrate the drilling fluid to the formation. Therefore, formation damage reduced at lower porosity, permeability and water saturation. Furthermore, at higher oil viscosities, due to the increased mobilization of oil through the porous media, filtrated mud saturation penetration through the core length would be increased slightly. Consequently, at the oil viscosity of 3 cP, the decrease rate of filtrated mud saturation is slower than other oil viscosities which indicated increased invasion of filtrated mud into the formation.

A Fundamental Model for Wormhole Formation Including Multiphase Flow

Summary The industry has been developing numerical models to simulate the wormholing phenomenon in carbonate matrix acidizing, both to save cost and time with experiments and to scale up the laboratory results to field scale. The two-scale continuum model is a fundamental model that has been successfully used for this end. Previous studies with this model only simulated single-phase flow: injection of acid into a water-saturated rock. However, significantly different behavior is observed experimentally by injecting acid into oil- or gas-saturated cores. In this work, we present a fundamental multiphase model for wormhole formation using the two-scale continuum approach, allowing the simulation of wormholing for acid injection into oil- or gas-bearing rocks, with different saturations. The two-scale continuum model represents the fluid flow and acid transport in the porous medium at the Darcy scale but calculates the acid-rock reaction with dissolution of the rock at the pore-scale. This model was implemented in transient 3D anisotropic form. Each phase occupies a volume fraction in each gridblock, defined by the porosity and fluid saturations. The fluid flow is calculated by solving the Darcy-Brinkman-Stokes equation, in which the relative permeabilities are functions of the saturations. The acid transport and reaction equations are solved, and as the acid injection proceeds and the rock is dissolved, the porosity increases in the gridblocks where dissolution occurs. The pore properties, such as permeability, pore radius, and specific surface area, are updated as porosity evolves, being scaled up from pore to Darcy scale. The simulation keeps track of the different fluid phases by calculating the saturations using the Implicit Pressure Explicit Saturation (IMPES) method. The developed model was implemented in an open-source computational fluid dynamics package and validated against experimental data. For the validation, the adjustable parameters in the model were calibrated so that the simulation results represent the different dissolution patterns and correctly reproduce the acid efficiency curves obtained experimentally. The same calibrated model was used to simulate the coreflooding experiments with water- and oil-saturated cores. The dissolution patterns (face dissolution, conical wormhole, dominant wormhole, etc.) and acid efficiency curves predicted by the new model match the experimental data. Other simulations presented include the shift in the acid efficiency curves observed for different oil viscosities, residual oil saturations, and different water saturations. To our knowledge, this is the first 3D two-scale continuum model to simulate wormhole propagation including multiphase flow. With adequate history match, it was shown to accurately predict the acidizing results for different fluid saturations, as observed in experiments.

A new physics-preserving IMPES scheme for incompressible and immiscible two-phase flow in heterogeneous porous media

In this work we consider a new efficient IMplicit Pressure Explicit Saturation (IMPES) scheme for the simulation of incompressible and immiscible two-phase flow in heterogeneous porous media with capillary pressure. Compared with the conventional IMPES schemes, the new IMPES scheme is inherently physics-preserving, namely, the new algorithm is locally mass conservative for both phases and it also enjoys another appealing feature that the total velocity is continuous in the normal direction. Moreover, the new scheme is unbiased with regard to the two phases and the saturations of both phases are bounds-preserving if the time step size is smaller than a certain value. The key ideas in the new scheme include that the Darcy flows for both phases are rewritten in the formulation based on the total velocity and an auxiliary velocity referring to as the capillary potential gradient, and the total discretized conservation equation is obtained by the summation of the discretized conservation equation for each phase. The upwind strategy is applied to update the saturations explicitly, and the upwind mixed finite element methods are used to solve the pressure–velocity systems which can be decoupled. We also present some interesting examples to demonstrate the efficiency and robustness of the new algorithm.

Non-classical flow modeling of spontaneous imbibition in spatially heterogeneous reservoirs

In this study, imbibition process is investigated in heterogeneous hydrocarbon reservoir using the non-classical (effective viscous finger) model equations proposed by Luo et al. (SPE Res. Eval. Engineering, 20(4), 1–16 2017). Counter-current spontaneous imbibition is analyzed in different heterogeneous porous media. The heterogeneity is mainly due to the spatial variation of the porosity and absolute permeability. The proposed mathematical model is solved considering the Implicit Pressure Explicit Saturation Scheme (IMPES) via finite difference method (FDM). Effect of the viscosity ratio, growth rate exponent of the fingering function, and the maximum fingering cross section on the normalized water saturation are analyzed. The present model is verified in the limiting case of classical model equations in homogeneous medium. It is noted that rate and extent of the counter-current imbibition are significantly controlled by the distribution of heterogeneous parameters like porosity and absolute permeability of the reservoir.