Method and System for Pore-Scale Modeling of a Multi-Phase Hydrocarbon Extraction Process

Abstract

A computer system models a hydrocarbon extraction process using a dynamic pore network model that is generated to represent a subterranean reservoir containing hydrocarbons as pores connected by throats. Solvent is injected into the subterranean reservoir to mobilize the hydrocarbons for extraction thereof. An iterative process may be repeated over time to determine hydrocarbon extraction based on changes in the molar balance of the components over time. A first set of characteristics for each pore is defined from which a second set of characteristics can be derived for two-phase pores. The molar balance of the components is determined based on the first and second sets of characteristics. The first set of characteristics is updated based on the molar balance of the components and the process can be repeated for subsequent times. Parameters for injection of the solvent may be adjusted and the iterative process repeated over time to identify preferred parameters.

Description

This application claims priority to U.S. Provisional Patent Application Serial No. 63/325,955 filed on Mar. 31, 2022, the contents of which are hereby incorporated by reference.

TECHNICAL FIELD

The present techniques relate to a method and system for pore-scale modeling to improve a multi-phased hydrocarbon extraction process.

BACKGROUND

This section is intended to introduce various aspects of the art, which may be associated with exemplary embodiments of the present techniques. This discussion is believed to assist in providing a framework to facilitate better understanding of particular aspects of the present techniques. Accordingly, it should be understood that this section should be read in this light, and not necessarily as admissions of prior art.

Hydrocarbons may be extracted from a subterranean reservoir using various in-situ methods that require drilling one or more wells into the subterranean reservoir. Heating agents or solvents may be injected into the subterranean reservoir through one of these wells in order to assist in the extraction of the hydrocarbons from the subterranean reservoir. Injection of such heating agents or solvents promotes pressure and gravity-assisted drainage of hydrocarbons enabling one of the wells to capture the draining hydrocarbons.

When such heating agents or solvents are employed, there is a resulting flow that includes more than one phase. When an aqueous heating agent is injected into the subterranean reservoir, there is a three-phase flow involving water, oil and gas during extraction. The use of a hydrocarbon-based solvent can be considered, for modeling purposes, to result in a two-phase flow involving oil and solvent vapour during extraction.

Computer systems for computer modeling to predict production rates for hydrocarbons is useful to improve production rates and operational parameters of the equipment used in an extraction process within the subterranean reservoir. Since the subterranean reservoir can be considered to be a porous media, computer modelling of the extraction process is often based on an application of Darcy’s law in a reservoir model with large numerical grids to represent the entire subterranean reservoir. However, the generalized Darcy’s law for two-phase flow in a porous media may not be applied when the relationship between the pressure gradient and flow rate is no longer linear (see, for example, Y. Gao, Q. Lin, B. Bijeljic, and M.J. Blunt, “Pore-scale dynamics and the multiphase Darcy law” Physical Review Fluids, 2020 5(1) and Y. Zhang, B. Bijeljic, Y. Gao, Q. Lin, and M.J. Blunt, “Quantification of Nonlinear Multiphase Flow in Porous Media” Geophysical Research Letters, 2021 48(5).)

Most existing computer models of hydrocarbon extraction are based on reservoir-scale models. However, such large-scale modelling of a solvent-based extraction process typically has difficulty accounting for mass transfer between a vapor chamber formed by an injected solvent and the hydrocarbons (such as heavy oil or bitumen) due to a mixing zone therebetween. In a real implementation, this mixing zone is measured in millimeters making it too thin to be effectively modeled with a large sized grid that is measured in meters.

The foregoing examples of the related art and limitations related thereto are intended to be illustrative and not exclusive. Other limitations of the related art will become apparent to those of skill in the art upon a reading of the specification and a study of the drawings.

SUMMARY

The following embodiments and aspects thereof are described and illustrated in conjunction with systems, tools and methods which are meant to be exemplary and illustrative, not limiting in scope. In various embodiments, one or more of the above-described problems have been reduced or eliminated, while other embodiments are directed to other improvements.

In an embodiment, a method of operating a computer system for determining hydrocarbon extraction from a subterranean reservoir, the method comprising: (a) generating a pore network model of the subterranean reservoir representing a structure of the subterranean reservoir and conditions in the subterranean reservoir, the structure of the subterranean reservoir being represented by the pore network model as a plurality of pores connected to each other by a plurality of throats, the plurality of pores and the plurality of throats containing components including hydrocarbons; (b)establishing operational parameters for injection of a solvent into the subterranean reservoir represented by the pore network model, wherein the solvent when injected becomes one of the components in at least one from the plurality of pores and the plurality of throats; (c) based on the operational parameters for injection of the solvent, performing: (i) defining, for each of the plurality of pores, a first set of characteristics based on a molar density of each component in a pore; (ii) deriving, for each of the plurality of pores based on phase conditions in the pore, a second set of characteristics for the pore from the first set of characteristics for the pore; (iii) determining, for each of the plurality of throats, phase conditions in a throat from the plurality of throats based on the phase conditions in ones of the plurality of pores that are connected to the throat; (iv) determining a molar balance of the components within the plurality of pores and the plurality of throats in the pore network model based on the first set of characteristics and the second set of characteristics for each of the plurality of pores and the phase conditions for each of the plurality of throats; (v) determining an updated first set of characteristics for each of the plurality of pores based on the molar balance of the components within the pore network model; (vi) repeating (ii) to (v) over a predetermined time at consecutive times to determine hydrocarbon extraction from the subterranean reservoir over time given the operational parameters for injection of the solvent based on changes in the molar balance of the components over time representing mass transfer of the components within the pore network model; and (d) adjusting the operational parameters for injection of the solvent and repeating (c) with the adjusted operational parameters to identify operational parameters for injection of the solvent that provide a predetermined hydrocarbon extraction.

In an embodiment, a computer system for determining hydrocarbon extraction from a subterranean reservoir, the computer system comprising: an interface for receiving information about properties of the subterranean reservoir; a processor configured to: (a) generate a pore network model of the subterranean reservoir based on the received information about the properties of the subterranean reservoir, the pore network model representing a structure of the subterranean reservoir and conditions in the subterranean reservoir, the structure of the subterranean reservoir being represented by the pore network model as a plurality of pores connected to each other by a plurality of throats, the plurality of pores and the plurality of throats containing components including hydrocarbons; (b) establish operational parameters for injection of a solvent into the subterranean reservoir represented by the pore network model, wherein the solvent when injected becomes one of the components in at least one from the plurality of pores and the plurality of throats; (c) based on the operational parameters for injection of the solvent, perform: (i) define, for each of the plurality of pores, a first set of characteristics based on a molar density of each component in a pore; (ii) derive, for each of the plurality of pores based on phase conditions in the pore, a second set of characteristics for the pore from the first set of characteristics for the pore; (iii) determine, for each of the plurality of throats, phase conditions in a throat from the plurality of throats based on the phase conditions in ones of the plurality of pores that are connected to the throat; (iv) determine a molar balance of the components within the plurality of pores and the plurality of throats in the pore network model based on the first set of characteristics and the second set of characteristics for each of the plurality of pores and the phase conditions for each of the plurality of throats; (v) determine an updated first set of characteristics for each of the plurality of pores based on the molar balance of the components within the pore network model; (vi) repeat (ii) to (v) over a predetermined time at consecutive times to determine hydrocarbon extraction from the subterranean reservoir over time given the operational parameters for injection of the solvent based on changes in the molar balance of the components over time representing mass transfer of the components within the pore network model; and (d) adjusting the operational parameters for injection of the solvent and repeating (c) with the adjusted operational parameters to identify operational parameters for injection of the solvent that provide a predetermined hydrocarbon extraction.

In an embodiment, a computer readable medium having stored thereon statements and instructions that when executed by a computer system operate the computer system to determine hydrocarbon extraction from a subterranean reservoir, the statements and instructions when executed by the computer system cause the computer system to: (a) generate a pore network model of the subterranean reservoir representing a structure of the subterranean reservoir and conditions in the subterranean reservoir, the structure of the subterranean reservoir being represented by the pore network model as a plurality of pores connected to each other by a plurality of throats, the plurality of pores and the plurality of throats containing components including hydrocarbons; (b) establish operational parameters for injection of a solvent into the subterranean reservoir represented by the pore network model, wherein the solvent when injected becomes one of the components in at least one from the plurality of pores and the plurality of throats; (c) based on the operational parameters for injection of the solvent, perform: (i) define, for each of the plurality of pores, a first set of characteristics based on a molar density of each component in a pore; (ii) derive, for each of the plurality of pores based on phase conditions in the pore, a second set of characteristics for the pore from the first set of characteristics for the pore; (iii) determine, for each of the plurality of throats, phase conditions in a throat from the plurality of throats based on the phase conditions in ones of the plurality of pores that are connected to the throat; (iv) determine a molar balance of the components within the plurality of pores and the plurality of throats in the pore network model based on the first set of characteristics and the second set of characteristics for each of the plurality of pores and the phase conditions for each of the plurality of throats; (v) determine an updated first set of characteristics for each of the plurality of pores based on the molar balance of the components within the pore network model; (vi) repeat (ii) to (v) over a predetermined time at consecutive times to determine hydrocarbon extraction from the subterranean reservoir over time given the operational parameters for injection of the solvent based on changes in the molar balance of the components over time representing mass transfer of the components within the pore network model; and (d) adjust the operational parameters for injection of the solvent and repeat (c) with the adjusted operational parameters to identify operational parameters for injection of the solvent that provide a predetermined hydrocarbon extraction.

In addition to the exemplary aspects and embodiments described above, further aspects and embodiments will become apparent by reference to the drawings and by study of the following detailed descriptions.

BRIEF DESCRIPTION OF THE DRAWINGS

The advantages of the present techniques will be better understood by referring to the following detailed description and the attached drawings. It is intended that the embodiments and figures disclosed herein are to be considered illustrative rather than restrictive. Embodiments will now be described, by way of example only, with reference to the attached figures, in which:

FIG. 1 illustrates a computer system for computer modeling of a hydrocarbon extraction process using a pore network model;

FIG. 2 illustrates a first embodiment of a flow for a computer system of computer modelling of a hydrocarbon extraction process using a pore network model;

FIG. 3 illustrates an example of a fluid configuration and cross section of a throat used in a pore network model;

FIG. 4 illustrates a second embodiment of a flow for a computer system of computer modelling of a hydrocarbon extraction process using a pore network model;

FIG. 5 illustrates a third embodiment of a flow for a computer system of computer modelling of a hydrocarbon extraction process using a pore network model;

FIG. 6 is a simplified exemplary embodiment of computer modelling of a hydrocarbon extraction process using a pore network model; and

FIGS. 7A to 7D illustrate exemplary embodiments of computer modelling of a hydrocarbon extraction process using a pore network model.

DESCRIPTION OF SELECTED EMBODIMENTS

In the following detailed description section, specific embodiments of the present techniques are described. However, to the extent that the following description is specific to a particular embodiment or a particular use of the present techniques, this is intended to be for exemplary purposes only and simply provides a description of the exemplary embodiments. Accordingly, the techniques are not limited to specific embodiments described below but rather, include all alternatives, modifications, and equivalents.

At the outset, for ease of reference, certain terms used in this application and their meaning as used in this context are set forth below. To the extent a term used herein is not defined below, it should be given the broadest definition persons in the pertinent art have given that term as reflected in at least one printed publication or issued patent. Further, the present techniques are not limited by the usage of the terms shown below, as all equivalents, synonyms, new developments and terms or techniques that serve the same or a similar purpose are considered to be within the scope hereof.

A “hydrocarbon” is an organic compound present in a subterranean reservoir that is being extracted from the subterranean reservoir as a result of a production or extraction process.

“Bitumen” is the hydrocarbon component found in oil sands.

“Heavy oil” includes bitumen as the hydrocarbon component, as well as lighter materials that may be found in a sand or carbonate reservoir. Heavy oil extraction generally relies on decreasing viscosity using an increased temperature or a solvent. The heavy oil will drain using gravity and/or through the use of a mobilizing fluid (e.g., solvent vapour), for example by way of a pressure gradient, to enable extraction of the heavy oil.

“Solvent” is an agent, heated or unheated, that may dilute or dissolve hydrocarbons in the subterranean reservoir, such as heavy oil, and may reduce a viscosity of the hydrocarbons or components containing the hydrocarbons. Injection of a solvent into a subterranean reservoir may assist in hydrocarbon extraction by reducing viscosity and/or through employment of a pressure gradient as the solvent is injected at a pressure. The solvents used in extraction of hydrocarbons in the subterranean reservoir may include hydrocarbon components but are referred to herein simply as solvents regardless of the presence of hydrocarbon components.

A “subterranean reservoir” is a subsurface rock or sand reservoir from which a resource such as hydrocarbon can be recovered.

“Reservoir matrix” refers to the solid porous material forming the structure of the subterranean reservoir and surrounding pore spaces within the subterranean reservoir. The pore spaces may contain hydrocarbon resources.

“Pore” refers to the pore spaces in the subterranean reservoir that are surrounded by the reservoir matrix. The pores may contain the hydrocarbon resources that are to be extracted, such as heavy oil.

“Throat” is a tube of any geometrical cross section that connects two pores.

“Model” is a computer model of a hydrocarbon extraction process relying on an interaction of hardware components, software components and specially developed components that cooperate together to form a computer system with an electronic model of a subterranean reservoir and actions and processes within the subterranean reservoir over time.

The terms “approximately,” “about,” “substantially,” and similar terms are intended to have a broad meaning. It should be understood that these terms provide for a description of certain features without restricting the scope of these features to the precise numeral ranges provided.

The terms “the”, “a” and “an” are not necessarily limited to mean only one, but rather are inclusive and open ended so as to include, optionally, multiple such elements.

“At least one,” is understood to mean at least one entity selected from any one or more of the entity in the list of entities, but not necessarily including at least one of each and every entity specifically listed within the list of entities and not excluding any combinations of entities in the list of entities. The expressions “at least one,” “one or more,” and “and/or” are open-ended expressions that are both conjunctive and disjunctive.

Background

A subterranean reservoir contains a reservoir matrix of a solid porous material that forms the structure for the subterranean reservoir. Within the reservoir matrix are pores that are largely surrounded by the reservoir matrix but that may be connected to other pores via throats. Within the subterranean reservoir, hydrocarbons reside in the pores and throats.

Computer modeling of hydrocarbon extraction processes has been used to investigate and predict hydrocarbon extraction from subterranean reservoirs. Such computer modeling has been used to determine placement of wells within a subterranean reservoir as well as optimal operating parameters for the wells and any substances injected into the subterranean reservoir (e.g., vapor, solvents, steam, etc.) during the extraction process. Information on the subterranean reservoir must be obtained to form a basis for computer modeling of the hydrocarbon extraction process. While computer modeling of a hydrocarbon extraction process can be used during all stages of the extraction process (e.g., planning, early and late) to adjust to changing conditions in the subterranean reservoir during hydrocarbon extraction, such computer modelling can be computationally expensive.

In a solvent-based hydrocarbon extraction process, the solvent injected into the subterranean reservoir through one well forms a vapor chamber. The boundary between the vapor chamber and the hydrocarbons in the subterranean reservoir can involve complex mechanisms at their interface, such as: diffusion, two-phase flow, dispersion and hydrocarbon precipitation. There is a mass transfer layer between the injected solvent and the hydrocarbons in the subterranean reservoir that consists of a dynamic gas-liquid capillary mixing zone and a single phase liquid zone where dispersive forces mix the injected solvent and the hydrocarbons in the subterranean reservoir. Within the two-phase mixing zone at the interface, mass transfer rate is enhanced by interface renewal. As the diluted hydrocarbons mobilize within the subterranean reservoir and drain along the mass transfer layer, a new interface between the injected solvent and the hydrocarbons in the subterranean reservoir is created resulting in surface renewal of the interface.

Computer modeling of a hydrocarbon extraction process using numerical methods have attempted to investigate mass transfer and convective flow mechanisms during solvent-based extraction processes. However, such computer modelling using numerical methods is macroscopic and has difficulty describing actions in the subterranean reservoir on a pore scale. In such numerical methods of modelling a hydrocarbon extraction process, determining a dispersion coefficient at the interface between the solvent and the hydrocarbons in the subterranean reservoir is complicated. This is because mass transfer depends on the microstructure of the subterranean reservoir, gravity, capillary effects, interface renewal and multiphase flow. In Q. Wang, X. Jia, Z. Chen, “Mathematical Modeling of the Solvent Chamber Evolution in a Vapor Extraction Heavy Oil Recovery Process”, Fuel, 186, (2016), 339-349 and Q. Wang, X. Jia, and Z. Chen, “Modelling of dynamic mass transfer in a vapor extraction heavy oil recovery process”, The Canadian Journal of Chemical Engineering, 2017. 95(6): p.1171-1180, a semi-analytical model was developed to model a solvent chamber boundary to predict the evolution of the vapor chamber during a vapor extraction process (VAPEX) and the production flow rate; however, the effective diffusion coefficient incorporated all mechanisms contributing to oil-solvent mixing. This coefficient was obtained by fitting production data which reduces the predictive capability of the model.

Computer modeling of hydrocarbon extraction processes using pore-scale modelling has been considered as a means to overcome the problems with large-scale numerical modelling. A dynamic pore network model may be used as an improved method to describe the complex mechanisms of solvent-based extraction processes of hydrocarbons from subterranean reservoirs.

S. Chen, C. Qin, B. Guo, “Fully Implicit Dynamic Pore-Network Modeling of Two-Phase Flow and Phase Change in Porous Media”, Water Resources Research, 10.1029/2020WR028510, (2020) describes a dynamic pore network model for two-phase flow; however, this process relies on the use of several (2(N_c-1)+2) primary parameters to represent two-phases for each pore. This increases the computational complexity due to a size of a matrix used in the calculations and computational cost of this model, making implementation reliant on larger processing power and storage devices. Further, there is no consideration in this model of diffusion in a liquid phase and as such there is no modelling of mass transfer by diffusion in the liquid phase.

Overview

Computer modelling of a hydrocarbon extraction process using a dynamic pore network model may be used to address concerns of computational complexity and representation of real-world phenomena during the hydrocarbon extraction process. The dynamic pore network model uses information about a subterranean reservoir to simulate the void space (pores) in the reservoir matrix (porous media) as a network of interconnected geometries of the pore space over time to enable an examination of the impact of time on phase distribution (e.g., liquid phase and gas phase) within the subterranean reservoir. Within a structure of the dynamic pore network model representing the subterranean reservoir, the geometries of the throats connecting the pores may be represented by simple geometries such that their hydraulic conductivities can be analytically or numerically estimated. The information about the subterranean reservoir may be obtained, for example, through capture devices (as well as other techniques) that may scan a subterranean reservoir to obtain information about the reservoir matrix as well as other information about the conditions within the subterranean reservoir.

The dynamic pore network model may describe diffusion of solvent into hydrocarbons in the subterranean reservoir and its impact on hydrocarbon viscosity with changing characteristics of the solvent such as solvent concentration, flow rate, pressure of the injected solvent, etc. When the hydrocarbons in the subterranean reservoir are sufficiently dilute as a result of the injected solvent, they become mobile and move in a dynamic two-phase flow through the reservoir matrix. The two-phase flow pattern in the reservoir matrix has an impact on mass transfer and hydrocarbon extraction.

Time dependent multiphase flow is simulated by the computer modelling using the dynamic pore network model with an evolution of the interface between the solvent and the hydrocarbons in the subterranean reservoir. Capillary and viscous forces at the interface between hydrocarbon and solvent as well as drainage, imbibition, movement of trapped phases and intermittent flow are taken into account during computer modelling of the hydrocarbon extraction process using the dynamic pore network model. Conservation of fluid components and thermodynamic equilibrium are considered for each pore. For a pore where two-phases coexist, mass conservation and thermodynamic equilibrium are used to determine characteristics (mole fraction of components in each phase, liquid or gas phase pressure, phase saturation) representing each pore at a given temperature. In order to reduce the computational cost of the pore network model, a first set of characteristics are used to define a state of a pore. A second set of characteristics may be determined for each pore using the first set of characteristics for that pore to further describe conditions in the pore.

By capturing the impact of transient pore-scale mechanisms on mass transfer in liquid and gas phases, the impact of operating parameters of the solvent injection during the hydrocarbon extraction process can be determined along with an identification of preferrable operating parameters and macroscopic properties that can be used in a reservoir scale model.

Direct numerical methods of modelling multiphase flow in the pore space of the subterranean reservoir are computationally too expensive to simulate a representative elementary volume large enough to represent macroscopic properties and study transient multiphase flow phenomena and transport mechanisms. Computer modelling using the dynamic pore network model can simulate large representative elementary volumes and study transient multiphase and transport phenomena in macroscopic pore scale porous media.

Detailed Description

FIG. 1 illustrates a system 10 for computer modeling of a hydrocarbon extraction process using a pore network model.

The system 10 includes a capture device 50 in communication with a computer system 70 that includes a processor 20, a memory 90. and specialized modules including a pore network generation module 22, a capture device interface module 24, an operational control module 26, a first characteristics determination module 28, a second characteristics determination module 30, a throat conditions module 32, a hydrocarbon extraction determination module 34, a time control module 36, a molar balance determination module 38, an iterative control module 40, and a pore conditions module 42.

The capture device 50 obtains information about the structure of a subterranean reservoir that has void spaces (e.g., pores) in a porous media connected together by throats. The structure of the subterranean reservoir may be extracted from scans of the subterranean reservoir that are captured by the capture device 50. The capture device 50 may be, for example, computerized tomography (CT), or any other mechanism for obtaining in-situ images and measurements of the subterranean reservoir.

In addition to information about the structure of the subterranean reservoir, information about conditions in the subterranean reservoir may be obtained by the capture device 50, a component of the capture device 50 or the capture device 50 may comprise multiple devices, one of which may obtain this information. For example, information about conditions in the subterranean reservoir may be obtained by drilling a core sample or some other invasive or non-invasive technique. The conditions in the subterranean reservoir represent the state of the pores and throats in the subterranean reservoir at a given point in time. These conditions may include, for example, saturation of liquid phases in the pores of the subterranean reservoir, liquid pressure in the subterranean reservoir, or molar density of components in regions of the subterranean reservoir, and, if the hydrocarbons extraction process has already started, a molar density of the solvent that has been injected into the subterranean reservoir. The capture device 50 may include mechanisms by which the conditions can be directly measured or may include mechanisms by which other information can be measured and the above conditions can be derived.

The capture device 50 is in communication with the computer system 70 through the capture device interface module 24 via a network (not shown) or by way of a direct connection that may be wired or wireless. The capture device interface module 24 obtains the information about the subterranean reservoir that has been captured by the capture device 50 and either stores this information in the memory 90 or provides this information directly to the pore generation network module 22 for the generation of a pore network model.

In addition to information from the capture device 50, the memory 90 may also contain previously obtained information about the subterranean reservoir.

The processor 20 may interact with the specialized modules including the pore network generation module 22, the capture device interface module 24, the operational control module 26, the first characteristics determination module 28, the second characteristics determination module 30, the throat conditions module 32, the hydrocarbon extraction determination module 34, the time control module 36, the molar balance determination module 38, the iterative control module 40, and the pore conditions module 42 to provide for computer modeling of a multiphase hydrocarbon extraction process using a dynamic pore network model as set out in the flows illustrated in FIGS. 2, 4 and 5.

The pore network generation module 22 obtains information from the capture device 50 either through the capture device interface module 24 or as previously stored in the memory 90. The obtained information is used to generate a structure of the pore network model by representing pores in the subterranean reservoir that are connected via throats as discussed in further details with respect to FIGS. 2, 4 and 5. Conditions in the subterranean reservoir are also derived from this information and are used to set the initial conditions for the pore network model.

The operational control module 26 sets and adjusts the operating parameters for injection of a solvent into the subterranean reservoir in the pore network model. These operating parameters may include at least some of the following: placement of injection wells and/or production wells within the structure of the pore network model, a flow rate for the solvent injected into the subterranean reservoir, a gas pressure for the injected solvent, solvent composition, and solvent additives. Results from the hydrocarbon extraction determination module 34 and the molar balance determination module 38 may be provided back to the operational control module 26 to adjust the operating parameters within the subterranean reservoir.

The first characteristics determination module 28 determines and updates a first set of characteristics that define a state of a pore. The characteristics that are included in the first set of characteristics may depend on whether the pore is in one phase (i.e., either liquid phase or gas phase) or two phases (i.e., both liquid phase and gas phase). As the extraction process is modeled over time, the values for the first set of characteristics for each pore will change as the hydrocarbons are extracted and may be updated with results from the molar balance determination module 38. Additionally, whether the pore is in one phase or two phases will also change throughout the extraction process over time and, as such, the characteristics that are included in the first set of characteristics may change as they are updated based on phase conditions of the pore as determined by the pore conditions module 42. The characteristics that make up the first set of characteristics may vary depending on constant parameters that are established to represent conditions used in thermodynamic calculations of the second set of characteristics as is discussed in respect of FIGS. 2, 4 and 5.

The pore conditions module 42 determines the conditions in each pore of the pore network module. The conditions determined for each pore may include the phase conditions of each pore, which describes whether a pore is in a single phase (liquid or gas) or two phases (gas and liquid). The conditions for each pore may be provided to both the first characteristics determination module 28 and the second characteristics determination module 30. Determination of the pore conditions is discussed in respect of FIGS. 2, 4 and 5.

The second characteristics determination module 30 determines a second set of characteristics for each pore that contains two phases. The second characteristics determination module 30 is in communication with the pore conditions module 42 to obtain the phase conditions for each pore to determine whether or not a pore contains two phases. The second characteristics determination module 30 determines the second set of characteristics for such pores. The second set of characteristics comprises characteristics or properties that can be derived from the first set of characteristics. Determination of the second set of characteristics is discussed in respect of FIGS. 4 and 5.

The throat conditions module 32 determines the conditions in each throat of the pore network module. The conditions determined for each throat are based on the conditions in pores connected to the throat and may include phase conditions or throat occupancy of whether a throat is occupied by a single phase (liquid or gas) or two phases (gas and liquid), and phase conductivities in the throat. Determination of the throat conditions is discussed in respect of FIGS. 2, 4 and 5.

The molar balance determination module 38 determines the molar balance for each component, including the solvent and hydrocarbon components, in the pore network model taking into consideration conservation of mass. This may then be used to examine mass transfer of the components within the pore network model. The molar balance for each component may be determined by a numerical analysis method, such as an iterative process using the Newton method with the result of multiple iterations being a completed molar balance for each component for a given time. The molar balance is determined by the molar balance determination module 38 for each time in the computer modelling of the hydrocarbon extraction process. Further discussion of the molar balance determination is found below in respect of FIGS. 2, 4 and 5.

The iterative control module 40 controls the iterative process of the molar balance determination module 38. As an iterative numerical analysis method may be used to determine the molar balance, the iterative control module 40 controls whether an iteration is complete and the next iteration can occur or whether all iterations for a given time t are complete and the molar balance for the given time t has been found. The iterative control module 40 may assess a change in a predetermined characteristic (e.g., the first set of characteristics) between iterations to determine if that change is greater than an acceptable limit which may indicate that the molar balance has been not found and iterations at the molar balance determination module 38 are to continue.

The time control module 36 controls a progression of time within the computer modelling of the hydrocarbon extraction process. The time control module 36 may interact with the molar balance determination module 38 and the iterative control module 40 to check on the iterations at the molar balance determination module 38 and check on whether or not a difference in a predetermined characteristic between times is sufficiently small or large. Within the iterations at the molar balance determination module 38, if the change in a predetermined characteristic is too great between consecutive iterations, then the time step Δt between times may be changed to ensure that the iterations at the molar balance determination module 38 eventually find the molar balance. Likewise, the time control module 36 may check on a change in a predetermined characteristic between consecutive times to ensure that the change is not too great, otherwise the time step Δt may be altered. The time control module 36 controls the progression of time from the beginning of the hydrocarbon extraction process in the pore network model to an end of the time for the hydrocarbon extraction process. The end time may be controlled by either a set time or may be based on results from the hydrocarbon extraction module 34 being obtained.

The hydrocarbon extraction module 34 determines the extraction of hydrocarbons during the various times by assessing the change in the molar balance of the hydrocarbons in the pore network model. Determination of the extraction of hydrocarbons is discussed further in respect of FIGS. 2, 4 and 5.

The interaction of the specialized modules in the computer system 70, the flow between the specialized modules and which modules are used will vary in the flows illustrated in respect of FIGS. 2, 4 and 5.

FIG. 2 illustrates an exemplary flow 100 for a computer system of computer modelling of a multiphase extraction process using a dynamic pore network model.

Information describing the subterranean reservoir is obtained in step 102. The information describing the subterranean reservoir may include a structure of the subterranean reservoir and conditions therein.

The structure of the subterranean reservoir may be represented by void spaces (e.g., pores) in the porous media of the subterranean reservoir that are connected together by throats. The structure of the subterranean reservoir may be extracted from scans of the subterranean reservoir that are captured by a capture device. The measurements obtained by the capture device may be used by the computer system in step 104 for the generation of a pore network model to be used in the flow 100.

The structure of the pore network model may also be randomly generated or may be generated based on known qualities of the porous media (e.g., subterranean reservoir) in consideration.

The conditions in the subterranean reservoir represent the state of the pores and throats in the subterranean reservoir at a given point in time. These conditions may include, for example, a molar density of the hydrocarbons and gas or liquid pressure in each pore and may also include liquid saturation in each pore. Additionally, phase conditions in each pore may be set or determined by the conditions in the subterranean reservoir to establish whether a pore is in one state (i.e., liquid or gas) or two states (i.e., liquid and gas) at the given point in time.

In a case where the hydrocarbon extraction process has not started in the subterranean reservoir, the conditions that are obtained may be initial conditions of the subterranean reservoir. However, conditions at other points in time such as after the hydrocarbon extraction process has started or towards the end of the hydrocarbon extraction process may also be used. The conditions in the subterranean reservoir, regardless of the given point in time when they are obtained, may be obtained from information measured by the capture device, associated devices or based on information obtained through other means.

A pore network model of the subterranean reservoir is generated by a computer system in step 104 based on the information about the subterranean reservoir that was obtained in step 102. The pore network model includes a structure representing the structure of the subterranean reservoir and the conditions in the subterranean reservoir.

The pore network model simulates the void space in the subterranean reservoir as a network of pores connected to each other by throats. The information about the structure of the subterranean reservoir may classified by the computer system to determine each area as pores representing larger void spaces, reservoir matrix representing solid material such as rock, and throats representing void spaces connecting the pores. After this classification, the pores and throats may be assembled by the computer system to form the structure of the pore network model.

The geometries of the throats and the pores may be determined based on measurements from the capture device or may be simplified by idealized geometries based on assumptions about the subterranean reservoir. For example, the pores may be modeled as cubes to reduce computational complexity of the flow 100. However, other geometries or a combination of different geometries may also be used to represent the pores. The throats may be represented by simple geometries to simplify computational complexity of the flow 100. Geometries and cross sections of the pores and throats are discussed in more detail in connection with FIG. 3.

In addition to the structure of pore network model being generated, the conditions in the subterranean reservoir are incorporated with the structure to form the pore network model used in the flow 100. These conditions may be considered to be the initial conditions for the pore network model and may include pore saturation and solvent mole fraction based on the conditions in the subterranean reservoir that are obtained in step 102. For example, if wells are either present or are being drilled in the subterranean reservoir to be modeled, pore saturation may be selected based on existing conditions that have been determined or estimated through known mechanisms. Likewise, the solvent mole fraction used for the conditions in the subterranean reservoir may be based on known solvent mole fractions that are being used for similar purposes. The initial conditions for the pore network model may also include the phase conditions for the pore that establish whether the pore is in one phase (i.e., liquid or gas) or two phases (i.e., liquid and gas).

The pore network model that is generated in step 104 simulates multiphase flow, and in particular two phase flow through the porous media that forms the structure of the pore network model. In the pore network model as used in FIG. 2, it is assumed that the gas phase in the two-phase flow is the non-wetting phase and the liquid phase is the wetting phase. However, it is possible for the gas phase to be the wetting phase and the liquid phase to be the non-wetting phase.

In step 106, a first set of characteristics for each pore in the pore network model is determined. The first set of characteristics for each pore in a single (or liquid) phase defines a primary set of characteristics for the pore and may be given by

$X={\left[x_i^{k=solv}{P}_i^l\right]}_{i=1,N_{pore}\cdot }$

The first set of characteristics for a pore in two phases (i.e., both gas phase and liquid phase) may be given by

$X={\left[S_i^l{P}_i^l\right]}_{i=1,N_{pore}}$

where

$x_i^k$

is the mole fraction of component k in liquid phase in pore i,

$P_i^l$

is the liquid pressure in pore i,

$S_i^l$

is the liquid phase saturation in pore i, N_pore is the number of pores in the pore network model. In FIGS. 4 and 5 different characteristics may form the first set of characteristics and a second set of characteristics. As the extraction process is modeled over time, values for the first set of characteristics for each pore will change as the hydrocarbons are extracted from the subterranean reservoir that is modeled by the pore network model. Further, during the extraction process, the pores will change from containing one phase to two phases and vice versa over time, which will alter the characteristics in the first set of characteristics that is used to describe the pore.

In using the above to describe a pore, assumptions may be made with respect to behavior of the components in the pores. For example, in the case of a subterranean reservoir containing hydrocarbons that are being removed with the assistance of a solvent, it is assumed that when the hydrocarbon is bitumen, the hydrocarbon is a single component that is present only in the liquid phase while the solvent is present in both the liquid and gas phases. Given that bitumen may be composed of multiple components, bitumen may be modeled as multiple pseudo-components to represent all of the components contained in bitumen in some embodiments. It may be assumed that the solubility of solvent in hydrocarbon is constant and independent of variations in pressure. It may also be assumed that the gas phase contains exclusively solvent.

Phase conditions of each throat is determined based on connected pores using pore capillary pressures, throat invasion and snap-off events. The phase conditions may include which phase occupancy indicating which phases (e.g., liquid phase and/or gas phase) are present in the throat and phase conductivities for the throat.

FIG. 3 illustrates example geometries of the throats in the pore network model. A throat 202 may be modeled as a cube with a wetting phase 204 residing in the corners and along the edges. The wetting phase 204 (e.g., liquid) wets the surface of the throat 202 whereas a nonwetting phase 206 (e.g., gas) does not and tends to reside near the center of the throat 202. While FIG. 3 illustrates the throat 202 (and by extension the throat) as having a square cross section to reduce computational complexity, other geometries may be used. In general, a shape of an interface between teh wetting phase and the nonwetting phase shape depends on a contact angle between them. Other geometries for the pores and throats may be determined by measured properties in the subterranean reservoir. While FIG. 3 illustrates an example geometry of a throat, a similar geometry between the interface of the wetting phase and the nonwetting phase may also apply to pores in the pore network model.

Throat occupancy reflects whether a throat is occupied by a single phase (liquid or gas) or two phases (gas and liquid). Like a pore, the throat may have a square cross section in which, if a throat is occupied by two phases, the wetting phase will be at the corners and along the walls while the nonwetting phase will be near the center of the throat. Initially, the throat occupancy for the throat depends on a liquid saturation of the pores connected to the throat.

A pressure difference, the pore capillary pressure, across a curved interface surface between the liquid phase and the gas phase in the subterranean reservoir can result in capillary action. The pore capillary pressure may depend on a radius 210 R_i of the pore and liquid saturation

$S_i^l.$

When the pore capillary pressure in the pores connected to a throat exceeds a threshold entry capillary pressure of the connected throat (i.e., threshold value for throat invasion), the throat becomes invaded with the gas phase. This invasion of gas into the throat reduces conductivity of the liquid in the throat. The gas phase invades the throat when the pore capillary pressure in the neighboring pores body exceeds an entry capillary pressure of the throat.

When the capillary pressure in a throat drops below a threshold value, liquid in the throat swells and snap-off event occurs to fill the throat with the liquid phase. The throat capillary pressure is the maximum capillary pressure in the connected pores. The capillary pressure in the throat is assumed to be the maximum capillary pressure in the connected pores.

Based on the above, the throat occupancy can be determined and may be used to represent the state of the throat as being in either a single phase or two phases.

In step 110, phase conductivities of the liquid phase and gas phase in the throats are determined.

When the throat is in a single phase (i.e., liquid phase, also valid for gas phase), the conductivity of the phase α (gas phase or liquid phase) may be determined by (Chen, C. Qin, B. Guo, Fully Implicit Dynamic Pore-Network Modeling of Two-Phase Flow and Phase Change in Porous Media. 2020):

$K_{ij}^{\alpha }=\frac{0.5623G_{ij}{\left(A_{ij}\right)}^2}{{\mu }^{\alpha }{L}_{ij}}$

where

$G_{ij}=\frac{A_{ij}}{{\left(C_{ij}\right)}^2}$

is a shape factor for the throat cross section, A_ij, C_ij, and L_ijare cross sectional area, perimeter and length of the throat connecting pores i and j, respectively, and µ is the viscosity of the phase α.

When a throat is invaded and two phases coexist, the gas phase may be considered to occupy a central portion of the throat and have a generally round geometry while the liquid phase may be considered to occupy corners of the throat surrounding the gas phase.

The conductivity of the gas phase may be determined by (Chen, C. Qin, B. Guo, Fully Implicit Dynamic Pore-Network Modeling of Two-Phase Flow and Phase Change in Porous Media. 2020):

$K_{ij}^{nw}=\frac{0.5623G_{ij}^{nw}{\left(A_{ij}^{nw}\right)}^2}{{\mu }^{nw}{L}_{ij}}$

is a shape factor of the gas phase region in the throat cross section,

$A_{ij}^{nw}$

is a cross sectional area of the gas phase region given by

$A_{ij}^{nw}=4R_{ij}^2-4r_{ij}^2\left[\frac{\sin \left(\frac{\pi }{4}-\theta \right)}{\sin \frac{\pi }{4}}\cos \theta -\left(\frac{\pi }{4}-\theta \right)\right]$

where R_ij is a radius of an inscribed circle of the throat cross section and r_ij is a radius of a fluid-fluid meniscus given by

$r{}_{ij}=\frac{p_{ij}^C}{{\sigma }^{nw}},$

θ is the contact angle of the liquid phase, and σ^nw is the interfacial tension between the gas phase and the liquid phase.

A perimeter of the gas phase region

$C_{ij}^{nw}$

is given by:

$C_{ij}^{nw}=8\left[R_{ij}+r_{ij}\left(-\frac{\sin \left(\frac{\pi }{4}-\theta \right)}{\sin \frac{\pi }{4}}+\left(\frac{\pi }{4}-\theta \right)\right)\right]$

The conductivity of the liquid phase in the corners of the throat may be determined by: (Chen, C. Qin, B. Guo, Fully Implicit Dynamic Pore-Network Modeling of Two-Phase Flow and Phase Change in Porous Media. 2020):

$K_{ij}^w=4\frac{2{\tilde{g}}_{ij}{l}_{ij}^4}{{\mu }_{ij}^w{L}_{ij}}$

$l_{ij}=r_{ij}\frac{\cos \left(\theta +\beta \right)}{\sin \beta }$

$\begin{array}{l}{\tilde{g}}_{ij}=\\ exp\left\{\frac{\left[-18.2066{\left({\tilde{G}}_{ij}^w\right)}^2+5.88287{\tilde{G}}_{ij}^w-0.351908+0.02\sin \left(\beta -\frac{\pi }{\sigma }\right)\right]}{\left(\frac{1}{4\pi }-{\tilde{G}}_{ij}^w\right)}+2ln{\tilde{A}}_{ij}^w\right)\end{array}$

${\tilde{A}}_{ij}^w={\left[\frac{\sin \beta }{\cos \left(\theta +\beta \right)}\right]}^2\left[\frac{\sin \beta \cos \left(\theta +\beta \right)}{\sin \beta }+\theta +\beta -\frac{\pi }{2}\right]$

${\tilde{G}}_{ij}^w=\frac{{\tilde{A}}_{ij}^w}{4{\left[1-\left(\theta +\beta -\frac{\pi }{2}\right)\sin \left(\beta \right)/\cos \left(\theta +\beta \right)\right]}^2}$

where

$K_{ij}^w$

is the conductivity of the liquid phase for the throat between pores i and j, g̃_ij is a dimensionless conductivity of the liquid phase at a half corner angle β, β is

$\frac{\pi }{4}$

fora square-tube pore throat,

${\tilde{G}}_{ij}^w$

is a shape factor of the liquid phase at a corner with a unit meniscus-apex distance,

${\tilde{A}}_{ij}^w$

is an area of the liquid phase region at a corner with a unit meniscus-apex distance, l_ij is a meniscus-apex distance along the wall at a corner.

Based on the above, the phase conductivities of the liquid phase and the gas phase in the throats can be determined. After steps 108 and 110 are completed, conditions in the throat are known, including throat occupancy as well as phase conductivities of the liquid phase and the gas phase.

In steps 112 to 122, the molar balance for each component, including the solvent and hydrocarbon components, in the pore network model are determined for a given time. In determining the molar balance for each component, conservation of mass within the pore network model is considered, which may then be used to examine mass transfer of the components within the pore network model. As the molar balance for each component can be represented by a series of nonlinear equations, an implicit numerical analysis method may be used to provide more stable computations in determining the molar balance for each component. The series of nonlinear equations that represent the molar balance for each component may be solved using an iterative method, such as the Newton method, to determine the molar balance for each component at a given time as such a method may be computationally more advantageous with faster time to a solution and reduced processor computations. Steps 112 to 122 are performed for each iteration in the Newton method. The result of all iterations of steps 112 to 122 is a completed molar balance for each component for a given time, starting from a time t. As time progresses in the flow 100, the multiple iterations of steps 112 to 122 are performed at each time and thus, the molar balance of each component is determined at each time. That is, each time in the flow 100 will include multiple iterations of steps 112 to 122. Time progresses in the flow 100 by a time step Δt that may be selected based on what is being modelled within the flow 100. The time step Δt may change based on results either between iterations within the iterative method of steps 112 to 122 or based on results from one time to a previous t-Δt (or subsequent t+Δt) time. By tracking a change in the molar balance of each component in each pore in the pore network model over time the mass transfer for each component over time can be traced.

In step 112, a representation for the molar balance for solvent and hydrocarbon components in the pore network model is established.

For each pore: i = 1, N_pore, the molar balance equation for the k components (where k may be hydrocarbons, solvent, or other components in the subterranean reservoir that are to be considered) to account for mass conservation can be represented by the residual vector F as determined by:

$\begin{array}{cc}\begin{array}{l}\left[F_i^k\right]i=1,N_{pore}=\\ V_i\frac{1}{\Delta \text{t}}\left[\left({\rho }_i^i{x}_i^k{S}_i^l+{\rho }_i^g{y}_i^k{S}_i^g\right)-{\left({\rho }_i^i{x}_i^k{S}_i^l+{\rho }_i^g{y}_i^k{S}_i^g\right)}^0\right]\\ +{\displaystyle {\sum }_{j=1}^{N_i}{\rho }_{ij}^l{x}_{ij}^k{K}_{ij}^l\left(P_i^l-P_j^l\right)}\\ +{\displaystyle {\sum }_{j=1}^{N_i}{\rho }_{ij}^g{y}_{ij}^k{K}_{ij}^g\left(P_i^g-P_j^g\right)}\\ +{\displaystyle {\sum }_{j=1}^{N_i}\frac{D_{ij}^l{A}_{ij}^l}{l_{ij}}}\left[{\rho }_i^l\max \left(x_i^k-x_j^k,0\right)+{\rho }_j^i\min \left(x_i^k-x_j^k,0\right)\right]\\ +{\displaystyle {\sum }_{j=1}^{N_i}\frac{D_{ij}^g{A}_{ij}^g}{l_{ij}}\left[{\rho }_i^g\max \left(y_i^k-y_j^k,0\right)+{\rho }_j^g\min \left(y_i^k-y_j^k,0\right)\right]=0}\end{array}& \text{­­­(1)}\end{array}$

where V_i is the pore volume, N_i is the number of pore throats connected to pore body

$i,x_i^k,y_i^k$

are the component k (solvent or hydrocarbon or other components) mole fraction in liquid and gas phases in pore

$\text{i,}{\text{S}}_i^l\text{and}{S}_i^g$

are the liquid phase saturation and gas phase saturation, respectively, in pore i, ρ^l, ρ^g are the liquid and gas molar densities (i.e., the sum of the molar densities of all components in the liquid or gas phase),

$D_{ij}^{l,g}$

is the diffusion coefficient in the liquid and gas phases between pores i and j,

$j\text{,}{P}_i^l\text{and}{P}_i^g$

are the liquid and gas pressures in pore i,

$K_{ij}^l\text{and}{K}_{ij}^g$

are the conductivities of liquid and gas phases in a throat connecting pores i and j,

$A_{ij}^{\alpha }$

is the cross sectional area of phase α in the throat connecting pores i and j

$x_{ij}^k,y_{ij}^k$

are the component k (solvent or hydrocarbon or other components ) mole fraction in liquid and gasphases in the throat connecting pores i and j, V_i is the volume of the pore i, and L_ij is the length of the throat between pores i and j.

The relationship between liquid and gas pressure in each pore may be determined by

$p_i^c\left(S_i^l\right)=P_i^g-P_i^l$

where

$p_i^c$

is the capillary pressure in a pore i. This may be used to determine the pressure of the gas phase given the saturation and pressure of the liquid phase.

The mole fractions in each phase satisfy

${\displaystyle {\sum }_{k=1}^{N_c}{x}_i^k=1},{\displaystyle {\sum }_{k=1}^{N_c}{y}_i^k=1}.$

The total volume of fluid fills each pore volume, therefore the liquid and gas saturations satisfies F_S = S^l + S^g - 1 = 0.

The molar balance of the solvent and hydrocarbon equations above are solved where the residual moles of the solvent and hydrocarbon within each pore are represented by F =

$\left[F_i^k\right]i=1,N_{pore}$

where k is the components in the subterranean reservoir.

In step 114, the rate of change of the molar balance of each component in each pore with respect to the first set of characteristics of the pore is determined.

A rate of change of the molar balance of each component in a pore i with respect to a rate of change of the first set of characteristics for a pore j is determined. This may be represented by a matrix, for example a Jacobian matrix that can be used in Newton iterations, given by

$Jac=\frac{\text{δ}{F}_i}{\text{δ}{X}_j}$

where the X represents the first set of characteristics as noted above, in particular,

$X={\left[S_i^l{P}_i^l\right]}_{i=1,N_{pore}}$

for a two phase pore and

$X={\left[x_i^{k=solv}{P}_i^l\right]}_{i=1,N_{pore}}$

for a single phase pore.

Using a Jacobian matrix with the Newton iterations can provide some computational advantages. Such a combination can result in reduced computational cost by providing a faster convergence of the solution of the molar balance of all components. Further, the properties of the Jacobian matrix may be exploited in multi-processor computer systems to provide for parallel processing.

In step 116, a change in the first set of characteristics for a pore, represented by X for the pore, is determined. The change in the first set of characteristics for the pore is determined using the previously determined molar balance and the previously determined change of the molar balance in with respect to the first set of characteristics. The change in the first set of characteristics of the pore (ΔX) is determined by Jac × ΔX = -F.

In step 118, the first set of characteristics for the pore is updated based on the change in the first set of characteristics for the pore (ΔX) that was determined in step 116. For example, the first set of characteristics for the pore may be updated based on X^t+Δt,iter+1 = X^t+Δt,iter + ΔX where iter represents the current iteration and t represents the current time. As step 118 is part of the iterative method used to determine the molar balance of the components at a given time, the first set of characteristics is updated for each iteration in the iterative method.

In step 120, phase conditions for each pore are determined. The phase conditions may include phase occupancy, that is whether a pore contains a single phase or two phases. In order to determine whether a pore is in a single phase or in two phases, the Rachford-Rice equation may be solved for the mole fraction of the gas phase β,

${\displaystyle {\sum }_{i=1}^{N_c}\frac{z_i\left(K_i-1\right)}{1-\beta +\beta K_i}}=0$

where

$K_i=\frac{y_i}{x_i}$

is an equilibrium factor for component i. x_i, y_i are the mole fraction of x_i components i in the gas and liquid phases, respectively. z_i is the overall mole fraction of component i in the pore. The above equation is solved for the variable mole fraction of the gas phase (β) assuming a constant K_i. If the result is negative, or β is less than zero, then there is no gas phase in the pore.

In the case where the pore network model represents a subterranean reservoir containing hydrocarbons into which a solvent is injected, there are some simplifications that may be made to reduce complexity in step 120. For example, an assumption might be made that the heavy hydrocarbons are not in the gas phase (particularly if the hydrocarbons are bitumen). Additionally, the solubility of the solvent injected into the subterranean reservoir may be considered to be constant if the numerical domain is small enough and the pressure does not change significantly.

In step 122, the phase saturation is determined for each pore. In the case where the pore network model represents a hydrocarbon-containing subterranean reservoir, the solubility of solvent in the hydrocarbon may be constant and independent of the variations of the pressure. When the hydrocarbon is bitumen, it may be present only in the liquid phase where the solvent is present in both gas and liquid phases.

As noted above, steps 108 to 122 are iterative and are performed for a given time. Decision points at steps 124, 126 or 128 assess the results of the iterative steps and whether those results are as expected based on a selected time step Δt between times. If a time step is selected that does not produce desired results, the time step Δt may change in steps 124 or 128. Steps 124 and 128 assess whether changes are occurring too quickly between iterations and between time. It may be desirable to control the amount of change as there may be multiple events occurring successively over a short period of time which may impact whether a molar balance result is able to be found or may result in divergence if the changes between iterations or between time are outside of an acceptable range. Step 126 assesses whether or not there has been a change from a previous iteration, or how much of a change has occurred.

In step 124, a maximum change in pore liquid saturation between consecutive iterations is used to determine if the time step Δt is to be reduced and the iteration is to be performed again. The maximum change in pore liquid saturation with respect the previous iteration may be determined by:

$\Delta S^{iter}=max{\left|\frac{S_i^{iter+1}-S_i^{iter}}{S_i^{iter}}\right|}_{i=1,N_{pore}}$

If the maximum change in pore saturation between iterations ΔS^iter exceeds an allowed value then this may represent a large change between iterations. Such large changes may reflect that assumptions made within the pore network model may be incorrect or invalid. To reduce the amount of change, a smaller time step may be determined. If a smaller time step is to be used then the time step Δt will be reduced and the iterative process will be repeated from step 108 with the revised time step.

The amount of maximum change in pore saturation between iterations ΔS^iter that is considered desirable or allowable will depend on the system being modeled among other factors. For example, a 20% change may be considered to be the maximum allowable change. In other situations the maximum allowance change may be different.

If the amount of maximum change in pore saturation between iterations ΔS^iter is considered to be acceptable then in step 126 the relative change of X (first set of characteristics of each pore) between iterations is examined. If the change in the first set of characteristics between iterations is below a certain amount, then this may indicate that there is little or no change in the first set of characteristics of each pore. The Newton method is considered to have converged, and the molar balance of the components is considered to have been determined, when there is little or no change in the first set of characteristics of each pore between iterations. In this case, the flow 100 continues to step 128. However, if the change between iterations of the first set of characteristics of each pore exceeds a given amount, this indicates that there is still change and the molar balance of the components as determined by the Newton method has not been determined (and the Newton method has not converged). In this case, the flow 100 returns to step 108 for the next iteration.

In step 128, the maximum change in pore saturation between consecutive time steps is used to determine if the time step is to be reduced and the process is to be performed again from step 108.

Saturation in each pore

${\left[S_i^{t+\Delta t}\right]}_{i=1,N_{pore}}$

may be updated after it is determined in step 126 that the molar balance for each component has been determined for the given time. After the saturation in the pores has been updated, a maximum change in pore saturation with respect to a previous time step is determined:

$\Delta S^t=max{\left|\frac{S_i^{t+\Delta t}-S_i^t}{S_i^i}\right|}_{i=1,N_{pore}}$

As with step 124, in step 128, if the maximum change in pore saturation between time steps ΔS^t exceeds an allowed value, then this may represent a large change between time steps. Such large changes may reflect that assumptions made within the pore network model may be incorrect or invalid. To reduce the amount of change, a smaller time step may be determined. If a smaller time step is to be used then the time step Δt will be reduced and the process will be repeated from step 108 with the revised time step.

The amount of maximum change in pore saturation between iterations ΔS^t that is considered desirable or allowable will depend on the system being modeled among other factors. For example, a 20% change may be considered to be the maximum allowable change. In other situations the maximum allowance change may be different.

In step 130, the next time step is estimated for the next time to be used for flow 100. This next time step may be based on a rate of local liquid drainage or imbibition in pores and a rate of change of solvent mole fraction or the rate of change of solvent mole fraction.

Local liquid drainage may be determined based on a change in saturation of a pore. For example, if liquid phase saturation of a pore is decreasing then the pore is draining and liquid is coming out of the pore. If liquid is coming out of the pore then a time step may be kept small enough to avoid unrealistic values of saturation in the updated time step. Similarly, if liquid phase saturation of a pore is increasing, then the pore is imbibed with liquid and gas is coming out of the pore. The time step should be small enough to avoid the saturation value exceeding 1.

The next time step Δt may be selected to be a minimum of an estimated local pore time step, At_i. The local time step to pore i, Δt_i, may be proportional to

$\left(1-S_i^l\right)$

or to

$\left(S_i^l-S_i^{lmin}\right)\text{where}{S}_i^l$

is the liquid saturation. (Journal of Fluid Mechanics (2010), Joekar-Niasar et al., Vol. 655, page 38-71, S.Chen, C. Qin, B. Guo, Fully Implicit Dynamic Pore-Network Modeling of Two-Phase Flow and Phase Change in Porous Media. 2020). For example, for a pore in imbibition the local time step may be determined by

$\Delta t_i=\frac{V_i}{q_n\left(i\right)}\left(1-S_i^l\right)$

whereas for a pore under drainage the local time step may be determined by

$\Delta t_i=\frac{V_i}{q_n\left(i\right)}\left(S_i^l-S_i^{lmin}\right)$

where q_n(i) is the net volumetric gas flow rate to pore i.

In step 132, the process from step 108 is repeated at a new time t+Δt until a final time is reached. The time selected should be long enough so that all, or at least most, pores become invaded with the gas phase. This may depend on how fast hydrocarbons are produced in the subterranean reservoir.

Extraction of hydrocarbons during the various times in flow 100 may be considered by assessing the change in the molar balance of the hydrocarbons in the pore network model. Any difference in the molar balance of hydrocarbons between different times is considered to be hydrocarbons that are extracted or produced by the production well.

In step 134, the results from the pore network are used to adjust operational parameters of wells within the subterranean reservoir. The operational parameters of the solvent injection may be adjusted in a constant or incremental manner to examine the impact of changes of a particular operational parameter(s) on the extraction of hydrocarbons. A flow pattern (e.g., viscous fingering or capillary fingering) in the subterranean reservoir that is represented by the pore network model may also be determined and used to adjust the operational parameters of the solvent injection. If a high flow rate for the solvent is used, a thin diluted layer of hydrocarbons at the interface between the solvent and the hydrocarbons in the subterranean reservoir may be moved away from the interface to expose fresh hydrocarbons to solvent. Such refreshing of this interface improves mass transfer of the hydrocarbons in the subterranean reservoir. For example, when the solvent flow rate is high enough, a thin diluted layer at the interface between the solvent and the hydrocarbons is moved away from the interface to expose fresh hydrocarbons to solvent. This results in mass transfer enhancement under conditions of the viscous fingering flow pattern, which tend to occur at higher injection rates. When solvent flow rate is low, solvent is transferred slowly into the hydrocarbons and the thickness of the layer of diluted hydrocarbons at the interface between solvent and hydrocarbon increases with time. A diluted hydrocarbon layer remains immobile until the capillary pore pressure exceeds the capillary entry pressure at the throat and the throat becomes invaded by capillary fingering flow pattern, which tends to occur when viscous forces are small compared to capillary forces. A capillary number (see for example, S. An, H. Erfani, O. E. Godinez-Brizuela, Vahid Niasar, “Transition from Viscous Fingering to Capillary Fingering: Application of GPU-Based Fully Implicit Dynamic Pore Network Modeling,” Water Resources Research, 10.1029/2020WR028149), which is a ratio of viscous force to capillary force, and the viscosity ratio, which is the viscosity ratio of the invading phase over the receding phase, may be used to determine whether there is a capillary fingering pattern or a viscous fingering pattern present in the pore network model. Based on whether there is a capillary fingering pattern or a viscous fingering pattern present, the operational parameters may be adjusted to obtain the desired pattern.

One such operational parameter that may be adjusted is a flow rate of the solvent. The flow rate of the solvent may be adjusted in a stepwise manner with changing increments (e.g., increasing or decreasing) or may be adjusted based on other factors.

Another such operational parameter that can be adjusted is the pressure difference between an injection well that injects solvent into the subterranean reservoir and a production well that extracts hydrocarbons that have been mobilized by the solvent from the subterranean reservoir. That is, the pore network model can also be used to design experimental and pilot scale set-ups, select distance between injector and producer wells, interpret results and understand the impact of scale up on the pore scale two phase flow pattern. Such an operating parameter may be adjusted prior to drilling of the wells into the subterranean reservoir based on results of the pore network model using conditions in the subterranean reservoir and performed prior to such drilling. The wells may then be drilled based on the results of the pore network model.

Another operational parameter that can be adjusted is gas pressure of the injected solvent. The gas pressure may also be adjusted so that there is a viscous fingering pattern or a capillary fingering pattern.

Further, different solvents or different additives to the solvents can be used to adjust operational parameters of the injected solvents in the flow 100 to determine improved extraction of the hydrocarbons.

FIG. 4 illustrates an alternate embodiment of an exemplary flow 300 for a computer system of computer modelling of a multiphase extraction process using a dynamic pore network model. The flow 300 in FIG. 4 uses a first set of characteristics to define a state of a pore based on specified constant parameters representing conditions used in thermodynamic calculations.. The pore network model that is developed and used in flow 300 is based on the constant parameters of constant volume and temperature for thermodynamic flash calculations to determine the composition of phases in equilibrium.

Information describing a subterranean reservoir is obtained in step 302. This information may include information on a structure of pores, matrix, and throats of the subterranean reservoir as well as conditions in the subterranean reservoir. Similar to step 102, the information describing the subterranean reservoir may be randomly generated, determined using known properties of a subterranean reservoir or may be obtained through measurements of a subterranean reservoir that are obtained by a capture device.

A structure for pore network model of a subterranean reservoir is generated by the computer system in step 304 based on the information about the structure of the subterranean reservoir that was obtained in step 302 in a manner similar to generation of the pore network model in step 102. The structure for the pore network model represents the structure of the subterranean reservoir including pores that are interconnected by throats.

In step 306, initial conditions in the pores and throats in the pore network model are established. These initial conditions may include pore saturation and solvent mole fraction based on the conditions in the subterranean reservoir that are obtained in step 302. For example, if wells are either present or are being drilled in the subterranean reservoir to be modeled, pore saturation may be selected based on existing conditions that have been determined or estimated through known mechanisms. Likewise, the solvent mole fraction used for the conditions in the subterranean reservoir may be based on known solvent mole fractions that are being used for similar purposes. The initial conditions for the pore network model may also include the phase conditions for the pore that establish whether the pore is in one phase (i.e., liquid or gas) or two phases (i.e., liquid and gas).

In step 308, the operating parameters for injection of the solvent are established. These will be the initial operating parameters and may include at least some of the following: placement of an injection well within the structure of the pore network model, a flow rate for the solvent, a gas pressure for the injected solvent, initial pressure in the pores and pressure at the production well, solvent composition, and solvent additives.

In step 310, a first set of characteristics for each pore in the pore network model is defined. The first set of characteristics may be characteristics or properties that are either measured or obtained from conservation equations. A second set of characteristics for each pore will be determined at a later point as the second set of characteristics will be those characteristics or properties that can be derived from the first set of characteristics.

In the flow 300 of FIG. 4, the first set of characteristics for each pore (whether the pore is in one phase or two phases) is

$X_p=\left[{\left\{C_i^{k=1,N_C}\right\}}_{i=1,Npore}\right]\text{where}{C}_i^k$

is the molar density of component k in all phases, N_c is the number of components within the subterranean reservoir that are being considered in the pore network model, and N_pore is the number of pores in the pore network model. The overall density of the components k in a pore is determined by

$C_i^k=\left(C_{gi}^k{S}_i^g+C_{li}^k{S}_i^l\right)\text{where}{C}_{\infty i}^k$

is the molar density in pore i of component k in phase α (i.e., gas (g) or liquid (l)), S^α is the saturation for the phase α, k is the component (i.e., solvent or hydrocarbon).

In step 312, the phase conditions of each pore are determined based on constant temperature and volume. The phase conditions of each pore describe whether a pore is in a single phase or two phases. Stability analysis based on tangent plane distance criterion can be used to determine whether a pore is in a single phase or two phases. For this, a minimum of the tangent plane distance function is determined based on overall components molar density in the pores at constant temperature and volume. If the function at this minimum is negative, it means that the phase in the pore is unstable and there is more than one phase in the pore.

The hydrocarbons that are to be extracted from the subterranean reservoir may be considered a mixture of components and, as such, phase equilibrium calculations based on an equation of state should be performed to obtain the composition of the phases in the pores.

The condition of thermodynamic equilibrium between liquid and gas phases in each pore may be determined by

$log{\widehat{f}}_{k,i}^g-log{\widehat{f}}_{k,i}^l=0,k=1,N_c,i=1,N_{pore}\text{where}{\widehat{f}}_{k,i}^{\alpha }$

is thefugacity of component k in phase α in pore i. Phase stability and two-phase flash calculations may be performed to obtain the compositions of phases in equilibrium in the pores.

In step 314, a second set of characteristics for each pore are derived from the first set of characteristics by solving equation Ψ = 0, presented below based on the phase conditions for the pore. For each pore that contains two phases (e.g., a two phase pore), as determined by the phase conditions in step 312, the second set of characteristics may be derived from the first set of characteristics. The second set of characteristics may include component molar density in each phase and liquid saturation and may be represented by

$X_{SV}=$

$\left[{\left\{C_{li}^{k=1,N_c},C_{gi}^{k=1,N_c},S_i^l\right\}}_{i=1,N_{pore}}\right]\text{where}{C}_{li}^{k=1,N_c}$

is the molar density of the component k in the liquid phase in pore i and

$C_{gi}^{k=1,N_c}$

is the molar density of the component k in the gas phase in pore i. The second set of characteristics for each pore may be determined using thermodynamic flash calculation based on constant temperature and constant volume using the first set of characteristics. This provides a representation of what is happening in a pore with a description of each phase in a pore.

The relationship between the first set of characteristics and the second set of characteristics may be given by a system of (2N_c + 1)N_pore equations:

${}^{\text{Ψ=}\left\{\begin{array}{c}{\text{Ψ}}_{F_i}^k={\left\{\log \left({\widehat{f}}_k^g\right)-\log \left({\widehat{f}}_k^l\right)\right\}}_i=0\\ {\text{Ψ}}_{C_i}^k=\left\{C^k-\left(C_g^k\left(1-S^L\right)\right)+C_l^l{S}^l\right\}=0\\ {\text{Ψ}}_{P_i}={\left\{\frac{1}{RT}\left(-P^g+P^l+P^c\right)\right\}}_i=0\end{array}\right)}$

where pore i varies from 1 to N_pore and component k varies from 1 to N_c in the first two equations of the above,

$P_i^C$

is the capillary pressure in pore i and is a function of the liquid phase saturation, phase pressures P^g and P^l can be determined using the equation of state such as, for example, cubic equations of state such as Peng-Robinson equation of state or SAFT like equation of state.

The implicit function theorem gives:

$\frac{\partial X_{sv}}{\partial X_P}=-{\left(\frac{\partial \text{Ψ}}{\partial X_{SV}}\right)}^{-1}\frac{\partial \text{Ψ}}{\partial X_P}$

The matrix representing the rate of change (or the derivative) of the secondary set of characteristics with respect to the first set of characteristics is obtained by

$\frac{\partial X_{sv}}{\partial X_P},$

which is used later in step 322.

In steps 316 and 318, phase conditions in each throat are determined based on connected pores using pore capillary pressures, throat invasion and snap-off events. The phase conditions may include phase occupancy in the throat representing which phases (e.g., liquid phase and/or gas phase) are present in the throat and phase conductivities in the throat. The phase occupancy in each throat is determined in step 316 may be determined in a manner similar to the phase occupancy in each throat that is determined in step 108 for flow 100.

In step 318, the phase conductivities of the liquid and gas phases in the throats are determined. The phase conductivities of the liquid and gas phases in the throats in flow 300 may be determined in a manner similar to the determination made in step 110 of flow 100.

In steps 320 to 324, the molar balance for each component, including the solvent and hydrocarbon components, in the pore network model are determined. In determining the molar balance for each component, conservation of mass within the pore network model is considered, which may then be used to examine mass transfer of the components within the pore network model. Steps 312 to 326 are similar to steps 108 to 122 in that these steps are performed for each iteration in the Newton method. The result of all iterations of steps 312 to 326 is a completed molar balance for each component for a given time, starting from a time t. As time progresses in the flow 300, the multiple iterations of steps 312 to 326 are performed at each time and thus, the molar balance of each component is determined at each time with time progressing by a time step Δt that may change based on results either between iterations within the iterative method of steps 312 to 326 or based on results from one time to a previous t-Δt (or subsequent t+Δt) time. By tracking a change in the molar balance of each component in each pore in the pore network model over time the mass transfer for each component over time can be traced.

In step 320, a representation for the molar balance for solvent and hydrocarbon components in the pore network model is established.

At each pore i = 1, N_pore, the molar balance equation for the k components (where k may be hydrocarbons, solvent, or other components in the subterranean reservoir that are to be considered) to account for mass conservation can be represented by the residual vector F as determined by:

$\begin{array}{l}{F}_i^k=V_t\frac{1}{\Delta \text{t}}\left[\left(C_{li}^k{S}_i^l+C_{gi}^k{S}_i^g\right)-{\left(C_{li}^k{S}_i^l+C_{gi}^k{S}_i^g\right)}^0\right]\\ +{\displaystyle {\sum }_{j=1}^{N_i}{C}_{lij}^k}{K}_{ij}^l\left(P_i^l-P_j^l\right)\\ +{\displaystyle {\sum }_{j=1}^{N_i}{C}_{gij}^k}{K}_{ij}^g\left(P_i^g-P_j^g\right)\\ +{\displaystyle {\sum }_{j=1}^{N_i}\frac{D_{ij}^l{A}_{ij}^l}{l_{ij}}\left[{\text{ρ}}_i^l\max \left(x_i^k-x_j^k,0\right)+{\text{ρ}}_j^l\min \left(x_i^k-x_j^k,0\right)\right]}\\ +{\displaystyle {\sum }_{j=1}^{N_i}\frac{D_{ij}^g{A}_{ij}^g}{l_{ij}}\left[{\text{ρ}}_i^g\max \left(y_i^k-y_j^k,0\right)+{\text{ρ}}_j^g\min \left(y_i^k-y_j^k,0\right)\right]}=0\end{array}$

The residual moles of the solvent and hydrocarbon within each pore are represented by F =

In step 322, a change in the molar balance for all components with respect to the first set of characteristics X_P is determined.

When the Newton method is used, a rate of change of the molar balance in each pore with respect to the first set of characteristics for the pore in flow 300 may be represented by a matrix, for example a Jacobian matrix used in newton iterations, given by

$Jac=\left[\frac{\text{δ}F}{\text{δ}{X}_P}\right]={\left[\frac{\text{δ}F}{\text{δ}{X}_P}\right]}_{X_{SV}}+{\left[\frac{\text{δ}F}{\text{δ}{X}_{SV}}\right]}_{X_P}\left[\frac{\text{δ}{X}_{SV}}{\text{δ}{X}_P}\right]$

The Jacobian matrix above will be a [N_pore*(N_c+1)] x [N_pore*(N_c+1)] matrix where N_pore is the number of pores in the pore network model and N_c is the number of components being considered in the pore network model. The residual vector F will be a [N_pore*(N_c+1)] vector. The Jacobian matrix may be assembled using submatrices

$\frac{\text{δ}{F}_i}{\text{δ}{X}_j}$

with the only non-zero submatrices being when pores i and j are connected by a throat. This results in a submatrix that is [N_c+1] × [N_c+1] that corresponds to the derivative of the residual vector in pore i with respect to primary characteristics of pore j.

In step 324, a change in the first set of characteristics for a pore, represented by X_P for the pore, is determined. The change in the first set of characteristics for the pore is determined using the previously determined molar balance and the previous determined change of the molar balance in with respect to the first set of characteristics for the pore. The change in the first set of characteristics of the pore (ΔX) is determined by solving the equation: Jac × ΔX = -F.

Separation of characteristics of each pore into the first set of characteristics and the second set of characteristics that can be derived from the first set of characteristics reduces the computational complexity of computer modelling of the hydrocarbon extraction process using the pore network model. Since the second set of characteristics can be derived from the first set of characteristics, the Newton method is applied to solve only for the primary set of characteristics. As such, the Newton method is solved only for the first set of characteristics which reduces the size of the Jacobian matrix and the computational complexity. This provides a computational advantage that enables molar balance for each component to be determined and thus enables mass transfer within the pore network model to be traced.

In step 326, the first set of characteristics is updated based on the change in the first set of characteristics for the pore (ΔX_p) that was previously determined. For example, the first set of characteristics for the pore may be updated based on X^t+Δt,iter+1 = X^t+Δt,iter + ΔX. This updating of the first set of characteristics for a pore in step 326 of flow 300 may be done in a manner similar to the determination made in step 118 of flow 100.

As noted above, steps 312 to 326 are iterative and are performed for a given time. Decision points at steps 328, 330 and 332 assess the results of the iterative steps and whether those results are as expected based on a selected time step Δt between times. If a time step is selected that does not produce desired results, the time step Δt may change in steps 328 or 332.

In step 328, a change in a predetermined characteristic is examined with respect to the value of that predetermined characteristics from a previous consecutive iteration. For example, the predetermined characteristic may be a maximum change in pore saturation, which is compared between consecutive iterations and is used to determine if the time step Δt is to be reduced and the flow 300 is to be performed again with this revised time step. That is, when the predetermined characteristic is a maximum change in pore saturation, the maximum change in pore saturation at the current iteration is compared to a maximum change in pore saturation from the previous iteration. The determination of a change in predetermined characteristic with respect to the previous iteration in step 328 of flow 300 may be performed in a manner similar to step 124 of flow 100.

If the amount of change in the predetermined characteristic between iterations is considered to be acceptable then in step 330, then the relative change of the first set of characteristics of each pore between iterations is examined. If the change in the first set of characteristics between iterations is below a certain amount, then this may indicate that there is little or no change in the first set of characteristics of each pore. The Newton method is considered to have converged, and the molar balance of the components is considered to have been determined, when there is little or no change in the first set of characteristics of each pore between iterations. In this case, the flow 300 continues to step 332. However, if the change between iterations of the first set of characteristics of each pore exceeds a given amount, this indicates that there is still change and the molar balance of the components as determined by the Newton method has not been determined (and the Newton method has not converged). In this case, the flow 300 returns to step 312 for the next iteration. Step 330 of flow 300 is similar to step 126 of flow 100.

In step 332, the change in a predetermined characteristic between consecutive times is assessed. The same predetermined characteristic considered in step 328 with respect to previous iterations may also be considered in step 332. However, different predetermined characteristics may also be considered for step 328 versus 332 depending on what is being modeled. In general, step 328 will check for a change in a predetermined characteristic between consecutive iterations of the flow 300 while step 332 will check for a change in a predetermined characteristic between consecutive times of the flow 300 (where each time in the flow 300 may include multiple iterations). For example, if the predetermined characteristic being considered in step 332 is a maximum change in pore saturation, then a maximum change in pore saturation is compared to a maximum change in pore saturation from the previous time step. If the change in the predetermined characteristics between different times (e.g. t versus t+Δt), exceeds a tolerance, then this may indicate that assumptions made may be incorrect or that the characteristics are changing too fast to illustrate the changes happening in the subterranean reservoir. Determination of a change in the predetermined characteristic compared to the previous time step in step 332 of flow 300 may be performed in a manner similar to step 128 of flow 100.

In step 334, the next time step is estimated for the next time to be used for flow 300. This next time step may be based on, for example, a rate of local liquid drainage or imbibition in pores. The next time step determined in step 334 of flow 300 may be determined in a manner similar to the determination made in step 130 of flow 100.

In step 336, the process from step 312 is repeated until a time ends. The time selected should be long enough so that all, or at least most, pores become invaded with the gas phase. This may depend on how fast hydrocarbons are produced in the subterranean reservoir. The end criteria used in step 336 of flow 300 may be similar to the end criteria used in step 132 of flow 100.

In step 338, the extraction of hydrocarbons during the various times in flow 300 may be considered by assessing the change in the molar balance of the hydrocarbons in the pore network model. Any difference in the molar balance of hydrocarbons between different times is considered to be hydrocarbons that are extracted or produced by the production well. While step 338 is illustrated as occurring after the time has ended, step 338 may also be performed within every time step before the time changes.

In step 340, the results from the pore network model are used to adjust operational parameters of wells within the subterranean reservoir. The application of the results from the pore network model in step 340 of flow 300 may be similar to the application in step 134 of flow 100.

FIG. 5 illustrates an alternate embodiment of an exemplary flow 400 for a computer system of computer modeling of a multiphase extraction process using a dynamic pore network model. The flow 400 in FIG. 5 uses a different first set of characteristics to define a state of a pore. The flow 400 in FIG. 5 uses a first set of characteristics to define a state of a pore based on specified constant parameters representing conditions used in thermodynamic calculations.. The pore network model that is developed and used in flow 400 is based on the constant parameters of constant temperature and constant pressure for thermodynamic phase stability and flash calculations.

Steps 402 to 408 in flow 400 may be performed in a manner similar to steps 302 to 308 of flow 400. Information describing a subterranean reservoir is obtained in step 402. A structure for pore network model of a subterranean reservoir is generated by a computer system in step 404 based on the information about the structure of the subterranean reservoir that was obtained in step 402. In step 406, initial conditions in the pores and throats in the pore network model are established. In step 408, the operating parameters for injection of the solvent are established.

In step 410, a first set of characteristics for each pore in the pore network model is defined. The first set of characteristics may be characteristics or properties that are either measured or obtained from conservation equations. A second set of characteristics for each pore will be determined at a later point as the second set of characteristics will be those characteristics or properties that can be derived from the first set of characteristics.

In the flow 400 of FIG. 5, the first set of characteristics for each pore are X_p = [X_Pi=1,Npore], where

$X_{pi}=\left\{C_i^{k=1,N_c},P_i^{l,g}\right\}$

for two phase pore (i.e., both liquid phase and gas phase). For single phase pores the first set of characteristics can be

$X_{pi}=\left\{x_i^{k=1,N_c},P_i^l\right\}$

for single phase pore in liquid phase or

$X_{pi}=\left\{y_i^{k=1,N_c},P_i^g\right\}$

for a single phase pore in gas phase, where

$p_i^{lorg}$

is the pressure in the gas or liquid phase (depending on the phase conditions) forpore i,

$x_i^{k=1,N_c}$

is the mole fraction of the component k in liquid phase,

$y_i^{k=1,N_c}$

is the mole fraction of the component k in the gas phase,

$C_i^k$

is the overall molar density of component k, N_c is the number of components within the subterranean reservoir that are being considered in the pore network model, and N_pore is the number of pores in the pore network model. Alternatively, the first set of characteristics for a pore in a single phase may be:

$X_{pi}=$

$\left\{C_i^{k=1,N_c},P_i^{lorg}\right\}\text{as}{C}_i^k\text{and}{x}_i^{k=1,N_c}\text{or}{y}_i^{k=1,N_c}$

or related such that x_i (or y_i) for a component k is determined by C_i, the overall molar density for a component k divided by

$C_i^{allcomponents}$

the overall molar density for all components in the pore network model. The overall molar density of the components k in a pore i is determined by:

$\left\{C_i^k=\left({\rho }_i^g{S}_i^g{y}_i^k+{\rho }_i^l{S}_i^l{x}_i^k\right)\right\}$

where ρ^α and S^α are the molar density and saturation for the phase α (i.e., gas (g) or liquid (l)), k is the components (i.e., solvent or hydrocarbon), y^k is mole fraction of component k in the gas phase, and x^k is the mole fraction of component k in the liquid phase. The above pressure

$P_i^{l,g}\text{in}$

the first set of characteristics can be either liquid or gas pressure.

In step 412, the phase conditions of each pore are determined based on constant temperature and pressure. The phase conditions of each pore describe whether a pore is in a single phase or two phases. Although the conditions of constant temperature and pressure in flow 400 are different from the conditions of flow 300, step 412 is performed in a manner similar to step 312 but using the tangent plane distance function at constant temperature and pressure.

In step 414, a second set of characteristics for the pore are derived from the first set of characteristics based on the phase conditions for the pore. For each pore that contains two phases (e.g., a two phase pore), the second set of characteristics may be derived from the first set of characteristics. The second set of characteristics in flow 400 may include component mole fraction in each phase and phase saturations and may be represented by

$X_{SV}=$

$\left[{\left\{x_i^{k=1,N_c},y_i^{k=1,N_c},S_i^l,S_i^g\right\}}_{i=1,N_{pore}}\right].$

The second set of characteristics for each pore may be determined from the first set of characteristics using equations based on constant temperatureand constant pressure by solving equation Ψ = 0, presented below. This provides a representation of what is happening in a pore with a description of each phase in a pore.

For each pore i, the relationship between the first set of characteristics and the second set of characteristics may be given by the following system of (2N_c + 2) equations:

$\Psi =\left\{\begin{array}{c}{\Psi }_{Fi}^k={\left\{\log \left({\widehat{f}}_k^g\right)-\log \left({\widehat{f}}_k^l\right)\right\}}_i=0\\ {\Psi }_{Ci}^k={\left\{C^k-\left({\rho }^g{S}^g{y}^k+{\rho }^l{S}^l{x}^k\right)\right\}}_i=0\\ {\Psi }_{xi}={\left\{{\displaystyle \sum _{k=1}^{N_c}{x}^k}-1\right\}}_i=0\\ {\Psi }_{yi}={\left\{{\displaystyle \sum _{k=1}^{N_c}{y}^k}-1\right\}}_i=0\end{array}\right)$

where pore i in two phase state and component k varies from 1 to N_c. This system of equations may be determined in a manner similar to the system of equations from step 314 of flow 300. This provides a matrix representing the rate of change (or the derivative) of the secondary set of characteristics with respect to the first set of characteristics is obtained by

$\frac{\partial X_{sv}}{\partial X_P},$

which is used later in step 422.

In steps 416 and 418, phase conditions in each throat are determined based on connected pores using pore capillary pressures, throat invasion and snap-off events. The phase conditions may include phase occupancy in the throat representing which phases (e.g., liquid phase and/or gas phase) are present in the throat and phase conductivities in the throat. The phase occupancy in each throat is determined in step 416 may be determined in a manner similar to the phase occupancy in each throat that is determined in step 108 for flow 100.

In step 418, the phase conductivities of the liquid and gas phases in the throats are determined. The phase conductivities of the liquid and gas phases in the throats in flow 400 may be determined in a manner similar to the determination made in step 110 of flow 100.

In steps 420 to 424, the molar balance for each component, including the solvent and hydrocarbon components, in the pore network model are determined. In determining the molar balance for each component, conservation of mass within the pore network model is considered, which may then be used to examine mass transfer of the components within the pore network model. Steps 412 to 426 are similar to steps 312 to 326 in that these steps are performed for each iteration in the Newton method. The result of all iterations of steps 412 to 426 is a completed molar balance for each component for a given time, starting from a time t. As time progresses in the flow 400, the multiple iterations of steps 412 to 426 are performed at each time and thus, the molar balance of each component is determined at each time with time progressing by a time step Δt that may change based on results either between iterations within the iterative method of steps 412 to 426 or based on results from one time to a previous t-Δt (or subsequent t+Δt) time. By tracking a change in the molar balance of each component in each pore in the pore network model over time the mass transfer for each component over time can be traced.

In step 420, a representation for the molar balance for solvent and hydrocarbon components in the pore network model is established.

At each pore i = 1, N_pore, the molar balance equation for the k components (where k may be hydrocarbons, solvent, or other components in the subterranean reservoir that are to be considered) to account for mass conservation can be represented by the residual vector F as determined by:

$\begin{array}{l}{F}_i^k=V_i\frac{1}{\Delta \text{t}}\left[\left({\text{C}}_i^k\right)-{\left({\text{C}}_i^k\right)}^o\right]+{\displaystyle {\sum }_{j=1}^{N_i}{\rho }_{ij}^i{x}_{ij}^k{K}_{ij}^l\left(P_i^l-P_j^l\right)}\\ +{\displaystyle {\sum }_{j=1}^{N_i}{\rho }_i^g{x}_i^k{K}_{ij}^g\left(P_i^g-P_j^g\right)}\\ +{\displaystyle {\sum }_{j=1}^{N_i}\frac{D_{ij}^l{A}_{ij}^l}{l_{ij}}}\left[{\rho }_i^l\max \left(x_i^k-x_j^k,0\right)+{\rho }_j^l\min \left(x_i^k-x_j^k,0\right)\right]\\ +{\displaystyle {\sum }_{j=1}^{N_i}\frac{D_{ij}^g{A}_{ij}^g}{l_{ij}}}\left[{\rho }_i^g\max \left(y_i^k-y_j^k,0\right)+{\rho }_j^g\min \left(y_i^k-y_j^k,0\right)\right]=0\end{array}$

where V_i is the pore volume, N_i is the number of pore throats connected to pore body i,

$x_i^k,y_i^k$

are the component k (solvent or hydrocarbon) mole fraction in liquid and gas phases in pore i,

$S_i^l\text{and}{S}_i^g\text{are}$

the liquid phase saturation and gas phase saturation, respectively, in pore i, ρ^l, ρ^g are the liquid and gas molar densities,

$D_{ij}^{l,g}$

is the diffusion coefficient in the liquid and gasphases between pores i and j,

$P_i^l\text{and}{P}_i^g$

are the liquid and gas pressures in pore i,

$K_{ij}^l{\text{and K}}_{ij}^g$

are the conductivities of liquid and gas phases in a throat connecting pores i and j,

$A_{ij}^{\alpha }$

is the cross sectional area for phase α flow through the throat connecting pores i and j

$x_{ij}^k,y_{ij}^k$

are the component k (solvent or hydrocarbon) mole fraction in liquid and gas phases in the throat connecting pores i and j, V_i is the volume of the pore, and L_ij is the length of the throat between pores i and j.

For two-phase flow of a multicomponent flow in porous media, there is a system of (N_c + 1) N_pore equations. A residual vector F with (N_c + 1) N_pore elements is represented by

$F=\left[{\left\{F_i^{k=1,N_c}{F}_{Si}\right\}}_{i=1,N_{pore}}\right].\text{Each element}{F}_i^k\text{is}$

determined using the molar balance equation above. The last element ^Fs_i is obtained from the summation of phase saturations to 1. The system of equations (N_c + 1) N_pore is solved for the first set of characteristics X_p.

In step 322, a change in the molar balance for all components with respect to the first set of characteristics X_p is determined.

When the Newton method is used, a rate of change of the molar balance in each pore with respect to the first set of characteristics for the pore in flow 300 may be represented by a matrix, for example a Jacobian matrix used in newton iterations, given by

$Jac=\left[\frac{\delta F}{\delta X_P}\right]={\left[\frac{\delta F}{\delta X_P}\right]}_{X_{SV}}+{\left[\frac{\delta F}{\delta X_{SV}}\right]}_{X_P}\left[\frac{\delta X_{SV}}{\delta X_P}\right]$

The properties of the Jacobian matrix in flow 400 are similar to the properties of the Jacobian matrix discussed in respect of flow 300.

In step 424, a change in the first set of characteristics for a pore, represented by X_p, is determined. The change in the first set of characteristics for the pore is determined using the previously determined molar balance and the previous determined change of the molar balance in with respect to the first set of characteristics for the pore. The change in the first set of characteristics of the pore (ΔX) is determined by Jac × ΔX = -F. This change in the first set of characteristics for a pore in step 424 of flow 400 may be determined in a manner similar to the determination made in step 324 of flow 300.

In step 426, the first set of characteristics for the pore are updated based on the change in the first set of characteristics for the pore (ΔX) that was previously determined. The phase conditions determined in step 412 are used to determine whether or not the pore is in one phase or two phases. Based on this, the first set of characteristics are updated and the characteristics that form the first set of characteristics for a pore may change if the phase conditions for that pore have changed based on the discussion of the first set of characteristics with respect to step 410. For example, the first set of characteristics for the pore may be updated based on X^t+Δt,iter+1 = X^t+Δt,iter + ΔX. This updating of the first set of characteristics for a pore in step 426 of flow 400 may be determined in a manner similar to the determination made in step 118 of flow 100.

As noted above, steps 412 to 426 are iterative and are performed for a given time. Decision points at steps 428, 430 and 432 assess the results of the iterative steps and whether those results are as expected based on a selected time step Δt between times. If a time step is selected that does not produce desired results, the time step Δt may change in steps 428 or 432.

Steps 428 to 440 may be performed in a manner similar to steps 328 to 340 of flow 400. In step 428, a change in a predetermined characteristic is examined with respect to the value of that predetermined characteristics from a previous consecutive iteration, which will result in flow 400 being repeated from step 412 with a smaller time step if the value is exceeded. In step 430, the relative change of the first set of characteristics of each pore between iterations is examined, which will result in the flow 400 being repeated from step 412 with a new time if the change is below a predetermined amount. In step 432, the change in a predetermined characteristic between consecutive times is assessed, which will result in flow 400 being repeated from step 412 with a smaller time step if the value is exceeded. In step 434, the next time step is estimated for the next time to be used for flow 400. In step 436, the flow 400 is repeated from step 412 until a predetermined time is reached. In step 438, the extraction of hydrocarbons during the various times in flow 100 may be considered by assessing the change in the molar balance of the hydrocarbons in the pore network model. In step 440, the results from the pore network model are used to adjust operational parameters of wells within the subterranean reservoir.

FIGS. 6 and 7A to 7D illustrate an example implementations of the computer model using the dynamic pore network model for a multi-phase extraction process as discussed in connection with FIGS. 2, 4 and 5. FIG. 6 illustrates a plurality of pores 602 connected by throats 604 that represent a structure of a subterranean reservoir in which hydrocarbons are contained in the pores 602 and throats 604. An injection well 606 is positioned at a top of the pore network model for the injection of a solvent. A production well 608 is positioned at a bottom of the pore network model for extraction of hydrocarbons from the subterranean reservoir. Although FIG. 6 illustrates the injection well 606 and the production well 608 as separate entities, the function of both of these wells can be performed by a single well. Further, the relative position and orientation, particularly in three-dimensional space, of both the injection well 606 and the production well 608 may vary. The pores 602 in the pore network model may be illustrated as containing only a liquid phase 610, containing only a gas phase 612 or containing both the liquid phase and the gas phase (e.g., two phases 614 in the pore). FIG. 6 is a simplified example to illustrate the general structure of the pore network model and may not represent a complex interconnection of pores and throats.

FIGS. 7A to 7D illustrate an example implementation of the computer model using the dynamic pore network model for a multi-phase extraction process over time. FIGS. 7A to 7D provide a simplified representation with simplified connections between the pores. This is a simplification for illustration purposes only. FIGS. 7A to 7D have been illustrated with a presumed configuration of an injection well near the top right of the model and a production well near the top left of the mode. This configuration has been used for illustration purposes only as other configurations may also be used. FIG. 7A is illustrated with a magnified portion illustrating the pores connected by throats. This is a simplified illustration that also applies to FIGS. 7B to 7D.

FIG. 7A illustrates a time near a beginning of the computer model in which the solvent has just been injected into the subterranean reservoir and is starting to move into the pores and throats in the pore network model. FIG. 7A illustrates liquid saturation in the subterranean reservoir modeled by the pore network model in which black circles represent liquid and white circles represent gas. The hydrocarbons in the subterranean reservoir considered to be only in the liquid phase whereas the solvent may be in both the liquid phase and the gas phase.

FIG. 7B illustrates a time near a middle of the computer model illustrating liquid phase saturation into the pore network model. As can be seen in this figure, the liquid phase has started to migrate through the pores down toward the production well. FIG. 7B uses the same illustration for liquid and gas in the pores as does FIG. 7A.

The computer model of flows 100, 300 and 400 may illustrate solvent saturation through the pores, as shown in FIG. 7B or may illustrate other properties such as a solvent molar fraction as illustrated in FIG. 7C. FIG. 7C illustrates the same time as FIG. 7B but shows a solvent molar fraction in the pore network model but with different properties being illustrated . FIG. 7B illustrates solvent molar concentration with the white circles representing pores containing no solvent, the black circles representing pores containing only solvent and the grey circles representing pores containing both solvent and hydrocarbons.

FIG. 7D illustrates liquid phase saturation at a time closer to the predetermined time for the computer model of flows 100, 300 and 400 in which some hydrocarbons have been extracted from the subterranean reservoir leaving the solvent. FIG. 7D uses the same illustration for liquid and gas in the pores as does FIG. 7A. In FIG. 7D, it can be seen that some of the pores that had high liquid phase saturation in FIG. 7A now contain only the gas phase. This is indicative of hydrocarbons that have been mobilized by the solvent and extracted by the production well.

The computer model using the pore network model can be applied to any process based on in-situ solvent injection in which there are two main phases of liquid and gas, such as, for example, NSolv, based on heated solvent with condensation on bitumen interface, warm vapex or solvent+ where solvent is also heated.

In the computer model using the pore network model above, it was found that increasing solvent flow rate caused the two phase flow pattern to transition from capillary fingering to viscous fingering. This resulted in intermittent renewal of the mass transfer layer between solvent and hydrocarbons which enhanced mass transfer. Increasing solvent flow rate will change the two-phase flow pattern from capillary fingering (capillary dominated flow) to viscous fingering (viscous dominated flow). In the latter the interface between hydrocarbons and solvent is renewed as diluted hydrocarbons forms clusters in the pores moving to the outlet. The result is intermittent flow and periodic renewal of the interface causing a steep solvent concentration gradient and enhanced mass transfer.

While FIGS. 2, 4, and 5 illustrate flows in a relatively linear manner, it will be understood by those skilled in the art that various steps may be performed simultaneously with or in combination with other steps. It will also be understood that the steps may be performed in a different order than represented in the figures. It will be understood that the illustration in FIGS. 2, 4, and 5 is merely representative of an example flow but that simultaneous processing cannot be represented within the confines of FIGS. 2, 4, and 5.

Although the above description uses metric units for measurement, it will be understood that any appropriate measurement unit and any appropriate measurement system may also be used. The use of a particular measurement unit in the above description does not limit the present techniques to only the use of the above units that were used for ease of explanation of the present techniques.

The techniques described herein are not limited to heavy oils, but may also be used with any number of other subterranean reservoirs using extraction process and porous media that involves a two phase flow with mass transfer.

All of the documents mentioned above are incorporated herein by reference.

While the present techniques may be susceptible to various modifications and alternative forms, the embodiments discussed above have been shown only by way of example. However, it should again be understood that the techniques are not intended to be limited to particular embodiments disclosed herein. Indeed, the present techniques include all alternatives, modifications and equivalents falling within the scope hereof.

Claims

1. A method of operating a computer system for determining hydrocarbon extraction from a subterranean reservoir, the method comprising:

a) generating a pore network model of the subterranean reservoir representing a structure of the subterranean reservoir and conditions in the subterranean reservoir, the structure of the subterranean reservoir being represented by the pore network model as a plurality of pores connected to each other by a plurality of throats, the plurality of pores and the plurality of throats containing components including hydrocarbons;

b) establishing operational parameters for injection of a solvent into the subterranean reservoir represented by the pore network model, wherein the solvent when injected becomes one of the components in at least one from the plurality of pores and the plurality of throats;

c) based on the operational parameters for injection of the solvent, performing: i. defining, for each of the plurality of pores, a first set of characteristics based on a molar density of each component in a pore; ii. deriving, for each of the plurality of pores based on phase conditions in the pore, a second set of characteristics for the pore from the first set of characteristics for the pore; iii. determining, for each of the plurality of throats, phase conditions in a throat from the plurality of throats based on the phase conditions in ones of the plurality of pores that are connected to the throat; iv. determining a molar balance of the components within the plurality of pores and the plurality of throats in the pore network model based on the first set of characteristics and the second set of characteristics for each of the plurality of pores and the phase conditions for each of the plurality of throats; v. determining an updated first set of characteristics for each of the plurality of pores based on the molar balance of the components within the pore network model; vi. repeating (ii) to (v) over a predetermined time at consecutive times to determine hydrocarbon extraction from the subterranean reservoir over time given the operational parameters for injection of the solvent based on changes in the molar balance of the components over time representing mass transfer of the components within the pore network model; and

d) adjusting the operational parameters for injection of the solvent and repeating (c) with the adjusted operational parameters to identify operational parameters for injection of the solvent that provide a predetermined hydrocarbon extraction.

2. The method of claim 1, wherein generating the pore network model comprises:

measuring, with at least one capture device, properties of subterranean reservoir representing the structure of the subterranean reservoir and the conditions in the subterranean reservoir;

classifying the properties of the subterranean reservoir representing the structure as representing pore, matrix or throat;

generating the pore network model using the classified properties of pore for the plurality of pores and the classified properties of throat for the plurality of throats; and

establishing the properties representing the conditions in the subterranean reservoir as initial conditions for the pore network model, the initial conditions including liquid saturation in each of the plurality of pores.

3. The method of claim 1, wherein defining, for each of the plurality of pores, the first set of characteristics comprises:

establishing a constant parameter to represent conditions used in thermodynamic flash calculations;

when the constant parameter is constant volume: determining the molar density of each component in the pore network model for each of the plurality of pores as the first set of characteristics;

when the constant parameter is constant pressure: determining the molar density of each component in the pore network model and one of a gas pressure and a liquid pressure for each of the plurality of pores for the first set of characteristics;

and wherein deriving the second set of characteristics comprises: determining the phase conditions in each of the plurality of pores by determining whether the components in each of the plurality of pores is in one phase or two phases where one phase is one of a liquid phase and a gas phase and two phases includes both the liquid phase and the gas phase; when the constant parameter is constant volume: deriving, from the first set of characteristics, molar density of each component in the liquid phase, molar density of each component in the gas phase and saturation of the liquid phase for each of the plurality of pores that are in two phases as the second set of characteristics; when the constant parameter is constant pressure: deriving, from the first set of characteristics, a mole fraction of each component in the liquid phase, a mole fraction of each component in the gas phase, saturation of the liquid phase, and saturation of the gas phase for each of the plurality of pores that are in two phases as the second set of characteristics.

4. The method of claim 3, wherein defining, for each of the plurality of pores, the first set of characteristics when the constant parameter is constant pressure further comprises:

determining the mole fraction of each component based on the molar density of each component for the first set of characteristics such that the first set of characteristics is the one of the gas pressure and the liquid pressure and the mole fraction of each component.

5. The method of claim 3, wherein determining, for each of the plurality of throats, the phase conditions in the throat comprises:

determining phase occupancy for each of the plurality of throats by determining whether the components in each of the plurality of throats is in one phase or two phases where one phase is one of a liquid phase and a gas phase and two phases includes both the liquid phase and the gas phase; and

determining phase conductivity of the liquid phase and the gas phase in each of the plurality of throats.

6. The method of claim 5, wherein determining the molar balance of the components comprises:

determining a representation for the molar balance for each component in the subterranean reservoir based on a geometry of the plurality of pores, a geometry of the plurality of throats, the first set of characteristics for each of the plurality of pores, the second set of characteristics, molar density of each component in the subterranean reservoir, the phase conductivity in each of the plurality of throats, and pressure of each phase in each of the plurality of pores;

employing a numerical analysis method to determine the molar balance for each component for each of the plurality of pores based on the representation of the molar balance for each component and the first set of characteristics for each of the plurality of pores;

determining a change in the first set of characteristics for each of the plurality of pores based on the molar balance from the numerical analysis method; wherein based on the operational parameters for injection of the solvent, performing further comprises: prior to repeating (ii) to (v) over the predetermined time at consecutive times, iteratively repeating (ii) to (v) for a current time of the consecutive times until the change in the first set of characteristics is less than a predetermined change with respect to a change in the first set of characteristics from a previous iteration at the current time; and wherein determining the updated first set of characteristics for each of the plurality of pores comprises: updating the first set of characteristics for each of the plurality of pores based on the change in the first set of characteristics determined based on the molar balance from the numerical analysis method.

7. The method of claim 6, wherein employing the numerical analysis method to determine the molar balance for each component comprises:

determining a Jacobian matrix representing a rate of change of the representation of the molar balance with respect to a rate of change of the first set of characteristics; and

employing a fully implicit numerical analysis method to determine the molar balance for each component for each of the plurality of pores using the Jacobian matrix.

8. The method of claim 1, wherein determining the updated first set of characteristics for each of the plurality of pores comprises:

updating the first set of characteristics for each of the plurality of pores based on a change of the molar balance of the components for a corresponding one of the plurality of pores.

9. The method of claim 1, wherein repeating (ii) to (v) over the predetermined time comprises:

before repeating (ii) to (v), adjusting a time step between the consecutive times if a change in a predetermined characteristic between consecutive times is greater than predetermined amount; and

repeating (ii) to (v) at a new time which is the time step from a current time until the new time is the predetermined time.

10. The method of claim 1, wherein repeating (ii) to (v) over the predetermined time comprises:

determining the hydrocarbon extraction from the subterranean reservoir based on a difference in the molar balance in the pore network model at a previous consecutive time to the molar balance in the pore network model at a current time.

11. A computer system for determining hydrocarbon extraction from a subterranean reservoir, the computer system comprising:

an interface for receiving information about properties of the subterranean reservoir;

a processor configured to: a) generate a pore network model of the subterranean reservoir based on the received information about the properties of the subterranean reservoir, the pore network model representing a structure of the subterranean reservoir and conditions in the subterranean reservoir, the structure of the subterranean reservoir being represented by the pore network model as a plurality of pores connected to each other by a plurality of throats, the plurality of pores and the plurality of throats containing components including hydrocarbons; b) establish operational parameters for injection of a solvent into the subterranean reservoir represented by the pore network model, wherein the solvent when injected becomes one of the components in at least one from the plurality of pores and the plurality of throats; c) based on the operational parameters for injection of the solvent, perform: i. define, for each of the plurality of pores, a first set of characteristics based on a molar density of each component in a pore; ii. derive, for each of the plurality of pores based on phase conditions in the pore, a second set of characteristics for the pore from the first set of characteristics for the pore; iii. determine, for each of the plurality of throats, phase conditions in a throat from the plurality of throats based on the phase conditions in ones of the plurality of pores that are connected to the throat; iv. determine a molar balance of the components within the plurality of pores and the plurality of throats in the pore network model based on the first set of characteristics and the second set of characteristics for each of the plurality of pores and the phase conditions for each of the plurality of throats; v. determine an updated first set of characteristics for each of the plurality of pores based on the molar balance of the components within the pore network model; vi. repeat (ii) to (v) over a predetermined time at consecutive times to determine hydrocarbon extraction from the subterranean reservoir over time given the operational parameters for injection of the solvent based on changes in the molar balance of the components over time representing mass transfer of the components within the pore network model; and d) adjust the operational parameters for injection of the solvent and repeat (c) with the adjusted operational parameters to identify operational parameters for injection of the solvent that provide a predetermined hydrocarbon extraction.

12. The system of claim 11, further comprising:

a capture device for measuring the properties of the subterranean reservoir, the properties representing the structure of the subterranean reservoir and the conditions in the subterranean reservoir, the capture device being in communication with the interface to provide the information about the properties of the subterranean reservoir; wherein the processor is further configured when generating the pore network model to: classify the properties of the subterranean reservoir representing the structure as representing pore, matrix or throat; generate the pore network model using the classified properties of pore for the plurality of pores and the classified properties of throat for the plurality of throats; and establish the properties representing the conditions in the subterranean reservoir as initial conditions for the pore network model, the initial conditions including liquid saturation in each of the plurality of pores.

13. The system of claim 11, wherein the processor is further configured when defining, for each of the plurality of pores, the first set of characteristics to:

establish a constant parameter to represent conditions used in thermodynamic flash calculations;

when the constant parameter is constant volume: determine the molar density of each component in the pore network model for each of the plurality of pores as the first set of characteristics;

when the constant parameter is constant pressure: determine the molar density of each component in the pore network model; and determine a mole fraction of each component based on the molar density of each component and one of a gas pressure and a liquid pressure for each of the plurality of pores as the first set of characteristics; and wherein the processor is further configured when deriving the second set of characteristics to: determine the phase conditions in each of the plurality of pores by determining whether the components in each of the plurality of pores is in one phase or two phases where one phase is one of a liquid phase and a gas phase and two phases includes both the liquid phase and the gas phase; when the constant parameter is constant volume: derive, from the first set of characteristics, molar density of each component in the liquid phase, molar density of each component in the gas phase and saturation of the liquid phase for each of the plurality of pores that are in two phases as the second set of characteristics; when the constant parameter is constant pressure: derive, from the first set of characteristics, a mole fraction of each component in the liquid phase, a mole fraction of each component in the gas phase, saturation of the liquid phase, saturation of the gas phase for each of the plurality of pores that are in two phases as the second set of characteristics.

14. The system of claim 13, wherein the processor is further configured when determining, for each of the plurality of throats, the phase conditions in the throat to:

determine phase occupancy for each of the plurality of throats by determining whether the components in each of the plurality of throats is in one phase or two phases where one phase is one of a liquid phase and a gas phase and two phases includes both the liquid phase and the gas phase; and

determine phase conductivity of the liquid phase and the gas phase in each of the plurality of throats.

15. The system of claim 14, wherein the processor is further configured when determining the molar balance of the components to:

determine a representation for the molar balance for each component in the subterranean reservoir based on a geometry of the plurality of pores, a geometry of the plurality of throats, the first set of characteristics for each of the plurality of pores, the second set of characteristics, molar density of each component in the subterranean reservoir, the phase conductivity in each of the plurality of throats, and pressure of each phase in each of the plurality of pores;

employ a numerical analysis method to determine the molar balance for each component for each of the plurality of pores based on the representation of the molar balance for each component and the first set of characteristics for each of the plurality of pores;

determine a change in the first set of characteristics for each of the plurality of pores based on the molar balance from the numerical analysis method; and wherein the processor is further configured when, based on the operational parameters for injection of the solvent, performing to: prior to repeating (ii) to (v) over the predetermined time at consecutive times, iteratively repeat (ii) to (v) for a current time of the consecutive times until the change in the first set of characteristics is less than a predetermined change with respect to a change in the first set of characteristics from a previous iteration at the current time; and wherein the processor is further configured when determining the updated first set of characteristics for each of the plurality of pores to: update the first set of characteristics for each of the plurality of pores based on the change in the first set of characteristics determined based on the molar balance from the numerical analysis method.

16. The system of claim 15, wherein the processor is further configured when employing the numerical analysis method to determine the molar balance for each component to:

determine a Jacobian matrix representing a rate of change of the representation of the molar balance with respect to a rate of change of the first set of characteristics; and

employing a fully implicit numerical analysis method, determine the molar balance for each component for each of the plurality of pores using the Jacobian matrix.

17. The system of claim 11, wherein the processor is further configured when determining the updated first set of characteristics for each of the plurality of pores to:

update the first set of characteristics for each of the plurality of pores based on the change of the molar balance of the components for a corresponding one of the plurality of pores.

18. The system of claim 11, wherein the processor is further configured when repeating (ii) to (v) over the predetermined time to:

before repeating (ii) to (v), adjust a time step between the consecutive times if a change in a predetermined characteristic between consecutive times is greater than predetermined amount; and

repeat (ii) to (v) at a new time which is the time step from a current time until the new time is the predetermined time.

19. The system of claim 11, wherein the processor is further configured when repeating (ii) to (v) over the predetermined time to:

determine the hydrocarbon extraction from the subterranean reservoir based on a difference in the molar balance in the pore network model at a previous consecutive time to the molar balance in the pore network model at a current time.

20. A computer readable medium having stored thereon statements and instructions that when executed by a computer system operate the computer system to determine hydrocarbon extraction from a subterranean reservoir, the statements and instructions when executed by the computer system cause the computer system to:

a) generate a pore network model of the subterranean reservoir representing a structure of the subterranean reservoir and conditions in the subterranean reservoir, the structure of the subterranean reservoir being represented by the pore network model as a plurality of pores connected to each other by a plurality of throats, the plurality of pores and the plurality of throats containing components including hydrocarbons;

b) establish operational parameters for injection of a solvent into the subterranean reservoir represented by the pore network model, wherein the solvent when injected becomes one of the components in at least one from the plurality of pores and the plurality of throats;

c) based on the operational parameters for injection of the solvent, perform: (i) define, for each of the plurality of pores, a first set of characteristics based on a molar density of each component in a pore; (ii) derive, for each of the plurality of pores based on phase conditions in the pore, a second set of characteristics for the pore from the first set of characteristics for the pore; (iii) determine, for each of the plurality of throats, phase conditions in a throat from the plurality of throats based on the phase conditions in ones of the plurality of pores that are connected to the throat; (iv) determine a molar balance of the components within the plurality of pores and the plurality of throats in the pore network model based on the first set of characteristics and the second set of characteristics for each of the plurality of pores and the phase conditions for each of the plurality of throats; (v) determine an updated first set of characteristics for each of the plurality of pores based on the molar balance of the components within the pore network model; (vi) repeat (ii) to (v) over a predetermined time at consecutive times to determine hydrocarbon extraction from the subterranean reservoir over time given the operational parameters for injection of the solvent based on changes in the molar balance of the components over time representing mass transfer of the components within the pore network model; and

d) adjust the operational parameters for injection of the solvent and repeat (c) with the adjusted operational parameters to identify operational parameters for injection of the solvent that provide a predetermined hydrocarbon extraction.