# Predicting Fluid Properties in the MUFITS Reservoir Simulator with User-Supplied Modules

Authors: Andrey Afanasyev; Ivan Utkin
Journal: Geofluids
 (2021)

---

## Abstract

We present a recent development of the MUFITS reservoir simulator aiming at modelling the transport of fluids whose properties and phase equilibria are calculated in a user-supplied external shared library. Both the explicit correlations and tabulated data for the fluid parameters can be implemented in the library that we name the EoS-module (Equation of State-module). An iterative approach—which, for example, is based on the phase equilibria calculation through the Gibbs energy minimisation (GEM) method, can also be used in the EoS-module. A considerable effort has been undertaken to minimise the number of program procedures exported by the shared library. This should facilitate and ease the usage of the developed software extension by the scientific community. Furthermore, we supplement the article with the source code of two simple EoS-modules that can serve as templates in other modelling and software development efforts. The EoS-modules are also useful for coupling MUFITS with other elaborate software for fluid property prediction. To demonstrate such a possibility, we supplement the article with the source code of a more complicated EoS-module that couples MUFITS with the geochemical code GEMS3K. This module is used in a simple 1-D benchmark study showing the capabilities of MUFITS for modelling reactive transport in porous media.


---

## Body

## 1. Introduction

### 1.1. Reservoir Simulations

The numerical modelling of multicomponent multiphase flows in porous media is in demand in many areas of subsurface exploration and natural process prediction. These areas involve petroleum reservoir exploitation, subsurface CO2 and natural gas storage, geothermal energy extraction, the prediction of interactions between hydrothermal and volcanic systems, and other applications [1, 2]. The related transport in porous media is often complicated by multiphase equilibria of reservoir fluids, nonisothermal processes, and reactions between fluid and rock [3]. Many other complications can be observed if the flows occur in fractured reservoirs, the transport in a porous medium is coupled with multiphase flows in wellbores and surface facilities, or there is a complicated, i.e., nonequilibrium, behaviour at the pore scale [4].Given that the aforementioned complications often lead to nonlinear equations solved in a 3-D domain, the flows can be predicted only by means of numerical modelling, i.e., reservoir simulations [4, 5]. A computer program accounting for all relevant physical phenomena and exploitation processes in a convenient integrated manner is called a reservoir simulator. Certainly, reservoir simulators differ in several respects, i.e., target applications or available modelling options. There are academic, not-for-profit simulators often applied in the investigation of processes in nature, e.g., TOUGH2 [6], TOUGHREACT [7], and MODFLOW [8], among others. Usually, such simulators have a limited functionality, allowing users without sufficient backgrounds in reservoir simulations to do easy and quick assessments of some simple transport phenomena. Normally, the commercial software includes a much wider set of options to allow for modelling with engineering accuracy [4, 5].Some simulators take an intermediate place between open-source and commercial software. Such reservoir simulation codes, e.g., MUFITS [9, 10] and AD-GPRS [11], are used to develop and test new methods for numerical modelling because their object-oriented architectures are closer to those of commercial software. If these methods prove their robustness, then they are applied in commercial software. Thus, such codes serve as platforms for the improvement of more elaborate simulators.


### 1.2. Governing Equations for Modelling Transport in Porous Media

If we neglect any specific complications, most reservoir simulators are aimed at solving the following governing equations [4–6]:
(1)∂∂tϕ∑i=1npρicijsi+∇·∑i=1npρicijui+ψj=0,(2)∂∂tϕ∑i=1npρieisi+1−ϕρrer+∇·∑i=1npρihiui+ψe=0,(3)ui=−KKriμi∇Pi−ρig,(4)ΩP,T,ct,whereΩ=np,ρi,ei,hi,μi,si,cij,ct=ct1,ct2,⋯,ctnc,(5)Kri=Kris,Pi−Pk=Pc,iks,wheres=s1,⋯,snp,(6)∑i=1npsi=1,∑j=1nccij=1,∑j=1ncctj=1.Here,nc is the number of fluid components; np is the number of fluid phases; ϕ is the porosity; ρ is the density; cij is the jth component mass fraction in the ith phase; s is the phase saturation (i.e., volume fraction of pore space occupied by the phase); u is the Darcy velocity; ψj and ψe are the fluxes of the jth component and energy caused by mechanical dispersion, molecular diffusion, or heat conduction; e is the specific energy; h is the specific enthalpy; K is the absolute permeability of the medium; Kr is the relative permeability; μ is the dynamic viscosity; P is the pressure; g is the gravity acceleration; ctj is the bulk mass fraction of the jth component; T is the temperature; Pc is the capillary pressure; the subscripts i, j, and r denote the parameters of the ith phase, the jth component, and rock, respectively.Equation (1) is the mass conservation equation for each fluid component, Equation (2) is the energy conservation equation, and Equation (3) is Darcy’s law. Equation (4) is a general formulation of relations used for the fluid property prediction. It means that the phase equilibrium depends on the average fluid pressure P, temperature T, and bulk composition ct. Equation (5) defines the relative permeability and capillary pressure as saturation functions. Equations in (6) are additional relations. For simplicity, we assume that ϕ, K, and ρr are constants and that erT is a given function of temperature. If the isothermal flow is simulated, then Equation (2) should be excluded from the system of governing equations, and the condition T=const should be assumed.In the case of mechanical dispersion, the fluxesψj take the following form [12]:
(7)ψj=−ϕ∑i=1npsiDi∇cij,(8)Di=1ϕsiλL,iux,i+λT,iuy,i+uz,i000λL,iuy,i+λT,iux,i+uz,i000λL,iuz,i+λT,iux,i+uy,i,where D is the tensor of mechanical dispersion; λL and λT are the longitudinal and transverse mean free path of a small volume of the ith phase caused by mixing in the porous medium; and ux, uy, and uz are the components of the Darcy velocity.The governing equations can be distinctively split into two sets of equations. The first set describing the fluid transport includes Equations (1)–(3), (7), and (8). The second set of Equations (4) and (5) is responsible for predicting the phase equilibria and transport properties of the fluid. Often only the equations of the second set are different, and the equations of the first set are the same in different applications (Figure 1).Figure 1
Architecture of MUFITS (modified after [9]).The described splitting of the governing equations is particularly useful in designing reservoir simulators. It results in a natural architecture of the software, involving two large program modules (Figure1). One module that we further refer to as the “kernel” involves all procedures not related to predicting phase equilibria and all fluid properties in general. Such a module includes finite-difference schemes, linear and nonlinear solvers, and other procedures related to the solution of Equations (1)–(8). The second module that we further refer to is the “EoS-module” that contains only the procedures for calculating fluid properties. Obviously, the kernel is larger and more elaborate than an EoS-module. However, an EoS-module can also be quite elaborate if it is based on the iterative prediction of phase equilibria. For example, such complications arise in compositional reservoir simulations [13] or reactive transport modelling [14, 15].The modular architecture is very convenient for extending the application area of reservoir simulators. The kernel can be linked with different EoS-modules. Thereby, the simulator can be used for modelling the transport of different types of fluids. If the development of an EoS-module is not very time-consuming, then such an extension is easy to achieve. Moreover, such an extension has all the modelling options already implemented in the kernel, e.g., the mechanical dispersion (see Equations (7) and (8)), coupled reservoir-wellbore flows, and other sophisticated options. The discussed modular architecture is explicitly implemented in TOUGH2 [6] and MUFITS [9] codes.


### 1.3. Goal of This Work

This article is aimed at presenting the recent extension of the MUFITS simulator that makes possible the development of user-supplied EoS-modules (Figure1). MUFITS is an academic noncommercial software that has continued to develop over the past decade [9, 10, 12]. The software’s kernel allows for parallel simulations, the coupled modelling of flows in porous media and wellbores, the modelling of the mechanical dispersion, and other options. For example, MUFITS has recently been extended for the automated prediction of and history matching to Earth’s surface displacement and gravity changes associated with subsurface flows [16].The simulator has several built-in EoS-modules (Figure1). For example, it includes EoS-modules for the black-oil modelling of petroleum reservoirs and the subsurface storage of natural gas or CO2 in saline aquifers. Although the software includes the standard modules for such applications, permanent demand exists in the development of new specific EoS-modules intended to be applied in nonstandard cases and the narrow area of research and engineering applications. To address such demand, we have implemented the possibility of linking the MUFITS kernel with user-supplied EoS-modules. Now, an EoS-module for a specific fluid can be developed by the user and linked with the software as an external shared library. Thus, from now on, the MUFITS extension for modelling the flows of other fluids can be done solely by the users. To ease and facilitate such extensions, we have reduced the number of procedures exported from the shared library to the kernel to just four.Organising the development of an EoS-module as shared library with a convenient calling interface constitutes the advancement and innovation of the presented software’s extension. Usually, the extension of reservoir simulator for the flows of a specific fluid requires significant programming efforts that often can be carried out based on only in-depth knowledge of software including its kernel. A rather complicated interface between the kernel and EoS-module requires that EoS-module is embedded in software as static library. For example, such static EoS-modules exist in MUFITS [9, 10] (Figure 1), TOUGH2 [6], and other codes. Usually, extending a reservoir simulator to reactive transport also assumes static compilation [7, 8], which limits the flexibility of software development and modification. To some degree, the calling interface and the dynamic linking proposed in this study solve the problem by allowing users to develop their own EoS-modules without in-depth knowledge on the whole software’s architecture.We supplement the article with a few simple open-source examples of the dynamically linked EoS-modules and benchmark them against 1-D analytical solutions. In the future, these examples can serve as templates for other more complicated modules. The discussed development is also useful for coupling MUFITS with other elaborate software for fluid properties prediction. Here, we present the result of MUFITS coupling with the geochemical software GEMS3K [17, 18], which extends the simulator for reactive transport modelling [1, 2, 14]. The software is based on the Gibbs energy minimisation (GEM) approach, which allows for finding unknown phase equilibria containing a considerable number of fluid and solid phases. We supplement the article with the source code of the EoS-module that dynamically links the MUFITS kernel with GEMS3K, compiled as a shared library. The coupling is validated against the 1-D benchmark study of reactive transport.This article is organised as follows. In Section2, we describe the architecture of a user-supplied EoS-module, provide an overview of the primary data flows, and discuss the application program interface (API) for all the procedures exported by the shared library. In Section 3, we present two simple EoS-modules, one for a mixture of ideal gases and another for a two-phase fluid. Section 4 is devoted to MUFITS coupling with GEMS3K. We end the article with Conclusions, where we briefly outline possible further extensions of development.


## 1.1. Reservoir Simulations

The numerical modelling of multicomponent multiphase flows in porous media is in demand in many areas of subsurface exploration and natural process prediction. These areas involve petroleum reservoir exploitation, subsurface CO2 and natural gas storage, geothermal energy extraction, the prediction of interactions between hydrothermal and volcanic systems, and other applications [1, 2]. The related transport in porous media is often complicated by multiphase equilibria of reservoir fluids, nonisothermal processes, and reactions between fluid and rock [3]. Many other complications can be observed if the flows occur in fractured reservoirs, the transport in a porous medium is coupled with multiphase flows in wellbores and surface facilities, or there is a complicated, i.e., nonequilibrium, behaviour at the pore scale [4].Given that the aforementioned complications often lead to nonlinear equations solved in a 3-D domain, the flows can be predicted only by means of numerical modelling, i.e., reservoir simulations [4, 5]. A computer program accounting for all relevant physical phenomena and exploitation processes in a convenient integrated manner is called a reservoir simulator. Certainly, reservoir simulators differ in several respects, i.e., target applications or available modelling options. There are academic, not-for-profit simulators often applied in the investigation of processes in nature, e.g., TOUGH2 [6], TOUGHREACT [7], and MODFLOW [8], among others. Usually, such simulators have a limited functionality, allowing users without sufficient backgrounds in reservoir simulations to do easy and quick assessments of some simple transport phenomena. Normally, the commercial software includes a much wider set of options to allow for modelling with engineering accuracy [4, 5].Some simulators take an intermediate place between open-source and commercial software. Such reservoir simulation codes, e.g., MUFITS [9, 10] and AD-GPRS [11], are used to develop and test new methods for numerical modelling because their object-oriented architectures are closer to those of commercial software. If these methods prove their robustness, then they are applied in commercial software. Thus, such codes serve as platforms for the improvement of more elaborate simulators.


## 1.2. Governing Equations for Modelling Transport in Porous Media

If we neglect any specific complications, most reservoir simulators are aimed at solving the following governing equations [4–6]:
(1)∂∂tϕ∑i=1npρicijsi+∇·∑i=1npρicijui+ψj=0,(2)∂∂tϕ∑i=1npρieisi+1−ϕρrer+∇·∑i=1npρihiui+ψe=0,(3)ui=−KKriμi∇Pi−ρig,(4)ΩP,T,ct,whereΩ=np,ρi,ei,hi,μi,si,cij,ct=ct1,ct2,⋯,ctnc,(5)Kri=Kris,Pi−Pk=Pc,iks,wheres=s1,⋯,snp,(6)∑i=1npsi=1,∑j=1nccij=1,∑j=1ncctj=1.Here,nc is the number of fluid components; np is the number of fluid phases; ϕ is the porosity; ρ is the density; cij is the jth component mass fraction in the ith phase; s is the phase saturation (i.e., volume fraction of pore space occupied by the phase); u is the Darcy velocity; ψj and ψe are the fluxes of the jth component and energy caused by mechanical dispersion, molecular diffusion, or heat conduction; e is the specific energy; h is the specific enthalpy; K is the absolute permeability of the medium; Kr is the relative permeability; μ is the dynamic viscosity; P is the pressure; g is the gravity acceleration; ctj is the bulk mass fraction of the jth component; T is the temperature; Pc is the capillary pressure; the subscripts i, j, and r denote the parameters of the ith phase, the jth component, and rock, respectively.Equation (1) is the mass conservation equation for each fluid component, Equation (2) is the energy conservation equation, and Equation (3) is Darcy’s law. Equation (4) is a general formulation of relations used for the fluid property prediction. It means that the phase equilibrium depends on the average fluid pressure P, temperature T, and bulk composition ct. Equation (5) defines the relative permeability and capillary pressure as saturation functions. Equations in (6) are additional relations. For simplicity, we assume that ϕ, K, and ρr are constants and that erT is a given function of temperature. If the isothermal flow is simulated, then Equation (2) should be excluded from the system of governing equations, and the condition T=const should be assumed.In the case of mechanical dispersion, the fluxesψj take the following form [12]:
(7)ψj=−ϕ∑i=1npsiDi∇cij,(8)Di=1ϕsiλL,iux,i+λT,iuy,i+uz,i000λL,iuy,i+λT,iux,i+uz,i000λL,iuz,i+λT,iux,i+uy,i,where D is the tensor of mechanical dispersion; λL and λT are the longitudinal and transverse mean free path of a small volume of the ith phase caused by mixing in the porous medium; and ux, uy, and uz are the components of the Darcy velocity.The governing equations can be distinctively split into two sets of equations. The first set describing the fluid transport includes Equations (1)–(3), (7), and (8). The second set of Equations (4) and (5) is responsible for predicting the phase equilibria and transport properties of the fluid. Often only the equations of the second set are different, and the equations of the first set are the same in different applications (Figure 1).Figure 1
Architecture of MUFITS (modified after [9]).The described splitting of the governing equations is particularly useful in designing reservoir simulators. It results in a natural architecture of the software, involving two large program modules (Figure1). One module that we further refer to as the “kernel” involves all procedures not related to predicting phase equilibria and all fluid properties in general. Such a module includes finite-difference schemes, linear and nonlinear solvers, and other procedures related to the solution of Equations (1)–(8). The second module that we further refer to is the “EoS-module” that contains only the procedures for calculating fluid properties. Obviously, the kernel is larger and more elaborate than an EoS-module. However, an EoS-module can also be quite elaborate if it is based on the iterative prediction of phase equilibria. For example, such complications arise in compositional reservoir simulations [13] or reactive transport modelling [14, 15].The modular architecture is very convenient for extending the application area of reservoir simulators. The kernel can be linked with different EoS-modules. Thereby, the simulator can be used for modelling the transport of different types of fluids. If the development of an EoS-module is not very time-consuming, then such an extension is easy to achieve. Moreover, such an extension has all the modelling options already implemented in the kernel, e.g., the mechanical dispersion (see Equations (7) and (8)), coupled reservoir-wellbore flows, and other sophisticated options. The discussed modular architecture is explicitly implemented in TOUGH2 [6] and MUFITS [9] codes.


## 1.3. Goal of This Work

This article is aimed at presenting the recent extension of the MUFITS simulator that makes possible the development of user-supplied EoS-modules (Figure1). MUFITS is an academic noncommercial software that has continued to develop over the past decade [9, 10, 12]. The software’s kernel allows for parallel simulations, the coupled modelling of flows in porous media and wellbores, the modelling of the mechanical dispersion, and other options. For example, MUFITS has recently been extended for the automated prediction of and history matching to Earth’s surface displacement and gravity changes associated with subsurface flows [16].The simulator has several built-in EoS-modules (Figure1). For example, it includes EoS-modules for the black-oil modelling of petroleum reservoirs and the subsurface storage of natural gas or CO2 in saline aquifers. Although the software includes the standard modules for such applications, permanent demand exists in the development of new specific EoS-modules intended to be applied in nonstandard cases and the narrow area of research and engineering applications. To address such demand, we have implemented the possibility of linking the MUFITS kernel with user-supplied EoS-modules. Now, an EoS-module for a specific fluid can be developed by the user and linked with the software as an external shared library. Thus, from now on, the MUFITS extension for modelling the flows of other fluids can be done solely by the users. To ease and facilitate such extensions, we have reduced the number of procedures exported from the shared library to the kernel to just four.Organising the development of an EoS-module as shared library with a convenient calling interface constitutes the advancement and innovation of the presented software’s extension. Usually, the extension of reservoir simulator for the flows of a specific fluid requires significant programming efforts that often can be carried out based on only in-depth knowledge of software including its kernel. A rather complicated interface between the kernel and EoS-module requires that EoS-module is embedded in software as static library. For example, such static EoS-modules exist in MUFITS [9, 10] (Figure 1), TOUGH2 [6], and other codes. Usually, extending a reservoir simulator to reactive transport also assumes static compilation [7, 8], which limits the flexibility of software development and modification. To some degree, the calling interface and the dynamic linking proposed in this study solve the problem by allowing users to develop their own EoS-modules without in-depth knowledge on the whole software’s architecture.We supplement the article with a few simple open-source examples of the dynamically linked EoS-modules and benchmark them against 1-D analytical solutions. In the future, these examples can serve as templates for other more complicated modules. The discussed development is also useful for coupling MUFITS with other elaborate software for fluid properties prediction. Here, we present the result of MUFITS coupling with the geochemical software GEMS3K [17, 18], which extends the simulator for reactive transport modelling [1, 2, 14]. The software is based on the Gibbs energy minimisation (GEM) approach, which allows for finding unknown phase equilibria containing a considerable number of fluid and solid phases. We supplement the article with the source code of the EoS-module that dynamically links the MUFITS kernel with GEMS3K, compiled as a shared library. The coupling is validated against the 1-D benchmark study of reactive transport.This article is organised as follows. In Section2, we describe the architecture of a user-supplied EoS-module, provide an overview of the primary data flows, and discuss the application program interface (API) for all the procedures exported by the shared library. In Section 3, we present two simple EoS-modules, one for a mixture of ideal gases and another for a two-phase fluid. Section 4 is devoted to MUFITS coupling with GEMS3K. We end the article with Conclusions, where we briefly outline possible further extensions of development.


## 2. Architecture of an EoS-Module

### 2.1. Overview

#### 2.1.1. Data Flows

A sketch of data flows associated with any user-supplied EoS-module is shown in Figure2. Further, we refer to such a module as USEREOS. Figure 2(a) outlines the four mandatory program procedures of the module. Three of them are the preprocessing routines called to choose the modelling options and specify the arrays’ sizes and the general (i.e., global) parameters of a fluid model. The fourth routine, named PhaseEquilibrium, supplies the kernel with fluid parameters for a given P, T, and ct. This routine is called many times for any cell (i.e., grid block) of the reservoir model during time stepping. The coefficients of the fluid model are loaded directly into USEREOS by reading the configuration file.Figure 2
Sketch of data flows associated with USEREOS. The green and red dashed boxes in the right panel show the global and local parameters, respectively.Figure2(b) summarises the parameters and arrays in the kernel that are changed by or loaded into USEREOS. The data are gathered in four main blocks. The first block contains all the constants controlling memory usage. The second block, “Fluid model,” contains the parameters characterising the fluid and its phase and chemical equilibria. Some of these parameters that are shown within the green dashed box are global, i.e., they are identical for all the grid cells. Other parameters of the fluid that are within the red dashed box are local, i.e. they may be different in different cells. The third data block is for an auxiliary array that can be used as additional memory provided in the kernel for every cell. Again, this data block has global and local parameters. The fourth data block contains modelling options.


#### 2.1.2. Constants Controlling Memory Usage

There are just three constants controlling memory allocation in the kernel.(i)
nc is the number of fluid components (see also Equations (1)–(8))(ii)
npmax is the maximum number of fluid phases. The number of fluid phases np depends on P, T, and ct. Thus, it can be different in different grid cells. However, for any P, T, and ct, the number of phases exported by USEREOS must not be larger than npmax, i.e., np≤npmax(iii)
na is the size of the auxiliary array. It specifies the number of additional 8-byte real numbers allocated in the kernel for every grid cell


#### 2.1.3. Parameters of the Fluid Model

Every fluid componentj=1,⋯,nc is characterised by two global parameters:
(i)
C‐Name(j) is a character name of the component (e.g., ‘H2O’ or ‘CO2’)(ii)
Mj is the molar weight of the component, which is used in the kernel for calculating molar quantitiesEvery fluid phasei=1,…,np is characterised by one global and 6+nc local parameters. The global parameter P‐Namei is a character phase name (e.g., “GAS” or “LIQ”). The local parameters correspond to the vector
(9)Φi=ρi,hi,μi,si,Kri,ΔPi,ci1,⋯,cinc,where
(10)ΔPi=Pi−P,is the difference between the phase fluid pressure Pi and the average pressure P. The variable ΔPi is used for modelling the influence of capillary pressure Pc. For example, if the fluid can split into just two phases, liquid (i=1) and gas (i=2), then one can choose P=P1, ΔP1=0, and ΔP2=P2–P1=Pc,21. For every Φi, the specific energy of every phase is calculated in the kernel as
(11)ei=hi−Piρi.


#### 2.1.4. Auxiliary Array

The auxiliary array of sizena is allocated in the kernel for every grid cell. This array can be used as additional memory provided for USEREOS. Some parameters can be stored in this array for later usage in the following time steps. For example, the auxiliary array can be used for modelling the reaction kinetics, the hysteresis phenomena, or calculating additional (secondary) parameters of the fluid not belonging to Φi (see Equation (9)). The latter is used when coupling MUFITS with GEMS (see Section 4).The elements of the auxiliary array can also be saved in the summary files in every report time period for later usage in the postprocessing of the simulation results. To make such reporting possible, every elementk=1,⋯,na of the array is characterised by two global parameters:
(i)
A‐Namek is a character name of the element that must begin with the symbol “#.” One can refer to the element by its name in the input data to MUFITS(ii)
Unitsk is a string specifying the units of the physical quantity stored in the element. These units are reported into the summary files. For dimensionless quantities, Unitsk should be an empty string or “NODIM.”Every elementk of the array is also characterised by one numeric quantity ak of type “double” (8-byte real number). The variable ak is local, i.e., it can be changed in USEREOS in every grid cell and at every time step.


#### 2.1.5. The Option Array

The option array is an integer array of length 8,optl, where l=1,⋯,8. If optl=1, then the corresponding option is enabled (on). If optl=0, then the option is disabled (off). By default, all the options are disabled. At present, only opt1 and opt2 can be used, whereas the other elements of the option array, optl,l=3,⋯,8, are reserved for future developments.The optionopt1 controls the calculation of the relative permeability Kri. If opt1=0, then the relative permeability is calculated in the kernel using standard input data to the simulator (e.g., specified with the SATTAB keyword). In this case, there is no need to export Kri from USEREOS, i.e., the 5th element of Φi can be an undefined quantity. If opt1=1, then the default treatment of relative permeability is overridden by the value Kri exported from USEREOS.The optionopt2 controls the calculation of the capillary pressure Pc or, to be strict, the relative phase pressure ΔPi (see Equation (10)). If opt2=0, then the relative pressure is calculated in the kernel using standard input data to the simulator (e.g., provided with the SATTAB keyword). In this case, there is no need to export ΔPi from USEREOS, i.e., the 6th element of Φi can be an undefined quantity. If opt2=1, then the default treatment of capillary pressure is overridden by the value ΔPi exported from USEREOS.


#### 2.1.6. Configuration File

The configuration file is an input data file to USEREOS. The path to this file is imported from the kernel. It is specified as the second argument of the USEREOS keyword that activates the presented modelling option (see Section2.3).The format of the configuration file can be specific for every particular shared library, because the kernel does not read this file. Every new USEREOS module can be supplemented with its own format of the file.It is implied that the configuration file contains coefficients of the fluid model that are not explicitly defined as static variables in USEREOS. The configuration file can contain paths to other external files loaded into USEREOS.


### 2.2. Subroutines Exported from USEREOS

#### 2.2.1. General Description

In this section, we give a brief overview of the main subroutines exported from USEREOS. The four subroutines placed in the order of their invoking by the kernel are as follows:(i)
ReadConfigurationFile. This routine is called by the kernel at the initialisation stage, immediately after the USEREOS keyword is encountered in the input data file to MUFITS (Section 2.3). The path to the configuration file is the second argument of the keyword that is imported by USEREOS as an argument of the subroutine. It is implied that the shared library reads the file and then exports the title of the module to the kernel. The title may contain a brief description of the module and its version. The title is reported to the LOG-file of MUFITS.(ii)
GetDimensions. This subroutine is called at the end of initialisation stage. It is aimed at importing the sizes of arrays that need to be allocated into the kernel. The subroutine exports the number of fluid components nc, the maximum number of fluid phases (npmax), and the auxiliary array size (na).(iii)
GetGlobalParameters. The procedure is called immediately after the previous one. It is aimed at loading the global parameters of the fluid model into the kernel. First, the 3-byte character name (C‐Name(j)) and the molar weight (Mj) of every component j=1,…,nc must be specified. Second, the 3-byte character name (P‐Namei) for each fluid phase i=1,…,npmax is exported by USEREOS. Third, the 8-byte character name (A-Namek) and units (Unitsk) for every element k=1,…,na of the auxiliary array must be specified. Any A-Namek must begin with the “#” character. Fourth, the options array, optl, is exported.(iv)
PhaseEquilibrium. Normally, this routine is called multiple times for every grid cell during linearisation of the governing Equations (1)–(8) and time stepping. It is aimed at reporting the phase equilibrium for a given P, T, and ct. It reports the number of phases np≤npmax, the vector Φi, and the integer variable Idi for every phase i=1,…,np. The quantity Idi (the “phase identifier”) is aimed at associating the ith phase with the phase names (i.e., phase types) specified by GetGlobalParameters. The value of Idi must be equal to the serial number of the P-Namek, k=1,…,npmax specified by GetGlobalParameters.A very important argument ofPhaseEquilibrium is Mode. If Mode=0, then the subroutine is supplied by the kernel with the initial guess for the phase equilibrium. The variables np, Idi, and Φi, i=1,…,np, are used to specify the initial guess. In the case of Mode=0, PhaseEquilibrium must calculate the equilibrium using the provided initial guess. The procedure must not change np and Idi. It can only change Φi for a given P, T, and ctj. Since the number of phases np is not changed, the reported saturations si and concentrations cij and ctj are allowed to take negative or larger-than-one values [19]. If Mode=1, then the initial guess is not specified. In this case, the subroutine must calculate the phase equilibrium—i.e., np, Idi, and Φi—from scratch. If Mode=1, then the reported si and cij must satisfy the inequalities 0≤si≤1 and 0≤cij≤1, respectively.


#### 2.2.2. Calling Interface for the USEREOS Subroutines

The USEREOS API is provided in the Fortran programming language. However, USEREOS can be written in any language that supports compilation to a shared library. To export the subroutines described in Section2.2.1, it may be necessary to change these subroutines’ signatures to accommodate differences between programming languages. Since GEMS3K is written in C++, the USEREOS module for coupling GEMS3K and MUFITS was written in C++ as well (Section 4).The Fortran API for the four primary subroutines exported by USEREOS is given in Tables1–4. In the tables, every argument of each subroutine is supplemented with a brief description, its datatype, and the intent. The intent values “in” or “out” mean that the corresponding variable is imported to or exported from USEREOS, respectively. The value “in/out” means that variable is used as both an input and output parameter. The datatypes are shown in a standard format used in Fortran. For example, the datatypes “integer(4)” and “real(8)” correspond to 4-byte integer and 8-byte real numbers, respectively. All physical quantities must be given in SI units. For example, the units of ρi, hi, μi, and ΔPi in the vector Φi must be kg/m3, J/kg, Pa∙s, and Pa, respectively.Table 1
API forReadConfigurationFile.
Interface: call ReadConfigurationFile (filename ,ierr, title)VariableIntentDatatypeDescriptionfilenameincharacterConfiguration file name (array of size 256)ierroutinteger(4)Error specifier (ifierr=0, then no error occurred during subroutine execution; if ierr=−1, then the configuration file is not found)titleoutcharacterEoS-module title, i.e., a brief description (array of size 80)Table 2
API forGetDimensions.
Interface: call GetDimensions(nc, npmax, na)VariableIntentDatatypeDescriptionncoutinteger(4)Number of fluid componentsnpmaxoutinteger(4)Maximum number of fluid phasesnaoutinteger(4)Auxiliary array size (can be 0)Table 3
API forGetGlobalParameters.
Interface: call GetGlobalParameters(C-Name, M, P-Name, A-Name, units, opt)VariableIntentDatatypeDescriptionC-Nameoutcharacter3-byte names of the fluid components (arrayofsize3×nc)Moutreal(8)Molar weights of the components (array of sizenc), kg/molP-Nameoutcharacter3-byte names of the fluid phases (arrayofsize3×npmax)A-Nameoutcharacter8-byte names of the auxiliary arrays (arrayofsize8×na); any such name must begin with the “#” characterunitsoutcharacterUnits of the auxiliary arrays (arrayofsize8×na)optoutinteger(1)Options (array of size 8)Table 4
API forPhaseEquilibrium.
Interface call PhaseEquilibrium(P, T, ct, np, Φ, Id, a, Mode)VariableIntentDatatypeDescriptionPinreal(8)Average fluid pressure, PaTinreal(8)Fluid temperature, degree Kctinreal(8)Bulk mass composition (array of size nc)npin/outinteger(4)Number of fluid phasesΦin/outreal(8)Parameters of the phases (2-D array of size (6 + nc) ×npmax); only the first np columns of the array are relevant, i.e. Φi,i=1,…,npIdin/outinteger(1)Phase identifier (array of sizenpmax); only the first np elements of the array are relevant, i.e. Idi,i=1,…,npain/outreal(8)Auxiliary variables (array of sizena)Modeininteger(1)Calculation mode


### 2.3. Keywords and Mnemonics Associated with USEREOS

The input data to MUFITS are specified by the keywords, i.e., the commands to the simulator entered in the input data file. There is just one additional keyword USEREOS associated with the developed modelling option. This keyword activates the developed option. Its syntax is given in Table5.Table 5
The USEREOS keyword syntax. It takes two arguments: the paths to the library and configuration files.
USEREOS‘Shared library’ ‘Configuration file’ /The keyword USEREOS must be specified in the first section, RUNSPEC, of the input file. It takes just two arguments: the paths to the shared library and the configuration file. By default, the path to the library is USEREOS.DLL on Windows and USEREOS.DYLIB on Mac OS. The default path to the configuration file is an empty string, i.e., the library does not require such a file. When the simulator encounters the keyword, it searches for the compiled shared library and checks that it exports the four subroutines. If the library or the subroutines cannot be found, then the simulation is terminated with error. Otherwise, the simulator calls theReadConfigurationFile procedure, passing the path to the configuration file.A mnemonic is a short abbreviation of a reservoir model parameter. In the input data file, one can refer to a parameter by its mnemonic. There are eight additional mnemonics associated with USEREOS. These mnemonics are listed in Table6, where “C-Name(j)” and “P-Namei” are the names of the jth component and ith phase specified by the subroutine GetGlobalParameters. For example, if a component name is “CO2” and a phase name is “LIQ,” then CCO2-LIQ is the mass fraction of CO2 in the phase LIQ, DLIQ is the density of LIQ, and ZFCO2 is the bulk mass fraction of CO2.Table 6
Additional MUFITS mnemonics associated with the USEREOS option.
ParameterMnemoniccijC ‘C-Name(j)’-‘P-Namei’ρiD‘P-Namei’hiENTHM‘P-Namei’PiP‘P-Namei’KriREL‘P-Namei’siS‘P-Namei’μiVIS‘P-Namei’ctjZF‘C-Name(j)’


## 2.1. Overview

### 2.1.1. Data Flows

A sketch of data flows associated with any user-supplied EoS-module is shown in Figure2. Further, we refer to such a module as USEREOS. Figure 2(a) outlines the four mandatory program procedures of the module. Three of them are the preprocessing routines called to choose the modelling options and specify the arrays’ sizes and the general (i.e., global) parameters of a fluid model. The fourth routine, named PhaseEquilibrium, supplies the kernel with fluid parameters for a given P, T, and ct. This routine is called many times for any cell (i.e., grid block) of the reservoir model during time stepping. The coefficients of the fluid model are loaded directly into USEREOS by reading the configuration file.Figure 2
Sketch of data flows associated with USEREOS. The green and red dashed boxes in the right panel show the global and local parameters, respectively.Figure2(b) summarises the parameters and arrays in the kernel that are changed by or loaded into USEREOS. The data are gathered in four main blocks. The first block contains all the constants controlling memory usage. The second block, “Fluid model,” contains the parameters characterising the fluid and its phase and chemical equilibria. Some of these parameters that are shown within the green dashed box are global, i.e., they are identical for all the grid cells. Other parameters of the fluid that are within the red dashed box are local, i.e. they may be different in different cells. The third data block is for an auxiliary array that can be used as additional memory provided in the kernel for every cell. Again, this data block has global and local parameters. The fourth data block contains modelling options.


### 2.1.2. Constants Controlling Memory Usage

There are just three constants controlling memory allocation in the kernel.(i)
nc is the number of fluid components (see also Equations (1)–(8))(ii)
npmax is the maximum number of fluid phases. The number of fluid phases np depends on P, T, and ct. Thus, it can be different in different grid cells. However, for any P, T, and ct, the number of phases exported by USEREOS must not be larger than npmax, i.e., np≤npmax(iii)
na is the size of the auxiliary array. It specifies the number of additional 8-byte real numbers allocated in the kernel for every grid cell


### 2.1.3. Parameters of the Fluid Model

Every fluid componentj=1,⋯,nc is characterised by two global parameters:
(i)
C‐Name(j) is a character name of the component (e.g., ‘H2O’ or ‘CO2’)(ii)
Mj is the molar weight of the component, which is used in the kernel for calculating molar quantitiesEvery fluid phasei=1,…,np is characterised by one global and 6+nc local parameters. The global parameter P‐Namei is a character phase name (e.g., “GAS” or “LIQ”). The local parameters correspond to the vector
(9)Φi=ρi,hi,μi,si,Kri,ΔPi,ci1,⋯,cinc,where
(10)ΔPi=Pi−P,is the difference between the phase fluid pressure Pi and the average pressure P. The variable ΔPi is used for modelling the influence of capillary pressure Pc. For example, if the fluid can split into just two phases, liquid (i=1) and gas (i=2), then one can choose P=P1, ΔP1=0, and ΔP2=P2–P1=Pc,21. For every Φi, the specific energy of every phase is calculated in the kernel as
(11)ei=hi−Piρi.


### 2.1.4. Auxiliary Array

The auxiliary array of sizena is allocated in the kernel for every grid cell. This array can be used as additional memory provided for USEREOS. Some parameters can be stored in this array for later usage in the following time steps. For example, the auxiliary array can be used for modelling the reaction kinetics, the hysteresis phenomena, or calculating additional (secondary) parameters of the fluid not belonging to Φi (see Equation (9)). The latter is used when coupling MUFITS with GEMS (see Section 4).The elements of the auxiliary array can also be saved in the summary files in every report time period for later usage in the postprocessing of the simulation results. To make such reporting possible, every elementk=1,⋯,na of the array is characterised by two global parameters:
(i)
A‐Namek is a character name of the element that must begin with the symbol “#.” One can refer to the element by its name in the input data to MUFITS(ii)
Unitsk is a string specifying the units of the physical quantity stored in the element. These units are reported into the summary files. For dimensionless quantities, Unitsk should be an empty string or “NODIM.”Every elementk of the array is also characterised by one numeric quantity ak of type “double” (8-byte real number). The variable ak is local, i.e., it can be changed in USEREOS in every grid cell and at every time step.


### 2.1.5. The Option Array

The option array is an integer array of length 8,optl, where l=1,⋯,8. If optl=1, then the corresponding option is enabled (on). If optl=0, then the option is disabled (off). By default, all the options are disabled. At present, only opt1 and opt2 can be used, whereas the other elements of the option array, optl,l=3,⋯,8, are reserved for future developments.The optionopt1 controls the calculation of the relative permeability Kri. If opt1=0, then the relative permeability is calculated in the kernel using standard input data to the simulator (e.g., specified with the SATTAB keyword). In this case, there is no need to export Kri from USEREOS, i.e., the 5th element of Φi can be an undefined quantity. If opt1=1, then the default treatment of relative permeability is overridden by the value Kri exported from USEREOS.The optionopt2 controls the calculation of the capillary pressure Pc or, to be strict, the relative phase pressure ΔPi (see Equation (10)). If opt2=0, then the relative pressure is calculated in the kernel using standard input data to the simulator (e.g., provided with the SATTAB keyword). In this case, there is no need to export ΔPi from USEREOS, i.e., the 6th element of Φi can be an undefined quantity. If opt2=1, then the default treatment of capillary pressure is overridden by the value ΔPi exported from USEREOS.


### 2.1.6. Configuration File

The configuration file is an input data file to USEREOS. The path to this file is imported from the kernel. It is specified as the second argument of the USEREOS keyword that activates the presented modelling option (see Section2.3).The format of the configuration file can be specific for every particular shared library, because the kernel does not read this file. Every new USEREOS module can be supplemented with its own format of the file.It is implied that the configuration file contains coefficients of the fluid model that are not explicitly defined as static variables in USEREOS. The configuration file can contain paths to other external files loaded into USEREOS.


## 2.1.1. Data Flows

A sketch of data flows associated with any user-supplied EoS-module is shown in Figure2. Further, we refer to such a module as USEREOS. Figure 2(a) outlines the four mandatory program procedures of the module. Three of them are the preprocessing routines called to choose the modelling options and specify the arrays’ sizes and the general (i.e., global) parameters of a fluid model. The fourth routine, named PhaseEquilibrium, supplies the kernel with fluid parameters for a given P, T, and ct. This routine is called many times for any cell (i.e., grid block) of the reservoir model during time stepping. The coefficients of the fluid model are loaded directly into USEREOS by reading the configuration file.Figure 2
Sketch of data flows associated with USEREOS. The green and red dashed boxes in the right panel show the global and local parameters, respectively.Figure2(b) summarises the parameters and arrays in the kernel that are changed by or loaded into USEREOS. The data are gathered in four main blocks. The first block contains all the constants controlling memory usage. The second block, “Fluid model,” contains the parameters characterising the fluid and its phase and chemical equilibria. Some of these parameters that are shown within the green dashed box are global, i.e., they are identical for all the grid cells. Other parameters of the fluid that are within the red dashed box are local, i.e. they may be different in different cells. The third data block is for an auxiliary array that can be used as additional memory provided in the kernel for every cell. Again, this data block has global and local parameters. The fourth data block contains modelling options.


## 2.1.2. Constants Controlling Memory Usage

There are just three constants controlling memory allocation in the kernel.(i)
nc is the number of fluid components (see also Equations (1)–(8))(ii)
npmax is the maximum number of fluid phases. The number of fluid phases np depends on P, T, and ct. Thus, it can be different in different grid cells. However, for any P, T, and ct, the number of phases exported by USEREOS must not be larger than npmax, i.e., np≤npmax(iii)
na is the size of the auxiliary array. It specifies the number of additional 8-byte real numbers allocated in the kernel for every grid cell


## 2.1.3. Parameters of the Fluid Model

Every fluid componentj=1,⋯,nc is characterised by two global parameters:
(i)
C‐Name(j) is a character name of the component (e.g., ‘H2O’ or ‘CO2’)(ii)
Mj is the molar weight of the component, which is used in the kernel for calculating molar quantitiesEvery fluid phasei=1,…,np is characterised by one global and 6+nc local parameters. The global parameter P‐Namei is a character phase name (e.g., “GAS” or “LIQ”). The local parameters correspond to the vector
(9)Φi=ρi,hi,μi,si,Kri,ΔPi,ci1,⋯,cinc,where
(10)ΔPi=Pi−P,is the difference between the phase fluid pressure Pi and the average pressure P. The variable ΔPi is used for modelling the influence of capillary pressure Pc. For example, if the fluid can split into just two phases, liquid (i=1) and gas (i=2), then one can choose P=P1, ΔP1=0, and ΔP2=P2–P1=Pc,21. For every Φi, the specific energy of every phase is calculated in the kernel as
(11)ei=hi−Piρi.


## 2.1.4. Auxiliary Array

The auxiliary array of sizena is allocated in the kernel for every grid cell. This array can be used as additional memory provided for USEREOS. Some parameters can be stored in this array for later usage in the following time steps. For example, the auxiliary array can be used for modelling the reaction kinetics, the hysteresis phenomena, or calculating additional (secondary) parameters of the fluid not belonging to Φi (see Equation (9)). The latter is used when coupling MUFITS with GEMS (see Section 4).The elements of the auxiliary array can also be saved in the summary files in every report time period for later usage in the postprocessing of the simulation results. To make such reporting possible, every elementk=1,⋯,na of the array is characterised by two global parameters:
(i)
A‐Namek is a character name of the element that must begin with the symbol “#.” One can refer to the element by its name in the input data to MUFITS(ii)
Unitsk is a string specifying the units of the physical quantity stored in the element. These units are reported into the summary files. For dimensionless quantities, Unitsk should be an empty string or “NODIM.”Every elementk of the array is also characterised by one numeric quantity ak of type “double” (8-byte real number). The variable ak is local, i.e., it can be changed in USEREOS in every grid cell and at every time step.


## 2.1.5. The Option Array

The option array is an integer array of length 8,optl, where l=1,⋯,8. If optl=1, then the corresponding option is enabled (on). If optl=0, then the option is disabled (off). By default, all the options are disabled. At present, only opt1 and opt2 can be used, whereas the other elements of the option array, optl,l=3,⋯,8, are reserved for future developments.The optionopt1 controls the calculation of the relative permeability Kri. If opt1=0, then the relative permeability is calculated in the kernel using standard input data to the simulator (e.g., specified with the SATTAB keyword). In this case, there is no need to export Kri from USEREOS, i.e., the 5th element of Φi can be an undefined quantity. If opt1=1, then the default treatment of relative permeability is overridden by the value Kri exported from USEREOS.The optionopt2 controls the calculation of the capillary pressure Pc or, to be strict, the relative phase pressure ΔPi (see Equation (10)). If opt2=0, then the relative pressure is calculated in the kernel using standard input data to the simulator (e.g., provided with the SATTAB keyword). In this case, there is no need to export ΔPi from USEREOS, i.e., the 6th element of Φi can be an undefined quantity. If opt2=1, then the default treatment of capillary pressure is overridden by the value ΔPi exported from USEREOS.


## 2.1.6. Configuration File

The configuration file is an input data file to USEREOS. The path to this file is imported from the kernel. It is specified as the second argument of the USEREOS keyword that activates the presented modelling option (see Section2.3).The format of the configuration file can be specific for every particular shared library, because the kernel does not read this file. Every new USEREOS module can be supplemented with its own format of the file.It is implied that the configuration file contains coefficients of the fluid model that are not explicitly defined as static variables in USEREOS. The configuration file can contain paths to other external files loaded into USEREOS.


## 2.2. Subroutines Exported from USEREOS

### 2.2.1. General Description

In this section, we give a brief overview of the main subroutines exported from USEREOS. The four subroutines placed in the order of their invoking by the kernel are as follows:(i)
ReadConfigurationFile. This routine is called by the kernel at the initialisation stage, immediately after the USEREOS keyword is encountered in the input data file to MUFITS (Section 2.3). The path to the configuration file is the second argument of the keyword that is imported by USEREOS as an argument of the subroutine. It is implied that the shared library reads the file and then exports the title of the module to the kernel. The title may contain a brief description of the module and its version. The title is reported to the LOG-file of MUFITS.(ii)
GetDimensions. This subroutine is called at the end of initialisation stage. It is aimed at importing the sizes of arrays that need to be allocated into the kernel. The subroutine exports the number of fluid components nc, the maximum number of fluid phases (npmax), and the auxiliary array size (na).(iii)
GetGlobalParameters. The procedure is called immediately after the previous one. It is aimed at loading the global parameters of the fluid model into the kernel. First, the 3-byte character name (C‐Name(j)) and the molar weight (Mj) of every component j=1,…,nc must be specified. Second, the 3-byte character name (P‐Namei) for each fluid phase i=1,…,npmax is exported by USEREOS. Third, the 8-byte character name (A-Namek) and units (Unitsk) for every element k=1,…,na of the auxiliary array must be specified. Any A-Namek must begin with the “#” character. Fourth, the options array, optl, is exported.(iv)
PhaseEquilibrium. Normally, this routine is called multiple times for every grid cell during linearisation of the governing Equations (1)–(8) and time stepping. It is aimed at reporting the phase equilibrium for a given P, T, and ct. It reports the number of phases np≤npmax, the vector Φi, and the integer variable Idi for every phase i=1,…,np. The quantity Idi (the “phase identifier”) is aimed at associating the ith phase with the phase names (i.e., phase types) specified by GetGlobalParameters. The value of Idi must be equal to the serial number of the P-Namek, k=1,…,npmax specified by GetGlobalParameters.A very important argument ofPhaseEquilibrium is Mode. If Mode=0, then the subroutine is supplied by the kernel with the initial guess for the phase equilibrium. The variables np, Idi, and Φi, i=1,…,np, are used to specify the initial guess. In the case of Mode=0, PhaseEquilibrium must calculate the equilibrium using the provided initial guess. The procedure must not change np and Idi. It can only change Φi for a given P, T, and ctj. Since the number of phases np is not changed, the reported saturations si and concentrations cij and ctj are allowed to take negative or larger-than-one values [19]. If Mode=1, then the initial guess is not specified. In this case, the subroutine must calculate the phase equilibrium—i.e., np, Idi, and Φi—from scratch. If Mode=1, then the reported si and cij must satisfy the inequalities 0≤si≤1 and 0≤cij≤1, respectively.


### 2.2.2. Calling Interface for the USEREOS Subroutines

The USEREOS API is provided in the Fortran programming language. However, USEREOS can be written in any language that supports compilation to a shared library. To export the subroutines described in Section2.2.1, it may be necessary to change these subroutines’ signatures to accommodate differences between programming languages. Since GEMS3K is written in C++, the USEREOS module for coupling GEMS3K and MUFITS was written in C++ as well (Section 4).The Fortran API for the four primary subroutines exported by USEREOS is given in Tables1–4. In the tables, every argument of each subroutine is supplemented with a brief description, its datatype, and the intent. The intent values “in” or “out” mean that the corresponding variable is imported to or exported from USEREOS, respectively. The value “in/out” means that variable is used as both an input and output parameter. The datatypes are shown in a standard format used in Fortran. For example, the datatypes “integer(4)” and “real(8)” correspond to 4-byte integer and 8-byte real numbers, respectively. All physical quantities must be given in SI units. For example, the units of ρi, hi, μi, and ΔPi in the vector Φi must be kg/m3, J/kg, Pa∙s, and Pa, respectively.Table 1
API forReadConfigurationFile.
Interface: call ReadConfigurationFile (filename ,ierr, title)VariableIntentDatatypeDescriptionfilenameincharacterConfiguration file name (array of size 256)ierroutinteger(4)Error specifier (ifierr=0, then no error occurred during subroutine execution; if ierr=−1, then the configuration file is not found)titleoutcharacterEoS-module title, i.e., a brief description (array of size 80)Table 2
API forGetDimensions.
Interface: call GetDimensions(nc, npmax, na)VariableIntentDatatypeDescriptionncoutinteger(4)Number of fluid componentsnpmaxoutinteger(4)Maximum number of fluid phasesnaoutinteger(4)Auxiliary array size (can be 0)Table 3
API forGetGlobalParameters.
Interface: call GetGlobalParameters(C-Name, M, P-Name, A-Name, units, opt)VariableIntentDatatypeDescriptionC-Nameoutcharacter3-byte names of the fluid components (arrayofsize3×nc)Moutreal(8)Molar weights of the components (array of sizenc), kg/molP-Nameoutcharacter3-byte names of the fluid phases (arrayofsize3×npmax)A-Nameoutcharacter8-byte names of the auxiliary arrays (arrayofsize8×na); any such name must begin with the “#” characterunitsoutcharacterUnits of the auxiliary arrays (arrayofsize8×na)optoutinteger(1)Options (array of size 8)Table 4
API forPhaseEquilibrium.
Interface call PhaseEquilibrium(P, T, ct, np, Φ, Id, a, Mode)VariableIntentDatatypeDescriptionPinreal(8)Average fluid pressure, PaTinreal(8)Fluid temperature, degree Kctinreal(8)Bulk mass composition (array of size nc)npin/outinteger(4)Number of fluid phasesΦin/outreal(8)Parameters of the phases (2-D array of size (6 + nc) ×npmax); only the first np columns of the array are relevant, i.e. Φi,i=1,…,npIdin/outinteger(1)Phase identifier (array of sizenpmax); only the first np elements of the array are relevant, i.e. Idi,i=1,…,npain/outreal(8)Auxiliary variables (array of sizena)Modeininteger(1)Calculation mode


## 2.2.1. General Description

In this section, we give a brief overview of the main subroutines exported from USEREOS. The four subroutines placed in the order of their invoking by the kernel are as follows:(i)
ReadConfigurationFile. This routine is called by the kernel at the initialisation stage, immediately after the USEREOS keyword is encountered in the input data file to MUFITS (Section 2.3). The path to the configuration file is the second argument of the keyword that is imported by USEREOS as an argument of the subroutine. It is implied that the shared library reads the file and then exports the title of the module to the kernel. The title may contain a brief description of the module and its version. The title is reported to the LOG-file of MUFITS.(ii)
GetDimensions. This subroutine is called at the end of initialisation stage. It is aimed at importing the sizes of arrays that need to be allocated into the kernel. The subroutine exports the number of fluid components nc, the maximum number of fluid phases (npmax), and the auxiliary array size (na).(iii)
GetGlobalParameters. The procedure is called immediately after the previous one. It is aimed at loading the global parameters of the fluid model into the kernel. First, the 3-byte character name (C‐Name(j)) and the molar weight (Mj) of every component j=1,…,nc must be specified. Second, the 3-byte character name (P‐Namei) for each fluid phase i=1,…,npmax is exported by USEREOS. Third, the 8-byte character name (A-Namek) and units (Unitsk) for every element k=1,…,na of the auxiliary array must be specified. Any A-Namek must begin with the “#” character. Fourth, the options array, optl, is exported.(iv)
PhaseEquilibrium. Normally, this routine is called multiple times for every grid cell during linearisation of the governing Equations (1)–(8) and time stepping. It is aimed at reporting the phase equilibrium for a given P, T, and ct. It reports the number of phases np≤npmax, the vector Φi, and the integer variable Idi for every phase i=1,…,np. The quantity Idi (the “phase identifier”) is aimed at associating the ith phase with the phase names (i.e., phase types) specified by GetGlobalParameters. The value of Idi must be equal to the serial number of the P-Namek, k=1,…,npmax specified by GetGlobalParameters.A very important argument ofPhaseEquilibrium is Mode. If Mode=0, then the subroutine is supplied by the kernel with the initial guess for the phase equilibrium. The variables np, Idi, and Φi, i=1,…,np, are used to specify the initial guess. In the case of Mode=0, PhaseEquilibrium must calculate the equilibrium using the provided initial guess. The procedure must not change np and Idi. It can only change Φi for a given P, T, and ctj. Since the number of phases np is not changed, the reported saturations si and concentrations cij and ctj are allowed to take negative or larger-than-one values [19]. If Mode=1, then the initial guess is not specified. In this case, the subroutine must calculate the phase equilibrium—i.e., np, Idi, and Φi—from scratch. If Mode=1, then the reported si and cij must satisfy the inequalities 0≤si≤1 and 0≤cij≤1, respectively.


## 2.2.2. Calling Interface for the USEREOS Subroutines

The USEREOS API is provided in the Fortran programming language. However, USEREOS can be written in any language that supports compilation to a shared library. To export the subroutines described in Section2.2.1, it may be necessary to change these subroutines’ signatures to accommodate differences between programming languages. Since GEMS3K is written in C++, the USEREOS module for coupling GEMS3K and MUFITS was written in C++ as well (Section 4).The Fortran API for the four primary subroutines exported by USEREOS is given in Tables1–4. In the tables, every argument of each subroutine is supplemented with a brief description, its datatype, and the intent. The intent values “in” or “out” mean that the corresponding variable is imported to or exported from USEREOS, respectively. The value “in/out” means that variable is used as both an input and output parameter. The datatypes are shown in a standard format used in Fortran. For example, the datatypes “integer(4)” and “real(8)” correspond to 4-byte integer and 8-byte real numbers, respectively. All physical quantities must be given in SI units. For example, the units of ρi, hi, μi, and ΔPi in the vector Φi must be kg/m3, J/kg, Pa∙s, and Pa, respectively.Table 1
API forReadConfigurationFile.
Interface: call ReadConfigurationFile (filename ,ierr, title)VariableIntentDatatypeDescriptionfilenameincharacterConfiguration file name (array of size 256)ierroutinteger(4)Error specifier (ifierr=0, then no error occurred during subroutine execution; if ierr=−1, then the configuration file is not found)titleoutcharacterEoS-module title, i.e., a brief description (array of size 80)Table 2
API forGetDimensions.
Interface: call GetDimensions(nc, npmax, na)VariableIntentDatatypeDescriptionncoutinteger(4)Number of fluid componentsnpmaxoutinteger(4)Maximum number of fluid phasesnaoutinteger(4)Auxiliary array size (can be 0)Table 3
API forGetGlobalParameters.
Interface: call GetGlobalParameters(C-Name, M, P-Name, A-Name, units, opt)VariableIntentDatatypeDescriptionC-Nameoutcharacter3-byte names of the fluid components (arrayofsize3×nc)Moutreal(8)Molar weights of the components (array of sizenc), kg/molP-Nameoutcharacter3-byte names of the fluid phases (arrayofsize3×npmax)A-Nameoutcharacter8-byte names of the auxiliary arrays (arrayofsize8×na); any such name must begin with the “#” characterunitsoutcharacterUnits of the auxiliary arrays (arrayofsize8×na)optoutinteger(1)Options (array of size 8)Table 4
API forPhaseEquilibrium.
Interface call PhaseEquilibrium(P, T, ct, np, Φ, Id, a, Mode)VariableIntentDatatypeDescriptionPinreal(8)Average fluid pressure, PaTinreal(8)Fluid temperature, degree Kctinreal(8)Bulk mass composition (array of size nc)npin/outinteger(4)Number of fluid phasesΦin/outreal(8)Parameters of the phases (2-D array of size (6 + nc) ×npmax); only the first np columns of the array are relevant, i.e. Φi,i=1,…,npIdin/outinteger(1)Phase identifier (array of sizenpmax); only the first np elements of the array are relevant, i.e. Idi,i=1,…,npain/outreal(8)Auxiliary variables (array of sizena)Modeininteger(1)Calculation mode


## 2.3. Keywords and Mnemonics Associated with USEREOS

The input data to MUFITS are specified by the keywords, i.e., the commands to the simulator entered in the input data file. There is just one additional keyword USEREOS associated with the developed modelling option. This keyword activates the developed option. Its syntax is given in Table5.Table 5
The USEREOS keyword syntax. It takes two arguments: the paths to the library and configuration files.
USEREOS‘Shared library’ ‘Configuration file’ /The keyword USEREOS must be specified in the first section, RUNSPEC, of the input file. It takes just two arguments: the paths to the shared library and the configuration file. By default, the path to the library is USEREOS.DLL on Windows and USEREOS.DYLIB on Mac OS. The default path to the configuration file is an empty string, i.e., the library does not require such a file. When the simulator encounters the keyword, it searches for the compiled shared library and checks that it exports the four subroutines. If the library or the subroutines cannot be found, then the simulation is terminated with error. Otherwise, the simulator calls theReadConfigurationFile procedure, passing the path to the configuration file.A mnemonic is a short abbreviation of a reservoir model parameter. In the input data file, one can refer to a parameter by its mnemonic. There are eight additional mnemonics associated with USEREOS. These mnemonics are listed in Table6, where “C-Name(j)” and “P-Namei” are the names of the jth component and ith phase specified by the subroutine GetGlobalParameters. For example, if a component name is “CO2” and a phase name is “LIQ,” then CCO2-LIQ is the mass fraction of CO2 in the phase LIQ, DLIQ is the density of LIQ, and ZFCO2 is the bulk mass fraction of CO2.Table 6
Additional MUFITS mnemonics associated with the USEREOS option.
ParameterMnemoniccijC ‘C-Name(j)’-‘P-Namei’ρiD‘P-Namei’hiENTHM‘P-Namei’PiP‘P-Namei’KriREL‘P-Namei’siS‘P-Namei’μiVIS‘P-Namei’ctjZF‘C-Name(j)’


## 3. Simple Application Examples

To demonstrate the software development and ease its usage, we further present two simple examples of USEREOS and their validation. The examples are supplemented with the source code and compiled executables [20, 21].


### 3.1. Modelling Ideal Gas Flow through a Porous Medium

The first example concerns the numerical modelling of nonisothermal single-phase flows of ideal gas. The corresponding module can account for gases consisting of an arbitrary number of components. The fluid properties are calculated using the following equations:(12)ρ=PMRT,M=∑j=1ncMjxj,(13)h′=∑j=1ncCj′xjT,μ=∑j=1ncμjxj,(14)s=1,cj=Mjxj∑k=1ncMkxk,where R is the universal gas constant, x is the molar concentration, h′ and C′ are the molar enthalpy, and heat capacity at P=const, respectively. For simplicity, we omit the subscript i≡1 of the only phase of the fluid.Using the developed module, we present a 1-D study for a three-component gas composed of N2, O2, and CH4. We consider a 1-D flow in a horizontal homogeneous reservoir X∈0,100 m, where X is the length. The porosity is ϕ=0.2, and the permeability is K=100 mD. Other relevant parameters are the rock density, ρr=2500 kg/m3, and the specific rock heat capacity, Cr=1 kJ/(kg∙K). At t=0, the porous medium is saturated with the gas of composition ct=0.1,0.0,0.9, i.e., it consists of 10% of N2 and 90% of CH4. The initial reservoir pressure and temperature are 200 bar and 50°C, respectively. The Dirichlet boundary conditions are imposed at X=100 m, i.e., the initial pressure (and temperature) are kept constant at the open boundary X=100 m. A gas of different composition ct=0.77,0.23,0.0, i.e., air, is injected into the reservoir through the boundary X=0 m. The injected gas temperature is 30°C at P=200 bar. The injection rate is kept constant at 0.2 m/day at reservoir pressure. The permeability K is so large that the pressure deviation from its initial value does not exceed 0.5 bar. Thus, the volume injection rate can be assumed at the initial reservoir pressure. Other relevant parameters of the fluid are summarised in Table 7.Table 7
Parameters of the ideal gas mixture.
N2O2CH4Mj, g/mol283216Cj′, J/(mol∙K)29.1229.4435.2μj, cP0.0160.0220.012Assuming that all dispersion effects can be neglected, i.e.,ψ=0, the formulated problem admits a simple analytical solution that includes two shocks Sc and ST (Figure 3). The concentrations ct are discontinuous at the leading shock Sc. The initial gas composition is preserved ahead of Sc, whereas the porous medium is saturated with the injected gas behind Sc. The reservoir is at the initial temperature T=50°C ahead of the trailing shock ST and at the injection temperature of 30°C behind ST. When constructing the analytical solution, one can derive that Sc and ST are propagating at constant velocities of 0.709 m/day and 0.0254 m/day, respectively. This piecewise constant distribution can be obtained as the solution to Equations (1)–(6) where ψ=0 and the fluid properties in Equation (4) are given by Equations (12)–(14). The velocities of the shocks are given by the Rankine-Hugoniot conditions corresponding to the balance Equations (1) and (2). The initial and boundary conditions are discussed in the previous paragraph.Figure 3
Distributions of (a)ct and (b) ρ and T at t=100 days. The solid and dotted curves show the numerical and analytical solutions, respectively.
(a)(b)As shown in Figure3, the simulated distributions agree with the analytical solution. However, because of numerical dispersion, the shocks are zones of continuous transition of finite extent. The results are obtained with a regular grid consisting of 1000 elements and the maximum time step of 0.25 days.


### 3.2. Modelling Two-Phase Immiscible Displacement

This article is also supplemented with an example of an external library designed for modelling two-phase immiscible displacement. Although this USEREOS module allows for nonisothermal modelling, we consider a simpler case of an isothermal flow to validate the module.In the module, the fluid properties are predicted using the following relations:(15)ρi=ρi,ref1+αiP−Pref−βiT−Tref,(16)hi=CiT−Tref,μi=const,(17)si=ρicti∑k=12ρkctk,(18)c11=1,c12=0,c21=1,c22=0,where α is the isothermal compressibility, β is the coefficient of thermal expansion, Pref and Tref are the reference values of pressure and temperature, respectively, and C is the specific heat capacity at P=const. According to Equation (18), the component j=k of the mixture forms the phase i=k, where k=1, 2. For example, the phase (component) i=1 is water, and the phase i=2 is oil.We consider a 1-D flow in a horizontal homogeneous reservoir atX∈0,100 m. The porosity and permeability of the medium are ϕ=0.2 and K=200 mD, respectively. The relative permeability is given by
(19)Kr1=s12,Kr2=s23.The capillary pressure is assumed to be 0. Att=0, the reservoir is saturated by the oil phase i=2, i.e., ct1=0, at P=Pref=200 bar and T=Tref. The reservoir temperature Tref is kept constant at all t>0. Therefore, its value as well as the values of βi and Ci does not influence the flow. The initial pressure is kept constant at the open boundary X=100 m. The water phase (i=1) is injected through the boundary X=0 m at constant injection rate of 0.1 m/day.The compressibilitiesαi are assumed to be small such that the phase densities ρi can be considered constant. Therefore, only the phase viscosities, μ1=0.75 cP and μ2=1.5 cP, are relevant parameters.The water injection results in the immiscible displacement of oil from the boundaryX=0 m to X=100 m (Figure 4). The saturation distribution is governed by the Buckley–Leverett solution [22]. This analytical solution consists of the leading displacement front S, at which s1 is discontinuous, and the trailing Riemann wave, in which s1 changes continuously, reaching s1=1 at X=0 m. One can show that for Equation (19) and specified phase viscosities μ1 and μ2, the shock propagates at a velocity of 0.8 m/day. Thus, S is exactly at X=80 m at t=100 days.Figure 4
Distribution of the water saturation att=100 days. The solid and dashed curves correspond to the numerical and analytical solutions, respectively.As shown in Figure4, the numerical solution agrees with the analytical solution to the Buckley–Leverett problem. The simulation was conducted with a regular grid consisting of 1000 cells. The time step was limited by 0.25 days to reduce truncation errors.


## 3.1. Modelling Ideal Gas Flow through a Porous Medium

The first example concerns the numerical modelling of nonisothermal single-phase flows of ideal gas. The corresponding module can account for gases consisting of an arbitrary number of components. The fluid properties are calculated using the following equations:(12)ρ=PMRT,M=∑j=1ncMjxj,(13)h′=∑j=1ncCj′xjT,μ=∑j=1ncμjxj,(14)s=1,cj=Mjxj∑k=1ncMkxk,where R is the universal gas constant, x is the molar concentration, h′ and C′ are the molar enthalpy, and heat capacity at P=const, respectively. For simplicity, we omit the subscript i≡1 of the only phase of the fluid.Using the developed module, we present a 1-D study for a three-component gas composed of N2, O2, and CH4. We consider a 1-D flow in a horizontal homogeneous reservoir X∈0,100 m, where X is the length. The porosity is ϕ=0.2, and the permeability is K=100 mD. Other relevant parameters are the rock density, ρr=2500 kg/m3, and the specific rock heat capacity, Cr=1 kJ/(kg∙K). At t=0, the porous medium is saturated with the gas of composition ct=0.1,0.0,0.9, i.e., it consists of 10% of N2 and 90% of CH4. The initial reservoir pressure and temperature are 200 bar and 50°C, respectively. The Dirichlet boundary conditions are imposed at X=100 m, i.e., the initial pressure (and temperature) are kept constant at the open boundary X=100 m. A gas of different composition ct=0.77,0.23,0.0, i.e., air, is injected into the reservoir through the boundary X=0 m. The injected gas temperature is 30°C at P=200 bar. The injection rate is kept constant at 0.2 m/day at reservoir pressure. The permeability K is so large that the pressure deviation from its initial value does not exceed 0.5 bar. Thus, the volume injection rate can be assumed at the initial reservoir pressure. Other relevant parameters of the fluid are summarised in Table 7.Table 7
Parameters of the ideal gas mixture.
N2O2CH4Mj, g/mol283216Cj′, J/(mol∙K)29.1229.4435.2μj, cP0.0160.0220.012Assuming that all dispersion effects can be neglected, i.e.,ψ=0, the formulated problem admits a simple analytical solution that includes two shocks Sc and ST (Figure 3). The concentrations ct are discontinuous at the leading shock Sc. The initial gas composition is preserved ahead of Sc, whereas the porous medium is saturated with the injected gas behind Sc. The reservoir is at the initial temperature T=50°C ahead of the trailing shock ST and at the injection temperature of 30°C behind ST. When constructing the analytical solution, one can derive that Sc and ST are propagating at constant velocities of 0.709 m/day and 0.0254 m/day, respectively. This piecewise constant distribution can be obtained as the solution to Equations (1)–(6) where ψ=0 and the fluid properties in Equation (4) are given by Equations (12)–(14). The velocities of the shocks are given by the Rankine-Hugoniot conditions corresponding to the balance Equations (1) and (2). The initial and boundary conditions are discussed in the previous paragraph.Figure 3
Distributions of (a)ct and (b) ρ and T at t=100 days. The solid and dotted curves show the numerical and analytical solutions, respectively.
(a)(b)As shown in Figure3, the simulated distributions agree with the analytical solution. However, because of numerical dispersion, the shocks are zones of continuous transition of finite extent. The results are obtained with a regular grid consisting of 1000 elements and the maximum time step of 0.25 days.


## 3.2. Modelling Two-Phase Immiscible Displacement

This article is also supplemented with an example of an external library designed for modelling two-phase immiscible displacement. Although this USEREOS module allows for nonisothermal modelling, we consider a simpler case of an isothermal flow to validate the module.In the module, the fluid properties are predicted using the following relations:(15)ρi=ρi,ref1+αiP−Pref−βiT−Tref,(16)hi=CiT−Tref,μi=const,(17)si=ρicti∑k=12ρkctk,(18)c11=1,c12=0,c21=1,c22=0,where α is the isothermal compressibility, β is the coefficient of thermal expansion, Pref and Tref are the reference values of pressure and temperature, respectively, and C is the specific heat capacity at P=const. According to Equation (18), the component j=k of the mixture forms the phase i=k, where k=1, 2. For example, the phase (component) i=1 is water, and the phase i=2 is oil.We consider a 1-D flow in a horizontal homogeneous reservoir atX∈0,100 m. The porosity and permeability of the medium are ϕ=0.2 and K=200 mD, respectively. The relative permeability is given by
(19)Kr1=s12,Kr2=s23.The capillary pressure is assumed to be 0. Att=0, the reservoir is saturated by the oil phase i=2, i.e., ct1=0, at P=Pref=200 bar and T=Tref. The reservoir temperature Tref is kept constant at all t>0. Therefore, its value as well as the values of βi and Ci does not influence the flow. The initial pressure is kept constant at the open boundary X=100 m. The water phase (i=1) is injected through the boundary X=0 m at constant injection rate of 0.1 m/day.The compressibilitiesαi are assumed to be small such that the phase densities ρi can be considered constant. Therefore, only the phase viscosities, μ1=0.75 cP and μ2=1.5 cP, are relevant parameters.The water injection results in the immiscible displacement of oil from the boundaryX=0 m to X=100 m (Figure 4). The saturation distribution is governed by the Buckley–Leverett solution [22]. This analytical solution consists of the leading displacement front S, at which s1 is discontinuous, and the trailing Riemann wave, in which s1 changes continuously, reaching s1=1 at X=0 m. One can show that for Equation (19) and specified phase viscosities μ1 and μ2, the shock propagates at a velocity of 0.8 m/day. Thus, S is exactly at X=80 m at t=100 days.Figure 4
Distribution of the water saturation att=100 days. The solid and dashed curves correspond to the numerical and analytical solutions, respectively.As shown in Figure4, the numerical solution agrees with the analytical solution to the Buckley–Leverett problem. The simulation was conducted with a regular grid consisting of 1000 cells. The time step was limited by 0.25 days to reduce truncation errors.


## 4. Coupling with GEMS3K

This more elaborate example concerns an EoS-module that can be used as a template for MUFITS coupling with a geochemical software. We couple the simulator with the numerical code GEMS3K to simulate reactive transport in porous media. GEMS3K is a C++ library implementing the efficient IPM-3 algorithm for GEM [17]. Previously, it was coupled with various transport codes such as CSMP++, OpenGeoSys, COMSOL, and other software packages [23–25].In the developed USEREOS, the representative elementary volume of a porous medium consists of the inert rock, whose properties are calculated in the kernel, and the pore fluid described by GEMS3K (Figure5). The fluid can split into different phases, including liquid and gas, that can move through the porous medium and several precipitated minerals, i.e., solid phases. Thus, the rock consists of the inert part, which does not participate in chemical reactions, and several minerals that can precipitate from or dissolve in the moving phases. We assume that the reactive transport occurs under a local chemical equilibrium between fluid and rock. Thus, the equilibrium is characterised by P, T, and the bulk concentrations ct of the primary chemical elements that the fluid is composed of. For a given P, T, and ct, GEMS3K calculates the phase equilibrium, including the number of phases np and the mole amounts of the chemical species in the phases. The latter quantities are translated by USEREOS to the mass concentrations cij. The chemical speciation is stored in the auxiliary array discussed in Section 2.1.4.Figure 5
Sketch of the representative elementary volume. The fluid, i.e., the moving phases and precipitated minerals, occupies the pore space. The rest of the volume is regarded as inert rock.The configuration files (see Section2.1.6) for the developed USEREOS are identical to those for GEMS3K. They can be prepared in GEM-Selektor, a program for geochemical modelling with a graphical user interface, which provides access to various thermodynamic databases and allows for exporting parameters of the fluid model and numerical controls in a format directly compatible with GEMS3K.To validate and demonstrate the developed USEREOS, we consider a 1-D benchmark problem of mineral dissolution–precipitation that is often used for the validation of reactive transport modelling [8, 23, 24, 26]. Here, we follow the problem statement presented in [23, 24]. We consider a 1-D porous column located at X∈0,0.5 m, where ϕ=0.32. At t=0, the column is saturated with an aqueous solution in equilibrium with calcite at P=1 bar and T=25°C. The initial pressure is kept constant at the open boundary X=0.5 m. The column is flushed through the boundary X=0 with a MgCl2 solution at a constant flow rate q=3∙10−6 m/s. The initial and boundary elemental concentrations are listed in Table 8. The flow is assumed to be isothermal at T=25°C. The pressure in [23, 24] is also fixed at 1 bar. In MUFTIS, the presence of a pressure gradient is necessary to produce a flow according to Equation (3), but for the thermodynamic equilibria calculations in GEMS3K, the pressure is set to 1 bar to achieve a better match to the reference results.Table 8
Fluid composition (mass fractions) at the inflow boundaryX=0 and at t=0.
BoundaryInitialC1.20e-123.96e-06Ca4.01e-121.32e-05Cl7.09e-057.09e-09H1.12e-011.12e-01Mg2.43e-052.43e-09O8.88e-018.88e-01The mechanical dispersion in this problem is modelled by Equations (7) and (8). Since the study is 1-D, only the longitudinal dispersion length, λL=0.00214 m, must be specified. This length corresponds to that used by Shao et al. [24]. However, because Shao et al. [24] use a different definition for the dispersion length, the quantity λL scales to that in [24] by the factor ϕ.To simulate the fluid flow, we use a fluid consisting of the six components listed in Table8. The simulated distributions of the concentrations of Cl, Mg, and Ca in the fluid and the concentrations of calcite and dolomite after 21000 s of MgCl2 injection are shown in Figure 6. The simulated distributions are in good overall agreement with the results presented in [24]. Herewith, dolomite is formed in a finite moving zone, and calcite dissolves at the leading reaction front. Minor differences between the concentration profiles can be the result of using different equations for the hydrodynamic part of the model and different coupling schemes. In [23, 24], the sequential noniterative approach (SNIA) is employed. This approach assumes that the solution of the transport equations is conducted separately from the chemical equilibrium calculations at each simulation time step. In this work, the transport and the chemical equilibria are calculated simultaneously at every time step (the global implicit approach). This ensures that both the mass balance and the assumption of the chemical equilibrium are satisfied at the end of each time step.Figure 6
Concentrations of the aqueous components and minerals att=21000 s. The solid curves show the results of Shao et al. [24], and the symbols are the results of the MUFITS–GEMS3K coupled simulation.


## 5. Conclusions

We demonstrate that the developed extension is robust and useful because it can be applied in various studies ranging from the petroleum reservoir simulations to nonisothermal flows and reactive transport modelling. By using the 1-D studies, we validate the correctness of the program implementation of numerical algorithms associated with USEREOS. Although the considered benchmarks are rather simple, USEREOS can be applied in more sophisticated 2-D and 3-D scenarios involving complicated phenomena and processes. These can include coupled reservoir-wellbore flows or transport in fractured–porous media, which can be modelled with the available capabilities of the software implemented in its kernel.The created capabilities for MUFITS coupling with other software for fluid property prediction deserve special attention. Aiming at reactive transport modelling, we consider only an example of coupling with GEMS3K. However, in the same way, the simulator can be coupled with other geochemical software, e.g., PHREEQC [27] or those implementing accelerated approaches [15, 28]. It is important that from now on, such coupling can be done independently by the software’s users. This provides the scientific community with an advanced and innovative opportunity for modelling subsurface flows, given that the developed API between MUFITS and USEREOS is relatively simple and easy to use.



---

*Source: 4290248-2021-09-29.xml*