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

Afanasyev, Andrey; Utkin, Ivan

doi:https://doi.org/10.1155/2021/4290248

Geofluids

On this page

Abstract Introduction Conclusions Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2021 | Article ID 4290248 | https://doi.org/10.1155/2021/4290248

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

Andrey Afanasyev¹and Ivan Utkin¹

Academic Editor: Jinze Xu

Received27 May 2021

Accepted16 Aug 2021

Published29 Sept 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.

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 CO₂ 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]:

Here, nc is the number of fluid components; np is the number of fluid phases; is the porosity; is the density; is the th component mass fraction in the th phase; is the phase saturation (i.e., volume fraction of pore space occupied by the phase); is the Darcy velocity; and are the fluxes of the th component and energy caused by mechanical dispersion, molecular diffusion, or heat conduction; is the specific energy; is the specific enthalpy; is the absolute permeability of the medium; is the relative permeability; is the dynamic viscosity; is the pressure; is the gravity acceleration; is the bulk mass fraction of the th component; is the temperature; is the capillary pressure; the subscripts , , and denote the parameters of the th phase, the th 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 , temperature , and bulk composition . Equation (5) defines the relative permeability and capillary pressure as saturation functions. Equations in (6) are additional relations. For simplicity, we assume that , , and are constants and that 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 should be assumed.

In the case of mechanical dispersion, the fluxes take the following form [12]: where is the tensor of mechanical dispersion; and are the longitudinal and transverse mean free path of a small volume of the th phase caused by mixing in the porous medium; and , , and 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).

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 (Figure 1). 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 (Figure 1). 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 (Figure 1). For example, it includes EoS-modules for the black-oil modelling of petroleum reservoirs and the subsurface storage of natural gas or CO₂ 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 Section 2, 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 Figure 2. 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 , , and . 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(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) is the maximum number of fluid phases. The number of fluid phases np depends on , , and . Thus, it can be different in different grid cells. However, for any , , and , the number of phases exported by USEREOS must not be larger than , i.e., (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 component is characterised by two global parameters: (i)C‐Name_(j) is a character name of the component (e.g., ‘H₂O’ or ‘CO₂’)(ii) is the molar weight of the component, which is used in the kernel for calculating molar quantities

Every fluid phase is characterised by one global and local parameters. The global parameter P‐Name_i is a character phase name (e.g., “GAS” or “LIQ”). The local parameters correspond to the vector where is the difference between the phase fluid pressure and the average pressure . The variable is used for modelling the influence of capillary pressure . For example, if the fluid can split into just two phases, liquid () and gas (), then one can choose , , and . For every , the specific energy of every phase is calculated in the kernel as

2.1.4. Auxiliary Array

The auxiliary array of size na 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 (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 element of the array is characterised by two global parameters: (i)A‐Name_k 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) 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, should be an empty string or “NODIM.”

Every element of the array is also characterised by one numeric quantity of type “double” (8-byte real number). The variable 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, , where . If , then the corresponding option is enabled (on). If , then the option is disabled (off). By default, all the options are disabled. At present, only and can be used, whereas the other elements of the option array, , are reserved for future developments.

The option controls the calculation of the relative permeability . If , 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 from USEREOS, i.e., the 5th element of can be an undefined quantity. If , then the default treatment of relative permeability is overridden by the value exported from USEREOS.

The option controls the calculation of the capillary pressure or, to be strict, the relative phase pressure (see Equation (10)). If , 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 from USEREOS, i.e., the 6th element of can be an undefined quantity. If , then the default treatment of capillary pressure is overridden by the value 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 Section 2.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 (), 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 () of every component must be specified. Second, the 3-byte character name (P‐Name_i) for each fluid phase is exported by USEREOS. Third, the 8-byte character name (A-Name_k) and units () for every element of the auxiliary array must be specified. Any A-Name_k must begin with the “#” character. Fourth, the options array, , 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 , , and . It reports the number of phases , the vector , and the integer variable for every phase . The quantity (the “phase identifier”) is aimed at associating the th phase with the phase names (i.e., phase types) specified by GetGlobalParameters. The value of must be equal to the serial number of the P-Name_k, specified by GetGlobalParameters.

A very important argument of PhaseEquilibrium is Mode. If , then the subroutine is supplied by the kernel with the initial guess for the phase equilibrium. The variables np, , and , , are used to specify the initial guess. In the case of , PhaseEquilibrium must calculate the equilibrium using the provided initial guess. The procedure must not change np and . It can only change for a given , , and . Since the number of phases np is not changed, the reported saturations and concentrations and are allowed to take negative or larger-than-one values [19]. If , then the initial guess is not specified. In this case, the subroutine must calculate the phase equilibrium—i.e., np, , and —from scratch. If , then the reported and must satisfy the inequalities and , 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 Section 2.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 Tables 1–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 , , , and in the vector must be kg/m³, J/kg, Pa∙s, and Pa, respectively.

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 Table 5.

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 the ReadConfigurationFile 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 Table 6, where “C-Name_(j)” and “P-Name_i” are the names of the th component and th phase specified by the subroutine GetGlobalParameters. For example, if a component name is “CO₂” and a phase name is “LIQ,” then CCO₂-LIQ is the mass fraction of CO₂ in the phase LIQ, DLIQ is the density of LIQ, and ZFCO₂ is the bulk mass fraction of CO₂.

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: where is the universal gas constant, is the molar concentration, and are the molar enthalpy, and heat capacity at , respectively. For simplicity, we omit the subscript of the only phase of the fluid.

Using the developed module, we present a 1-D study for a three-component gas composed of N₂, O₂, and CH₄. We consider a 1-D flow in a horizontal homogeneous reservoir m, where is the length. The porosity is , and the permeability is mD. Other relevant parameters are the rock density, kg/m³, and the specific rock heat capacity, kJ/(kg∙K). At , the porous medium is saturated with the gas of composition , i.e., it consists of 10% of N₂ and 90% of CH₄. The initial reservoir pressure and temperature are 200 bar and 50°C, respectively. The Dirichlet boundary conditions are imposed at m, i.e., the initial pressure (and temperature) are kept constant at the open boundary m. A gas of different composition , i.e., air, is injected into the reservoir through the boundary m. The injected gas temperature is 30°C at bar. The injection rate is kept constant at 0.2 m/day at reservoir pressure. The permeability 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.

Assuming that all dispersion effects can be neglected, i.e., , the formulated problem admits a simple analytical solution that includes two shocks and (Figure 3). The concentrations are discontinuous at the leading shock S_c. The initial gas composition is preserved ahead of , whereas the porous medium is saturated with the injected gas behind . The reservoir is at the initial temperature °C ahead of the trailing shock and at the injection temperature of 30°C behind . When constructing the analytical solution, one can derive that and 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 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.

(a)

(b)

As shown in Figure 3, 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: where is the isothermal compressibility, is the coefficient of thermal expansion, and are the reference values of pressure and temperature, respectively, and is the specific heat capacity at . According to Equation (18), the component of the mixture forms the phase , where , 2. For example, the phase (component) is water, and the phase is oil.

We consider a 1-D flow in a horizontal homogeneous reservoir at m. The porosity and permeability of the medium are and mD, respectively. The relative permeability is given by

The capillary pressure is assumed to be 0. At , the reservoir is saturated by the oil phase , i.e., , at bar and . The reservoir temperature is kept constant at all . Therefore, its value as well as the values of and does not influence the flow. The initial pressure is kept constant at the open boundary m. The water phase () is injected through the boundary m at constant injection rate of 0.1 m/day.

The compressibilities are assumed to be small such that the phase densities can be considered constant. Therefore, only the phase viscosities, cP and cP, are relevant parameters.

The water injection results in the immiscible displacement of oil from the boundary m to m (Figure 4). The saturation distribution is governed by the Buckley–Leverett solution [22]. This analytical solution consists of the leading displacement front , at which is discontinuous, and the trailing Riemann wave, in which changes continuously, reaching at m. One can show that for Equation (19) and specified phase viscosities and , the shock propagates at a velocity of 0.8 m/day. Thus, is exactly at m at days.

As shown in Figure 4, 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 (Figure 5). 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 , , and the bulk concentrations of the primary chemical elements that the fluid is composed of. For a given , , and , 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 . The chemical speciation is stored in the auxiliary array discussed in Section 2.1.4.

The configuration files (see Section 2.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 m, where . At , the column is saturated with an aqueous solution in equilibrium with calcite at bar and °C. The initial pressure is kept constant at the open boundary m. The column is flushed through the boundary with a MgCl₂ solution at a constant flow rate m/s. The initial and boundary elemental concentrations are listed in Table 8. The flow is assumed to be isothermal at °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.

The mechanical dispersion in this problem is modelled by Equations (7) and (8). Since the study is 1-D, only the longitudinal dispersion length, 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 scales to that in [24] by the factor .

To simulate the fluid flow, we use a fluid consisting of the six components listed in Table 8. 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 MgCl₂ 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.

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.

Data Availability

The source code of the developed USEREOS (v1.0) modules, the compiled executables, and all input data associated with the presented development and benchmark studies can be downloaded at [21, 22]. The source code of GEMS3K can be downloaded at http://gems.web.psi.ch/GEMS3K/. MUFITS executable files are freely available at [21] and http://www.mufits.imec.msu.ru/download.html.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

The authors acknowledge funding from the Russian Science Foundation under Grant # 20-17-00184.

References

S. Geiger, R. Haggerty, J. H. Dilles, M. H. Reed, and S. K. Matthäi, “New insights from reactive transport modelling: the formation of the sericitic vein envelopes during early hydrothermal alteration at Butte, Montana,” Geofluids, vol. 2, no. 3, p. 201, 2002.
View at: Publisher Site | Google Scholar
S. A. Bea, U. K. Mayer, and K. T. B. MacQuarrie, “Reactive transport and thermo-hydro-mechanical coupling in deep sedimentary basins affected by glaciation cycles: model development, verification, and illustrative example,” Geofluids, vol. 16, no. 2, p. 300, 2016.
View at: Publisher Site | Google Scholar
J. Bear, “Porous media,” in Theory and Applications of Transport in Porous Media, pp. 1–98, Springer, 2018.
View at: Google Scholar
J. R. Fanchi, Principles of Applied Reservoir Simulation, Elsevier, 2006.
K. Aziz and A. Settari, Petroleum Reservoir Simulation, Applied Science Publishers, 2002.
K. Pruess, C. Oldenburg, and G. Moridic, “TOUGH2 User’s Guide Version 2.0,” Tech. Rep., Lawrence Berkeley National Lab.(LBNL), Berkeley, CA (United States), 1999.
View at: Google Scholar
T. Xu, N. Spycher, E. Sonnenthal, G. Zhang, L. Zheng, and K. Pruess, “TOUGHREACT Version 2.0: a simulator for subsurface reactive transport under non-isothermal multiphase flow conditions,” Computational Geosciences, vol. 37, no. 6, pp. 763–774, 2011.
View at: Publisher Site | Google Scholar
H. Prommer, D. A. Barry, and C. Zheng, “MODFLOW/MT3DMS-based reactive multicomponent transport modeling,” Ground Water, vol. 41, no. 2, pp. 247–257, 2003.
View at: Publisher Site | Google Scholar
A. Afanasyev, “Hydrodynamic modelling of petroleum reservoirs using simulator MUFITS,” Energy Procedia, vol. 76, pp. 427–435, 2015.
View at: Publisher Site | Google Scholar
A. Afanasyev, T. Kempka, M. Kühn, and O. Melnik, “Validation of the MUFITS reservoir simulator against standard CO₂ storage benchmarks and history-matched models of the Ketzin pilot site,” Energy Procedia, vol. 97, pp. 395–402, 2016.
View at: Publisher Site | Google Scholar
Y. Wang, D. Voskov, M. Khait, and D. Bruhn, “An efficient numerical simulator for geothermal simulation: a benchmark study,” Applied Energy, vol. 264, article 114693, 2020.
View at: Publisher Site | Google Scholar
A. Afanasyev, “Numerical modelling of solute flow dispersion in porous media using simulator MUFITS,” Journal of Physics Conference Series, vol. 1129, article 012002, 2018.
View at: Publisher Site | Google Scholar
K. H. Coats, “An equation of state compositional model,” Society of Petroleum Engineers Journal, vol. 20, no. 5, pp. 363–376, 1980.
View at: Publisher Site | Google Scholar
C. I. Steefel, C. A. J. Appelo, B. Arora et al., “Reactive transport codes for subsurface environmental simulation,” Computational Geosciences, vol. 19, no. 3, pp. 445–478, 2015.
View at: Publisher Site | Google Scholar
M. De Lucia, T. Kempka, and M. Kühn, “A coupling alternative to reactive transport simulations for long-term prediction of chemical reactions in heterogeneous CO₂ storage systems,” Geoscientific Model Development, vol. 8, no. 2, pp. 279–294, 2015.
View at: Publisher Site | Google Scholar
A. Afanasyev and I. Utkin, “Modelling ground displacement and gravity changes with the MUFITS simulator,” Advances in Geosciences, vol. 54, pp. 89–98, 2020.
View at: Publisher Site | Google Scholar
D. A. Kulik, T. Wagner, S. V. Dmytrieva et al., “GEM-Selektor geochemical modeling package: revised algorithm and GEMS3K numerical kernel for coupled simulation codes,” Computational Geosciences, vol. 17, pp. 1–24, 2013.
View at: Publisher Site | Google Scholar
A. Gysi, Y. Mei, and T. Driesner, “Advances in numerical simulations of hydrothermal ore forming processes,” Geofluids, vol. 2020, Article ID 7649713, 4 pages, 2020.
View at: Publisher Site | Google Scholar
H. Salimi, K.-H. Wolf, and J. Bruining, “Negative saturation approach for non-isothermal compositional two-phase flow simulations,” Transport in Porous Media, vol. 91, no. 2, pp. 391–422, 2011.
View at: Publisher Site | Google Scholar
A. Afanasyev and I. Utkin, “USEREOS v1.0,” 2021, https://github.com/utkinis/USEREOS/.
View at: Google Scholar
A. Afanasyev and I. Utkin, “MUFITS. User-supplied EoS-modules,” 2021, http://www.mufits.imec.msu.ru/example-usereos.html/.
View at: Google Scholar
S. E. Buckley and M. C. Leverett, “Mechanism of fluid displacement in sands,” Transactions of AIME, vol. 146, no. 1, pp. 107–116, 1942.
View at: Publisher Site | Google Scholar
A. Yapparova, T. Gabellone, F. Whitaker, D. A. Kulik, and S. K. Matthäi, “Reactive transport modelling of dolomitisation using the new CSMP++GEM coupled code: governing equations, solution method and benchmarking results,” Transport in Porous Media, vol. 117, no. 3, pp. 385–413, 2017.
View at: Publisher Site | Google Scholar
H. Shao, S. V. Dmytrieva, O. Kolditz, D. A. Kulik, W. Pfingsten, and G. Kosakowski, “Modeling reactive transport in non-ideal aqueous-solid solution system,” Applied Geochemistry, vol. 24, no. 7, pp. 1287–1300, 2009.
View at: Publisher Site | Google Scholar
V. J. Azad, C. Li, C. Verba, J. H. Ideker, and O. B. Isgor, “A COMSOL-GEMS interface for modeling coupled reactive-transport geochemical processes,” Computational Geosciences, vol. 92, pp. 79–89, 2016.
View at: Publisher Site | Google Scholar
L. H. Damiani, G. Kosakowski, M. A. Glaus, and S. V. Churakov, “A framework for reactive transport modeling using FEniCS–Reaktoro: governing equations and benchmarking results,” Computational Geosciences, vol. 24, no. 3, pp. 1071–1085, 2020.
View at: Publisher Site | Google Scholar
D. L. Parkhurst and C. A. J. Appelo, PHREEQC (Version 3.0.4)--a computer program for speciation, batch speciation, one-dimensional transport, and inverse geochemical calculations, U.S. Geological Survey Techniques and Methods, book, 2013.
J. Jatnieks, M. De Lucia, D. Dransch, and M. Sips, “Data-driven surrogate model approach for improving the performance of reactive transport simulations,” Energy Procedia, vol. 97, pp. 447–453, 2016.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2021 Andrey Afanasyev and Ivan Utkin. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

221

Downloads

674

Citations