# 1. Introduction
GPSFLOW is a General Purpose Subsurface Flow Simulator. It solves the governing equations for non-isothermal, multiphase, multicomponent flows in porous or fractured media. The code is designed to run on hybrid massively parallel computing architectures as well as workstations and PC computers which are equipped with multicore CPUs and GPUs. Parallelization is achieved through domain decomposition and multiple thread programming with the support of OpenCL or CUDA for running on GPUs. The codes are developed in C++ with Object-oriented programming (OOP) technology.
The GPSFLOW simulation capabilities are very diverse. Its main application areas cover CO2 geological sequestration, geothermal reservoir engineering, nuclear waste isolation studies, environmental assessment and remediation, gas/energy subsurface storage, unconditional oil/gas production, and flow and transport in variably saturated media and aquifers.
The simulator fully accounts for movement of different mobile phases, solid phase, heat flow, and the transition of components between the available phases, which may appear and disappear depending on the changing thermodynamic state of the system. The flowing processes under pressure, viscous, capillary and gravity forces can be described as Darcy flow using the multiphase extension of Darcy's law, non-Darcy flow using the Forchheimer Equation, or non-Newtonian flow and treated as power law or Bingham fluid. In addition, diffusive mass transport is allowed to occur in all mobile phases, decay/biodegradation and adsorption are accounted for mass components. Modeling capabilities of the simulator is under expanding and improving. Users may refer to the most updated simulator development technical report for details.
GPSFLOW uses the integral finite difference method (IFDM) (Edwards, 1972; Narasimhan and Witherspoon, 1976) for space discretization, which is similar to the approaches adopted in TOUGH family codes (Pruess et al. 1999). The IFDM allows using regular or irregular model domain discretizations in one, two, and three dimensions. Time discretization can be fully implicitly as a first-order backward finite difference or uses adaptive implicit method (AIM) (ref.). The resulting nonlinear algebraic equations are solved using Newton-Raphson iterations for each time step, which involves the calculation of a Jacobian matrix and the solution of a set of linear equations.
# 2. Compilation, Installation and Running Simulations
# 2.1 Software compilation requirement
Most users are not expected to compile the codes by themselves. If you get the executable and do not interested in compiling the codes, you may skip this section and directly go to the section for running the simulator.
The simulator can be compiled under different platforms. The following software packages are need for different platform and options:
· C++ compilers, compilation for some libraries may need FORTRAN compiler.
· CMAKE (for Linux, Cygwin or MacOS only).
· MPI, METIS/ParMETIS (for distributed memory parallel computing only)
· VISUAL Studio 2015 or later (for Windows only)
· OpenCL or CUDA (for using GPUs only)
· AMGCL linear solver library (default) (Demidov, 2018).
· Additional linear solvers (the interface of linear solvers for following libraries have been implemented in the simulator, user can use them if they are compiled with.)
-- PETSc
-- Trilinos
-- FASP
-- AMGX (for multiple GPUs simulations, CUDA)
-- rocALUTION (for multiple GPUs simulations, OpenCL)
-- ViennaCL (for single GPU simulations)
· GPSFLOW
# 2.2 Compiling GPSFLOW
Current version GPSFLOW does not require specific installation process. User only needs to copy the executable to a folder where it can be easily accessed.
For running on Windows, some DLL files may be required, if they are not existed in your computer. These DLL files can be found either in the Visual Studio, compilers, or Windows. It is also possible to find them on online and download them to your computer.
For MPI parallel simulations, it is required that the MPI has been installed in the computer where the software will be run.
# 2.3 Installation
Current version GPSFLOW does not require specific installation process. User only needs to copy the executable to a folder where it can be easily accessed.
For running on Windows, some DLL files may be required, if they are not existed in your computer. These DLL files can be found either in the Visual Studio, compilers, or Windows. It is also possible to find them on online and download them to your computer.
For MPI parallel simulations, it is required that the MPI has been installed in the computer where the software will be run.
# 2.4 Running GPSFLOW
GPSFLOW has four running modes, namely stand-alone running, OPENMP multithreaded parallel running, MPI parallel running, and hybrid parallel running.
User can choose the best way to carry out their simulations based on the model size and the computer structure that will be used for the simulations. If linear equation systems for the simulation are to be solved using GPU, the user must select the linear equation solver "ViennaCL" (suitable for single GPU calculation), "G_AMGX" or "rocAlution" (for single or multi-GPU on different GPU platform) under the keyword SOLVER_ in the main input file (INFILE). Regardless of the mode of running, the linear equation will be solved by GPU if the solver uses ViennaCL, AMGX, or rocAlution. If you choose a solver at runtime, make sure its corresponding library was included when compiling the source code.
\1) Stand-alone running
Generally, it is used for small-scale model, i.e. single-threaded simulation calculation for testing or debugging. The following instruction will start the stand-alone running:
./Gpsflow-t 1
If the solver selection is ViennaCL, AMGX, or RocAlution, the linear equation will be solved by a GPU accelerator.
\2) OPENMP multithreaded parallel running
GPSFLOW allows multithreaded computation, which has a significant acceleration effect on models of a larger size. The number of computing cores for the default operation is 8 (8 threads). When the running computer has less than 8 computing cores, all of its cores will be used:
./gpsflow
The above instruction will run the simulation with 8 cores or all the cores that the computing node has, if the available core number is less than 8.
./gpsflow-t 16
The above instruction will run the simulation with 16 computing cores
If the solver selection is ViennaCL, AMGX, or RocAlution, the linear equation will be solved by a GPU accelerator.
\3) MPI parallel running
MPI parallel model simulation is performed by domain decomposition, which requires the installation of MPI and METIS/ ParMETIS. In most cases, MPI parallel running is more efficient than OPENMP multithreaded parallel running, especially for more complex model computations. In general, MPI parallel simulation is run on for distributed memory computers, and it takes up more memory than OPENMP multithreaded parallel running. Depending on the size of the model, users can choose different amounts of CPU/ core for calculation by:
mpirun –np m ./gpsflow
where m is the number of CPUs/cores requested for the calculation.
If the solver selection is ViennaCL, AMGX, or RocAlution, the linear equation will be solved by one or multiple GPU accelerators. Requesting the number of GPU for the simulation is by the computer job scheduler. Requested GPU number must be equal or less than m.
\4) Hybrid running
At present, most supercomputers have hybrid memory structures. Usually a supercomputer is consist of computing nodes. Each node has a number of CPUs and GPUs, and each CPU has several computing cores. All computing cores on a node share the memory in the node. For this kind of computer structure, the hybrid parallel computing is able to achieve good efficiency. GPSFLOW can perform hybrid parallel computing using the following commands:
mpirun –np m ./gpsflow –t n
where m is the number of computing unit (CPU, node or other) participating in the computation, and n is the number of threads/cores participating in the computation within each computing unit.
If the solver selection is ViennaCL, AMGX, or RocAlution, the linear equation will be solved by one or multiple GPU accelerators. Requesting the number of GPU for the simulation is by the computer job scheduler. Requested GPU number must be equal or less than m.
The basic idea of hybrid parallel computing is to divide the large-scale model computation into m subdomains (the number of computing unit used) through domain decomposition, and each computing unit is responsible for one of them. The computation tasks within each unit is completed jointly by multiple cores (n) and/or GPUs.
When running GPSFLOW, the command line allows you to specify some parameters, and the following parameters are optional in the current version:
-t m
Specify the number of threads participating in the calculation (m);
-q
Fast computation mode. The software automatically selects the coarse convergence criteria of linear and nonlinear iterations, and selects the most suitable pre-condition method and the linear equation solution method to achieve the goal of rapid calculation. This method is not stable at present, and sometimes there is no convergence phenomenon.
-f fileN
This parameter specifies the name of the model main input file. The main input file is “fileN.in” and the output file is “fileN.out”. If this file name is not specified, the program uses the default input and output file names, i.e. INFILE and output.dat.
-i
The model will be run in isothermal condition, no matter what condition is requested in the main input file.
-v m
Volume weighted average of the porosity and permeability for the grids near wells (for solving the problem of poor convergence due to strong heterogeneity). *m* is the level to be averaged (starting from the well grid), i.e. the parameters for grids within *m*-layer connections are averaged. If *m* is a negative number, the porosity and permeability of these grids will be replaced by the corresponding parameters of the rock types of these grids, i.e. the heterogeneous parameters inputted via svparam.dat or INCON file will be ignored.
-w m
Volume weighted average of the porosity and permeability of each grid in the model (for solving the problem of poor convergence due to strong heterogeneity). *m* is the level to be averaged (starting from current grid), i.e. the grids within the m-layer connections are averaged with volume weighted.
-n m
*m* is the number of nodes involved in the calculation. It only works if you use the rocAlution and AMGX libraries and need multi-node computation. At present, OpenCL can only find all GPUs on one node. In order to perform GPU calculation with more than one node, it is required the input of the number of nodes involved in the calculation. Through -*n* input, the total number of GPUs involved in the calculation is the number of nodes multiplied by the number of GPUs on each node. For example:
mpirun –np 8 ./gpsflow –t 8 –n 2
The above instruction requests 8 MPI tasks and 8 threads for each task for this simulation run. The number of cores involved in the calculation is 8*8. There is a total of 2 nodes. If each node has 4 GPUs, a total of 8 GPUs participates in the calculation.
# 3. Preparation for Input Files
GPSFLOW input is through one or several input files. These files may include main input file (default name INFILE), mesh file (MESH), file for spatial distribution of rock property parameters (svparams.dat), initial condition file (INCON), source/sink or well information file (GENER), and dynamic simulation run control file (Run_Parameters.dat). Except the main input file, all other files can be either included in main input file or optional.
The main input of GPSFLOW is through keywords. The keywords are divided into two groups, the main keywords and secondary keywords. A main keyword may have several secondary keywords. Each secondary keyword must belong to a main keyword. Text lines or any locations with // are treated as the beginning of comments. A line does not associated with any keyword will be neglected by the program. Table 3-1 lists the main keywords.
Table 3-1 Main keywords
| Keywords | Description |
|---|---|
| TITLE_ | The first line of the main input file always provides the title of the model |
| ROCKS_ | Define rock properties |
| INCON_ | Input initial conditions |
| PARAM_ | Input Computation parameters |
| GENER_ | Source/sink or well |
| ELEME_ | Define gridblock geometry |
| CONNE_ | Define connection between neighboring gridblocks |
| TIMES_ | Specify output time points |
| FLUID_ | Define fluid properties |
| PVTDA_ | Input PVT table |
| SOLVE_ | Provide linear solver parameters |
| ENDCY | End of the input, the last effect main keyword |
| MODEL_ | Provide model definition parameters |
| FOFTS_ | Output time series for gridblocks |
| COFTS_ | Output time series for connections |
| GOFTS_ | Output time series for wells or source/sinks |
| ENDFI_ | End the input without run the simulation. |
| RPCAP_ | Provide default relative permeability and capillary pressure functions |
| TIMBC_ | Define time dependent boundary condition |
| INDOM_ | Input domain specified initial conditions |
| DIFFU_ | Input parameters for diffusion flux calculations |
| INGEQ_ | Input parameters for the gravity equilibrium calculation |
| ADVAN_ | |
| SVPAR_ | Specify porosity and permeability spatial variation condition |
| TRCON_ | Define changes of connection areas between two gridblocks |
| OUTCT_ | Output control |
| PREME_ | Mesh preprocessing |
| INICO_ | Initial condition, same as INCON_ |
| WELLS_ | Same as GENER_ |
| COMPO_ | Input property parameters for model mass components. |
| GMESH_ | Mesh generation |
| GRIDB_ | Same as GMESH_ |
| PVT_D_ | Same as PVTDA_ |
# 3.1 Conventions and Notation
Main keywords always have an underscore at the end and secondary keywords are not a real English word, this is for avoiding potential occurrence a keyword in the comments or prompting statements. Both uppercase and lowercase are treated as the same for keywords. The main keywords do not request a value input for itself, but introduces a section for input of a group parameters through secondary keywords. Each secondary keyword requests input of one or several parameters.
In general, SI units are required for all input parameters, unless otherwise specified. Table 3-2 lists the most common units used in the codes.
Table 3-2 Units
| Quantity | Units |
|---|---|
| Pressure | Pascal (Pa) |
| Temperature | Celcius (oC) |
| Distance | Meter (m) |
| Volume | Cubic Meter (m3) |
| Time | Second (s) |
| Velocity | Meter/Second (m/s) |
| Density | Kilogram/Cubic Meter (kg/m3) |
| Enthalpy | Joule/Kilogram (J/kg) |
| Mass | Kilogram (kg) |
| Rate | Kilogram/Second (kg/s) or Cubic Meter/Day (m3/d) |
| permeability | Square Meter (m2) |
| Viscosity | pascal-second (Pa.s) |
| Concentration | molarity (m) |
The model inputs are through text lines. Each text line may contain prompting statement, keyword, corresponding values, comment message, and separating symbol (😃. Figure 3-1 shows a typical input format. The prompting statements and comments are optional. The input value must follow the separating symbol (😃. Spaces are allowed between any parts of the inputs. Comments are always started with”//”. In the figure 3-1, only the keywords (ROCK_, RName, Rdens), values (sand, 2600.0) and separating symbol ”:” must be presented. All others are optional. Many inputs are optional. If no input is provided for a parameter, a default value will be adopted. Unless specified otherwise the sequence for parameter input can be random.
Figure 3-1 Typical input format for GPSFLOW simulation
# 3.2 ROCKS_, input for rock properties
Specifies material properties for different reservoir domains. You can define as more domains as you need.
Format: description (keyword): values. The description is optional, which can be specified by user. The first line for each rock must be the medium name.
Keywords:
RName: <string>
Rock name, the length of the name must be equal or less than 5 characters. This keyword always be the first line of inputs for each rock.
RDens: <float> [kg/m3] （default is 2650）
Rock grain density.
Rporos: <float> （default is 0.2）
Rock porosity for all elements belonging to this rock for which porosity has not been specified in the initial condition inputs or in the input for spatial variable parameters.
permInX: <float> [m2] (default is 1.0e-12）
Absolute permeability along X axis for all elements belonging to this rock for which absolute permeability in X direction has not been specified in the initial condition inputs or in the inputs for spatial variable parameters. User may be able to input the permeability values in darcy or millidarcy, if the comments in main keyword “ROCKS_” input line contains "KUNIT=D" or "KUNIT=M" (this also applies to the input for permInY and permInZ). If the permeabilities in X, Y, Z direction are the same, only one of them needs to be inputted.
permInY: <float> [m2] (default is 1.0e-12）
Absolute permeability along Y axis for all elements belonging to this rock for which absolute permeability in Y direction has not been specified in the initial condition inputs or in the input for spatial variable parameters.
permInZ: <float> [m2] (default is 1.0e-12）
Absolute permeability along Z axis for all elements belonging to this rock for which absolute permeability in Z direction has not been specified in the initial condition inputs or in the input for spatial variable parameters.
KThrW: <float> [W/m°C]
Rock heat conductivity under fully liquid-saturated conditions.
SpcHt: <float> [J/kg°C]
Rock grain specific heat.
Cbeta: <float> (default is 0.0, treated as Darcy’s flow)
Beta Coefficient for Non-Darcy Flow.
RCompress: <float> [1/Pa] (default is 0.0)
Rock compressibility.
Expansivity: <float> [1/°C]
Rock Expansivity.
KthrD: <float> [W/m2] (default is the value inputted through KThrW)
Rock heat conductivity under desaturated conditions
Tortu: <float> (default is 0.0)
Tortuosity factor for binary diffusion. If it is 0.0, a porosity and saturation-dependent tortuosity will be calculated internally from the Millington and Quirk (1961) model.
Klink: <float> [1/Pa]
Klinkenberg parameter b (Klinkenberg, 1941) for enhancing gas phase permeability according to the relationship kgas = kliq * (1 + b/P).
CritSat:<float>
Critical Mobile Saturation.
permExpon:<float>
Permeability reduction exponent for solid phases.
Gbeta:<float>
Empirical parameter for geomechanical porosity reduction.
Gama:<float>
Empirical parameter for geomechanical permeability reduction.
RelPermEqNum: <int>
The function number chosen from a list of relative permeability functions implemented in the simulator. Available functions and required parameters are discussed at the late part of this section.
RP parameters: <float>…
Relative permeability function parameters for the selected function. The number of parameters is different for different functions. The input parameters need to be separated by a comma “,”.
PcapEqNum: <int>
The function number chosen from a list of capillary pressure functions implemented in the simulator. Available functions and required parameters are discussed at the late part of this section.
CP_parameters: <float>…
Capillary pressure function parameters for the selected function. The number of parameters is different for different functions. The inputted parameters need to be separated by a comma “,”.
S_vug_ini: <>
The starting saturation in the case of non-zero water phase percolation when the water phase flows out of the large cave unit.
S_vug_end:<>
When the aqueous phase of the large cave unit flows out, the saturation of the aqueous phase increases to 1.0.
Sw_Kw_Ko: <float>…
Reading 4 parameters in the sequence of water saturation (sw) in an oil-water two phase system, relative permeability of aqueous phase (krw) at the corresponding sw, relative permeability of oil phase (kro) at the corresponding sw, and oil-water capillary pressure at the corresponding sw. The inputted parameters need to be separated by a comma “,”. Defaults are 0.0. Repeat this keyword for a multiple line table inputs. This keyword input is required by RelPermEqNum=3 and PcapEqNum=3.
Sw_Kw_Ko: <float>…
Reading 5 parameters in the sequence of gas saturation (sg) in an oil-gas two phase system, relative permeability of gas phase (krw) at the corresponding sg, relative permeability of oil phase (kro) at the corresponding sg, oil-gas capillary pressure at the corresponding sg, and water-gas capillary pressure at gas saturation=sg in a gas-water two phase system. The inputted parameters need to be separated by a comma “,”. Defaults are 0.0. Repeat this keyword for a multiple line table inputs. This keyword input is required by RelPermEqNum=3 and PcapEqNum=3.
FracMAb: <float> (default is 0.0)
Fraction of organic carbon present in the rock; used for calculating amount of VOC adsorbed. It will take effect only in the simulation where the absorption of VOC is considered.
BasedOn: <int>
This keyword is for simplifying the input of rock parameters. Parameters of current rock are assigned by copying all values of parameters from previous defined rock. Any followed inputs will replace the assigned values. The BasedOn must be larger than 0 and less than the number of already defined rocks.
Relative Permeability Functions:
The relative permeabilities are assumed to be functions of fluid saturations only and not to be affected by non-Newtonian behavior when simulating non-Newtonian fluid flow. The relative permeability to water phase is described.
where
In addition, two functional forms of three-phase relative permeability are also provided as alternatives. The Brooks-Corey type of functions may be extended into three-phase flow as,
And
Where
And,
with
The extended three-phase relative permeability of the van Genuchten model is given by (Kaluarachchi and Parker, 1989),
And,
RelPermEqNum = 1 linear functions.
krl increases linearly from 0 to 1 for aqueous saturation in the range RP_parameters1 and RP_parameters2.
krg increases linearly from 0 to 1 for gas saturation in the range RP_parameters3 and RP_parameters4.
krg increases linearly from 0 to 1 for oil/aqu2 saturation in the range RP_parameters5 and RP_parameters6.
Default values for RP_parameters1, range RP_parameters3 and range RP_parameters5 are 0.0, for RP_parameters2, range RP_parameters4 and range RP_parameters6 are 1.0.
（1）
RelPermEqNum = 2
Three-phase functions of Parker et al. (1987).
where
Four parameters are required:
RP_parameters1, residual gas saturation
RP_parameters2, residual water saturation
RP_parameters3, residual oil saturation
RP_parameters4, parameter
(2)
RelPermEqNum = 3
Tabulated two-phase relative permeability, together with Stone II function for three-phase relative permeability.
Four parameters are required:
RP_parameters1, residual gas saturation.
RP_parameters2, residual water saturation.
RP_parameters3, residual oil saturation.
RP_parameters4, oil phase relative permeability (
For RelPermEqNum = 3, additional inputs by reading tabulated data for two-phase relative permeabilities and capillary pressure are required through keyword “Sw_Kw_Ko” and “Sg_Kg_Ko”.
(3)
RelPermEqNum = 4
Relative permeability function for large cave or vug.
This function is for caves or vugs which are supposed to be multiphase gravity separation and instant equilibrium. In the model gridblocks with this type of rock, the multiphase fluids are distributed from top to bottom in the sequence of gas, oil and water.
Two parameters are required:
RP_parameters1, minimum water saturation for relative permeability of aqueous phase is non zero. At this saturation water flux from cave/vug to surrounding appears.
RP_parameters2, initial water saturation for relative permeability of aqueous phase is 1.0.
(4)
(4) RelPermEqNum = 5
Grant's curves (Grant, 1977).
where
with
Restrictions:
Two parameters are required:
RP_parameters1, residual gas saturation
RP_parameters2, residual water saturation
RP_parameters1+ RP_parameters2 must be less than 1.0. If this function applies to three-phase flow, the oil phase relative permeability equals to
RelPermEqNum = 5 is different from RelPermEqNum = 10 by simplified calculation for
(5)
(5) RelPermEqNum = 6
three-phase Brook-Corey function.
A modified version of the Brooks-Corey model (Luckner et al., 1989) has been implemented to prevent the capillary pressure from decreasing towards negative infinity as the effective saturation approaches zero. The modified Brooks-Corey model is invoked by setting both IRP and ICP to 10. $$ k_{rl}=S_{ek}\frac{2+3\lambda}{\lambda} $$
$$ k_{rg}=\left{\begin{matrix}(1-S_{ek})^2(1-S_{ek}\frac{\frac{2+\lambda}{\lambda}}{\lambda}\quad if\quad RP(3)=0\1-k_{rl}\qquad \qquad \qquad \quad \ \ if\quad RP(3)\neq 0\end{matrix}\right. $$
Where $$ S_{ek}=\frac{S_l-S_{lrk}}{1-S_{lrk}-S_{gr}} $$ Four parameters are required:
RP_parameters1, residual gas saturation.
RP_parameters2, residual water saturation.
RP_parameters3, residual oil saturation.
RP_parameters4, the exponential constant
(6)
RelPermEqNum = 7
two-phase van Genuchten function (van Genuchten, 1980).
$$
k_{rl}=\left{\begin{matrix}\sqrt{S^}\left{\begin{matrix}(1-[S^]^{1/\lambda)^\lambda}\end{matrix}\right}^2\quad ifS_1<S_{1s}\\qquad \quad 1\qquad \qquad \quad \quad \ \ ifS_1\geqslant S_{1s}\end{matrix}\right.
$$
Gas relative permeability can be chosen as one of the following two forms, the second of which is due to Corey (1954)
$$
k_{rg}=\left{\begin{matrix}1-k_{rl}\qquad \qquad \quad \ ifS_{gr}=0\(1-\hat{S})^2(1-\hat{S}^2)\quad ifS_{gr}>0\end{matrix}\right.
$$
subject to the restriction
Here,
Four parameters are required:
RP_parameters1, van Genuchten parameter
RP_parameters2, residual water saturation.
RP_parameters3, saturated water saturation.
RP_parameters4, residual gas/oil saturation.
This function can be used for gas-water or oil-water two-phase flow.
(7)
(7) RelPermEqNum = 8
two-phase relative function of Pickens et al. (1979)
One parameter is required:
RP_parameters1, the exponential constant.
(8)
RelPermEqNum = 9
all phases perfectly mobile
No parameter is required. Relative permeability for all phases equals to 1.0.
(9)
RelPermEqNum = 10
Grant's curves (Grant, 1977). $$ k_{rl}=\hat{S}^4\k_{rg}1-k_{rl} $$ where $$ \hat{S}=(S_l-S_{lr})/(1-S_{lr}-S{gr}) $$ With $$ S_{lr}=RP(1);S_{gr}=RP(2) $$ Restrictions: $$ RP(1)+RP(2)<1 $$ Two parameters are required:
RP_parameters1, residual gas saturation.
RP_parameters2, residual water saturation.
RP_parameters1+ RP_parameters2 must be less than 1.0. If this function applies to three-phase flow, the oil phase relative permeability equals to
(10)
RelPermEqNum = 11, 12, 13
three-phase Stone’s method modified by Aziz (1979).
Four parameters are required:
RP_parameters1, residual water saturation.
RP_parameters2, residual oil saturation.
RP_parameters3, residual gas saturation.
RP_parameters4, the exponential constant
If RelPermEqNum = 12, a correction factor is applied to
If RelPermEqNum = 13, the original
(11)
RelPermEqNum = -1, -2, -3
water-oil two-phase relative permeability functions for Buckley Leverett solution.
Capillary Pressure Functions:
(1)
PcapEqNum = 1
Capillary pressure functions are generally to be supplied through table lookup, and in addition one functional form of capillary pressures from Parker et al. (1987) is included with modification for use. In this modified capillary function,
$$
P_{cow}^{ow}=\frac{\rho_wg}{\alpha_{ow}}\left{\begin{matrix}(\overline{S}w)^{-1/\gamma}-1\end{matrix}\right}^{1/\beta} \
P{cgo}^{go}=\frac{\rho_wg}{\alpha_{go}}\left{\begin{matrix}(1-\overline{S}g)^{-1/\gamma}-1\end{matrix}\right}^{1/\beta}\
P{cgw}^{gw}=\frac{\rho_wg}{\alpha_{go}}\left{\begin{matrix}(1-\overline{S}g)^{-1/\gamma}-1\end{matrix}\right}^{1/\beta}
$$
where
Seven parameters are required:
CP_Parameters1, residual water saturation.
CP_Parameters2, van Genuchten parameter,
CP_Parameters3, maximum capillary pressure (Pa).
CP_Parameters4, the oil threshold saturation above which gas-water interfaces cease to exist.
CP_Parameters5,
CP_Parameters6,
CP_Parameters7,
(2)
PcapEqNum = 2
three-phase capillary functions from Parker et al. (1987) $$ m=1-1/n\ \overline{S}w=(S_w-S_m)/(1-S_m)\ \overline{S}I=(S_w+S_n-S_m)/(1-S_m)\ P{cgn}=-\frac{\rho_wg}{\alpha{gn}}[(\overline{S}1)^{-1/m}-1]^{1/n}\ P{cgw}=-\frac{\rho_wg}{\alpha_{nw}}[(\overline{S}w)^{\frac{-1}{m}}-1]^{\frac{1}{n}}-\frac{\rho_wg}{\alpha{gn}}[(\overline{S}1)^{\frac{-1}{m}}-1]^{\frac{1}{n}} $$ where $$ S_m=CP(1); n=CP(2); \alpha{gn}=CP(3);\alpha_{nw}=CP(4) $$ Four parameters are required:
CP_Parameters1, residual water saturation.
CP_Parameters2, van Genuchten parameter,
CP_Parameters3,
CP_Parameters4,
These functions have been modified so that the capillary pressures remain finite at low aqueous saturations using linear interpolation.
(3)
PcapEqNum = 3
table look up for capillary pressures
In general, PcapEqNum = 3 is used together with RelPermEqNum = 3. The input requires one optional parameter:
CP_Parameters1, the oil threshold saturation above which air-water interfaces cease to exist. This is an optional parameter in current version.
For PcapEqNum = 3, additional inputs by reading tabulated data for two-phase relative permeabilities and capillary pressure are required through keyword “Sw_Kw_Ko” and “Sg_Kg_Ko”. This table input is shared with the relative permeability calculation with RelPermEqNum = 3.
(4)
PcapEqNum = 4
linear function $$ P_{cap}=\left{\begin{matrix}-CP(1)\quad \ for\ S_1\leqslant CP(2)\0\qquad \qquad \ for\ S_1\leqslant CP(2)\-CP(1)\frac{CP(3)-S_1}{CP(3)-CP(2)}\ \ for\ CP(2)<S_1<CP(3)\end{matrix}\right. $$ Three parameters are required:
CP_Parameters1, the maximum capillary pressure.
CP_Parameters2, initial water saturation from which capillary pressure starts to reduce.
CP_Parameters3, ending water saturation at which capillary pressure reduces to 0.
If water saturation is between CP_Parameters2 and CP_Parameters3, its corresponding capillary pressure is linearly interpolated from the range of 0 to CP_Parameters1.
(5)
PcapEqNum = 5
function of Pickens et al. (1979)
P_{cap}=-P_0\left\{\begin{matrix}In[\frac{A}{B}(1+\sqrt{1-\frac{B^2}{A^2}})]\end{matrix}\right\}^\frac{1}{X}
withCP_Parameters1, the maximum capillary pressure.
CP_Parameters2, residual water saturation.
CP_Parameters3,
CP_Parameters4, the exponential constant x.
（6）
PcapEqNum = 6
Milly’s function (Milly, 1982) $$ P_{cap}=-97.783\times10^A $$ with $$ A=2.26\left(\frac{0.371}{S_l-S{lr}}\right)^{1/4} $$ where $$ S_{lr}=CP(1) $$ Restriction $$ CP(1)\geqslant 0 $$ One parameter is required:
CP_Parameters1, residual water saturation.
(7)
PcapEqNum = 7
van Genuchten function (van Genuchten, 1980) $$ P_{cap}=-P_0\left([S^]^{-1/\lambda}-a\right)^{1-\lambda} $$ subject to the restricition $$ -P_{max}\leqslant P_{cap}\leqslant 0 $$ Here, $$ S^=(S_1-S_{lr})/(S_{ls}-S_{lr}) $$ Parameters: $$ CP(1)=\lambda=1-1/n \CP(2)=S_{lr}(\rm should\ be\ chosen\ smaller\ than\ the\ corresponding\ \\rm parameter\ in\ the\ relative\ permeability\ function; see\ note\ blow.)\ CP(3)=1/P_0=\alpha/\rho_{w}g\ CP(4)=P_{max}\ CP(5)=S_{ls} $$ Five parameters are required:
CP_Parameters1, van Genuchten parameter,
CP_Parameters2, residual water saturation.
CP_Parameters3,
CP_Parameters4, the maximum capillary pressure, P_{max}
CP_Parameters5,
(8)
PcapEqNum = 0
No capillary pressure
No input parameter is required. Capillary pressure is 0 for all saturations.
Example:
ROCKS_
The first medium name (RName): rock1
Rock density (RDens): 2600.0
Rock Porosity (Rporos): 0.1
Permeability in X (permInX): 9.869e-15
Permeability in Y (permInY): 9.869e-15
Permeability in Z (permInZ): 9.869e-15
Heat conductivity,liquid-saturated (KThrW):
Rock grain specific heat (SpcHt):
Tortuosity factor (Tortu):
Medium pore compressibility (RCompress): 1.00E-11
Relative permeability function number (RelPermEqNum): 3
Relative permeability function parameters (RP_parameters): 0.02, 0.23, 0.3, 0.9
Capillary pressure function number (PcapEqNum): 3
Capillary pressure function parameters (CP_parameters): 0.1
Table input of RelP and PC for water-oil system (Sw_Krw_Kro): 0.23000, 0.00000, 1.00000
Sw_Kw_Ko: 0.25000, 0.00300, 0.90000
Sw_Kw_Ko: 0.30000, 0.01200, 0.70000
Sw_Kw_Ko: 0.35000, 0.02600, 0.52000
Sw_Kw_Ko: 0.40000, 0.04500, 0.37000
Sw_Kw_Ko: 0.45000, 0.07600, 0.25000
Sw_Kw_Ko: 0.50000, 0.11500, 0.16000
Sw_Kw_Ko: 0.55000, 0.16500, 0.08200
Sw_Kw_Ko: 0.60000, 0.22000, 0.02000
Sw_Kw_Ko: 0.62000, 0.24800, 0.00500
Sw_Kw_Ko: 0.63000, 0.26000, 0.00100
Sw_Kw_Ko: 0.64000, 0.27000, 0.00000
Table input of RelP and PC for oil-gas system (Sg_Krg_Kro): 0.0, 0.0,1.000
Sg_Kg_Ko: 0.020, 0.000, 0.860
Sg_Kg_Ko: 0.030, 0.001, 0.800
Sg_Kg_Ko: 0.050, 0.002, 0.690
Sg_Kg_Ko: 0.075, 0.005, 0.580
Sg_Kg_Ko: 0.100, 0.010, 0.480
Sg_Kg_Ko: 0.125, 0.017, 0.400
Sg_Kg_Ko: 0.150, 0.026, 0.330
Sg_Kg_Ko: 0.200, 0.048, 0.230
Sg_Kg_Ko: 0.300, 0.106, 0.140
Sg_Kg_Ko: 0.400, 0.200, 0.040
# 3.3 PARAM_, input for computational parameters
This keyword is for input of computational parameters, time stepping information, default initial conditions, and other computation related parameters.
Keywords:
MaxNI: <int> (default is 8)
Maximum number of Newtonian iterations per time step.
**MaxTSN**: <int> (default is 1000000)
Maximum number of time steps for current simulation run.
**OutputN:** <int> (optional, default is 1)
Output amount control.
**CPUTime:**<float>
Maximum CPU time requested for this simulation.
OutputFQ**: <int> (default is 10000)
Output frequency, write output at every OutputFQ time steps.
**NRIterationInfo:**<int> (optional, default is 0)
Print Newton Iteration Information.
**TSInfo<int>**
Time Stepping Information.
**JacobianInfo**<int> (optional, default is 0)
Jacobian Matrix Information.
**WellInfo** <int> (optional, default is 0)
Well Information
**SolidEF<int>** (optional, default is 0)
Selection of Solid Phase Effects.
**EOSInfo <int>** (optional, default is 0)
EOS Information.
**SolverInfo** <int> (optional, default is 0)
Solver Information.
**UWeighting**: <int> (default is 10)
Weighting methods for evaluation of mobility and permeability at interfaces.
UWeighting=
10: harmonic weighting for permeability, upstream weighting for mobility.
11: harmonic weighting for both permeability and mobility.
12: harmonic weighting for permeability, arithmetic average for mobility.
0: upstream weighting for both permeability and mobility.
1: upstream weighting for permeability, harmonic weighting for mobility.
2: upstream weighting for permeability, arithmetic average for mobility.
**Option_FluidComposition**: <int> (default is 0)
**Option_ThermalConductivity**: <int> (default is 0)
Chooses the interpolation formula for heat conductivity as a function of liquid saturation (Sl)
Option_ThermalConductivity=
0: C(S1)=CDRY+SQRT(S1*[CWET-CDRY])
1: C(S1)=CDRY+S1*(CWET-CDRY)
2: same as 0, but neglect the contribution from gas phase.
Where Cdry is the rock heat conductivity under desaturated conditions. The gas contribution for thermal conductivity is:
C(Sg)=porSgCg
C=C(Sl)+C(Sg)
Option_Aim:** <int> (default is 1)
Define the non-linear equation solution method, Option Aim=
0: FIM
1: AIM
2: IMPES
**SSRateInterpolation:** <int> (optional, default is 2)
Selection of Interpolation Procedure for Time-Dependent Sink/source Data.
**PermeabilityAdjustment:** <int> (optional, default is 0)
Selection of Permeability Adjustment Approach.
**Option_GasSolubility:** <int> (optional, default is 0)
Selection of Gas Solubility Method.
**DoublingDtCriterion:** <int> (optional, default is 4)
Criterion for Time-Step Size Increasing.
**Option_InterfaceDensity:** <int> (optional, default is 0)
Interface Fluid Density Calculation Method.
**BaseDiffusionCoef:** <float> [m2]
The base gas diffusion coefficient.
**DiffusionExpon:** <float>
Temperature dependence of gas phase diffusion coefficient.
**DiffusionStrength:** <float>
Temperature-related exponent of base gas diffusion coefficient**.**
**AIM_SwitchM:** <int> (default is 0)
Choose method for FIM and IMPES switching when AIM is selected.
0: saturation change. If saturation change is larger than given criterion, FIM will be used.
1: Courant–Friedrichs–Lewy (CFL) condition. If C<1.0 using IMPES; if C>1 using FIM.
**AimFullyUpdate:** <int> (default is 0)
When AIM is used, the program will perform
0: no update for secondary variables at explicit elements.
1: fully secondary variables update for all elements.
**SSRateInterpolation:** <int> (default is 0)
Determines interpolation procedure for time dependent well production rate or sink/source data (flow rates and enthalpies).
0: triple linear interpolation; tabular data are used to obtain interpolated rates and enthalpies for current time.
1: step function option; rates and enthalpies are taken as averages of the table values corresponding to the beginning and end of the time step.
The time in simulation time-stepping is enforced to be at any points of well rate change provided in the table.
**DoublingDtCriterion:** <int> (default is 4)
Provides automatic time step control based on model solution convergence behavior. Time step size will be doubled if convergence occurs within DoublingDtCriterion Newton-Raphson iterations. It is recommended to set this value in the range of 3- 5.
**SAVE_frequency:** <int> (default is 200)
Frequency (time steps) of writing the SAVE files.
**TSPrint_frequency**: <int> (default is 1)
Frequency (time steps) of writing time-series data.
**TimeUnit**: <int> (default is 0)
Defines the time unit for inputted times under the main keyword PARAM_. It will apply to all the following keywords “TimeShif”, “SimulationTimeEnd”, “initialTimeStep” “MaxTimeStep”, and “UserTimeStep”.
TimeUnit=
0: second
1: hour
2: day
3: year
TimeShif**: <float> (default is 0)
Beginning time of the simulation
**SimulationTimeEnd**: <float> (default is 1.0e15)
Time at the end of the simulation.
**InitialTimeStep**: <float> (default is 1000.0)
Initial time step size.
**MaxTimeStep**: <float> (default is 1.0e13)
The maximum time step size.
**TrackElemName:** <char>
Name of a tracked gridblock.
**UserTimeStep**: <float>
User Specified Timesteps
Input user specified time steps. User can define beginning time-step sizes of the simulation. The inputted list of time steps (any number) must be separated by comma “,”.
**gravityV**: <float> (default is 9.8)
Acceleration of gravity. If it is zero, no gravity is considered in this simulation. If it is a negative number, the z coordinates for model elements must be elevations (smaller value at bottom). If it is a positive number, the z coordinates for model elements must be depths (smaller value at top).
**DtReductionFactor**: <float> (default is 4.0)
Reduction factor for cutting time-step size in case of Newton-Raphson iteration convergence failure or other problems.
**MScale**: <float> (default is 1.0)
Scale factor to change the size of the mesh.
**ReOrderM**: <int> (default is 0)
Selects method for reordering element sequence of the model mesh.
ReOrderM=
0: no reordering.
1: Tarjan.
2: narrowest bandwidth.
**priv_convergence_crit**: <float>… (defaults are 100.0, 0.005, 0.005,……)
Convergence criteria for absolute error of primary variables. If changes of primary variables between two NR iterations are less than the given criteria, the solutions are considered to be converged. The number of data provided through this keyword can be as more as the number of primary variables for current simulation model. Sequence of the criterion list is corresponding to the sequence of primary variables. The inputted data must be separated by comma “,”.
rel_abs_error**: <float, float> (defaults are 1.0e-5, 1.0)
Convergence criteria for relative and absolute mass balance error.
**W_upstream:**
W_NRIteration**: <float> (default is 1.0)
Weighting factor for increments in Newton/Raphson – iteration (default = 1.0 is recommended). This parameter is in the range of 0 to 1.
**W_implicitness**: <float> (default is 1.0)
The factor for computing the implicitness - weighted time step increment.
**derivative_increment**: <float> (default is 1.0e-8)
Increment factor for numerically computing derivatives.
**P_overshoot:**<int> (default is 0)
The pressure overshooting allowed in triggering phase changes.
**T_overshoot:** <int> (default is 0)
The temperature overshooting allowed in triggering phase changes.
**S_overshoot:** <int> (default is 0)
The saturation overshooting allowed in triggering phase changes.
**default_initial_cond**: <float>…
Introduces a set of primary variables which are used as default initial conditions for all elements that are not assigned at elsewhere. The number of primary variables can be difference for different simulation modules. The list of primary variables must be separated by comma “,”.
If the corresponding phase status of this set of primary variables is known, it can be specified in two ways, by index or by name. For example, in a CO2 module simulation, an element has a liquid CO2 and aqueous two phases status. It is phase status is "L_A" and corresponding index is n=6. Suppose the element has CO2 pressure Pl=12MPa, CO2 saturation Sl=0.3, NACL mass fraction xs=0, and temperature T=66, the input can be in following ways:
By index
default_initial_cond: 1.2e7, 600.3, 0.0, 66 (the second variable has a value of sl+n*100)
By name
default_initial_cond: 1.2e7, 0.3, 0.0, 66, L_A
**Diff_FluxM:** <int> (optional, default is 0)
Selection Calculate diffusive flux method.
Example:
PARAM_
Maximum number of Newtonian iterations per time step, MaxNI: 6
Maximum number of time steps for the simulation run, MaxTSN: 4000
Write Output File Frequency, OutputFQ: 100
Numerical Method Selection (0 FIM, 1 AIM, 2 IMPES), Option_Aim: 0
Interpolation method for Time-Dependent Sink/source Data, SSRateInterpolation: 0
Criterion for Time-Step Size Increasing, DoublingDtCriterion: 4
Time at the beginning of the simulation, TimeShift: 0
Time at the end of the simulation, SimulationTimeEnd: 0.382798e+08
Initial time step size, InitialTimeStep: 100.0
Maximum time step size, MaxTimeStep: 1.0e9
Acceleration of gravity (gravityV): -9.8066
Reduction factor for cutting Dt, DtReductionFactor: 2.0
Convergence criteria for primary variables error, priv_convergence_crit:1000.0, 0.01, 0.01, 0.5
Convergence criterion for relative and absolute error, rel_abs_error: 1.0e-5, 1.0
AIM will perform fully secondary variables update, AimFullyUpdate: 0
Default initial values of the primary variables, default_initial_cond: 5.3e6, 0.0, 0.0, 0.0,65.0
# 3.4 GENER_, input for Well or Source/Sink data (buchong)
Source/Sink or well injection or production data are inputted through main keyword “GENER_”. Under the “GENER_” keyword, users are able to define any number of wells or source sink terms through following keywords.
Keywords:
**WellName**: < string>
Define the well or source/sink name. Length of the name is flexible. This keyword is always at the first line of inputs for each well and treats as the starting for a new well.
**ElementName:** < char>
Name of the mesh element containing the sink/source or the well. This element must be an active element of the models.
**WType:** < char>
Gives the well type which specifies different options for fluid or heat production and injection. For example, different fluid components may be injected, the nature of which depends on the EOS module being used. Different options for considering wellbore flow effects may also be specified.
| WType | Explanation |
|---|---|
| HEAT | Heat source/sink |
| COM1/Gas | Gas source/sink |
| COM2/Water | Water source/sink |
| COM3 | Nonaqueous phase/the second liquid phase source/sink |
| QLOW | Pumping well with total liquid pumping rate specified |
| QOIL | Pumping well with oil pumping rate specified |
| QWAT | Pumping well with water pumping rate specified |
| QGAS | Pumping well with gas pumping rate specified |
| QJQW | Injection well with water injection rate specified |
| QJQO | Injection well with oil injection rate specified |
| QJQG | Injection well with gas injection rate specified |
| DELV | Pumping well with a wellbore pressure specified |
| QJPW | Water injection well with injection pressure specified |
| QJPO | Gas injection well with injection pressure specified |
| MASS | Specifies mass production rate |
| CLOS | Close the well |
| SHUT | Close the well temporarily |
| GENE | Hybrid well |
**IPRATE:** < float>
Constant volumetric injection (q>0)/production (q<0) rates with sink/source generation terms. The units are usually kg/s or m3/d. Well volumetric pumping or injection rate at standard conditions (m3 /d) of fluids.
**PINDEX:** < float>
Production Index for specified bottom hole pressure. Its value is equal to the fluids volume for per bottom hole pressure.
**SENTHALPY:** < float>
Specific enthalpies of injected fluid will be read [J/kg]. The fluid type can be referred from WType.
**BPRESSURE:** < float>
Bottom hole pressure
Specified maximum well injection pressure (Pa) for an injection well with rate specified (default = infinite);
pumping or injecting pressure (Pa) for a pumping or injection well with wellbore pressure specified.
**LTHICKNESS:** < float>
The thickness of the well layer [m], length of well is related to the WType and the pumping well and injecting well is usually different.
**MPRESSURE:** < float>
The maximum (for injection) or minimum (for production) allowable pressure at the well. Specified maximum well injection pressure (Pa) for an injection well with rate specified (default = infinite); Specified minimum well flowing pressure for a pumping well with rate specified (default = reference pressure from data block).
**WSINDICATOR:** < float>
Water saturation indicator for well closure.
**WELLNOTE:** < char>
Detail notes for the well. it provides more detailed information about wells.
**VMFactor:** < float>
Volume-increasing factor for wellbore element (default = 1).
**TLength:** < int>
Table Length of time-dependent well operation changes.
**WRUnit:** < char>
Well production/injection rate unit length. Unit often is ‘kg/s’ or ‘kg/d’.
**TUint:** < char>
The time unit for well operation.it can be 'S'-Second, 'H'-hour, 'D'-day, 'M'-Month, 'Y'-year; default is 'S'
**TimeType:** < char>
The input time series type for well operation. It can be 'S' or 'C'; S-input time step, 'C'-input increasing time series, default is 'C'.
**MeasuredData:** < char>
The input production data are real measured data:
**MeasuredBHP:** < char>
The input well bottle pressure data are real measured data.
PHASEFRAC:
The mass fraction of each component of the mixture.
PHASECOMP:
Mix injection gas components. Like hydrogen, carbon dioxide, etc.
WOTS:
Well operation time series (WOTS): 3.100E+01, QLOW, 8.870E+00, 1.000E+05
*Note:* 31-monthor pumping with total liquid pumping rate specified of 8.87kg/s, the pumped fluid enthalpy is 1.00e5 J/kg.
Well operation time series (WOTS): 3.000E+01, QLOW, 2.670E+00, 1.000E+05
Well operation time series (WOTS): 3.000E+01, QLOW, 1.500E+00, 1.000E+05
Well operation time series (WOTS): 2.900E+01, SHUT, 0.000E+00, 1.000E+05
*Note:* shut-in for 29-month
Example:
GENER
The Well Name (WellName): WELL1
Element Name for the Well (ElementName): A1 1
Type of Well or Source/Sink (WType): QWAT
Injection/production rate (IPRate): 5.0
Production Index for specified bottom hole pressure (PIndex): 4.0e-12
Specific enthalpy of injected fluid (SEnthalpy): 2.5
Bottom hole pressure (BPressure): 1.25e6
Well layer thickness (LThickness): 10
The maximum (for injection) or minimum (for production) allowable pressure at the well (MPressure): 5.0e7
Water saturation indicator for well closure (WSIndicator):
Detail notes for the well (WellNote):
Volume multiplication factor for the well element (VMFactor):
Table Length of time-dependent well operation changes (TLength):
Well production/injection rate unit (WRUnit): KGS
The time unit for well operation (TUint): S
The input time series type for well operation (TimeType): S
The input production data are real measured data (MeasuredData): F
The input well bottle pressure data are real measured data (MeasuredBHP): F
PHASEFRAC
Well operation time series (WOTS): 3.100E+01, QLOW, 8.870E+00, 1.000E+05
Well operation time series (WOTS): 3.000E+01, QLOW, 2.670E+00, 1.000E+05
Well operation time series (WOTS): 3.000E+01, QLOW, 1.500E+00, 1.000E+05
Well operation time series (WOTS): 2.900E+01, SHUT, 0.000E+00, 1.000E+05
# 3.5 TIMES / Specifies the printout times.
This keyword permits the user to obtain printout at specified times (optional). This printout
will occur in addition to printout specified in record PARAM.
Keywords:
**OutputTime: <**float**>**
List of times (in ascending order) at which printout is desired.
**MaxPrintTimes:** <int>
Total number of times desired in OutputTime.
OPTimeIncrement:
Time increment for times of output.
TimeUnit:
Unit of the list times:0-seconds,1-hours, 2-days, 3-years
Example:
TIMES**
Write output at times (OutputTime):3.15576E7,1.2559E8,3.15576E8
Total number of times (MaxPrintTimes): 50
Time increment of output (OPTimeIncrement): 1.0E6
Unit of the list times (0-seconds,1-hours, 2-days, 3-years, TimeUnit): 0
# 3.6 ELEME / Defines mesh element.
This block mainly introduces the element (grid block) information. It is used for pre-processing of the model.
Keywords:
Eleme Name:** <string>
Code name of an element. It consists of n characters (3≤n≤13).
**Rock Name:** <string>
Rock name, the length of the name must be equal or less than 5 characters. This keyword always be the first line of inputs for each rock.
**Volume:** <int>
The volume of the element mentioned above, the unit is stere [m3].
**X, Y, Z: <int>**
Cartesian coordinates of grid block centers. These may be included in the ELEME data to make subsequent plotting of results more convenient.
Activity:
Four cases are contained: ‘A’, ‘I’, ‘V’ or blank. ‘A’ is default, namely active grid. ‘I’ is inactive grid. ‘V’ is Dirichlet boundary gird whose conditions can vary with time.
Example:
ELEME**
8AAAA, sand1, 2.8846E+07, -1.5000E+4, -1.0000E+4, -3.2084E+3, A
# 3.7 CONNE / Defines mesh connection information.
This block introduces information for the connections (interfaces) between elements.
Keywords:
Element1:
Code name of the first element in the connection.
Element2:
Code name of the second element in the connection.
ISO:
set equal to 1, 2, or 3; specified absolute permeability to be PER(ISO) for the materials in elements (Element1) and (Element2), where PER(ISO) is read in block ROCKS. PER(1)= permInX, PER(2)= permInY, and PER(3)= permInZ.
D1, D2:
The distance (m) from first and second element, respectively, to their common interface.
Area:
Interface area [m2]
BETA:
cosine of the angle between the gravitational acceleration vector and the line between the two elements. BETAX> 0 (<0) corresponds to first element being above (below) the second element.
Example:
CONNE**
AA1 6AB1 6 10.2500E+000.2500E+000.2000E+020.1000E+01
# 3.8 INCON/ Specifies initial conditions for specific elements.
This block introduces element-specific initial conditions.
Keywords:
ELEMENTNAME:
Code name of an element. It consists of n characters (3≤n≤13).
X1, X2, X3, X4:
Specifies the initial value for the primary variables
Example:
INCON**
Element name (**ELEMENTNAME**): AB1 3
**X1, X2, X3, X4:** 1.013e+05,0.0,1.0007181548e+01,2.5e+01
# 3.9 SOLVE/ Specifies related parameters for solving linear equations.
This block specifies parameters used by linear equation solvers.
Keywords:
SLibrary: (Default is AMGCL)
This variable selects the linear equation solver from among the following options:
AMGCL, PETSC, AMCL, FASP, TRILINOS, AMGX, VIENNACL, rocALUTION)
Pre_C:
The type of the Preconditioner selection. The impact of the selection on solution speed can be obvious. Different libraries have different preconditioned methods, and the following options are allowed in the current version.
AMGCL: CPR, CPD_DRS, AMG, ILU and SCHUR (Default is CPR and ILU when the numerical method is FIM and AIM, respectively.)
TRILINOS: JACOB, DOMDC, AMG and ILU (Default is DOMDC)
PETSC: JACOB, ILU, DOMDC, AMG, EUCLID and PARASAILS (Default is JACOB)
FASP: ILU and AMG (Default is ILU)
Note that AMG may not work for off-diagonal dominance, CPR may not work for AIM, and DOMDC may not work for large scale models (The number of grids is over millions).
LSolver: (Default is BICGSTAB)
Linear solver selection. Options include BICGSTAB and GMRES.
Maxit: (Default is 300)
Maximum iteration number.
c_crit: (Default is 1.0e-5)
Linear solver convergence criteria.
AIM_sa
AIM switching saturation criteria.
DoMultiD:
Example:
SOLVE**
Select Solver Library (SLibrary-AMGCL, PETSC, AMCL, FASP): A
PreConditioner Selection (Pre_C): ILU
Linear Solver Selection (LSolver): GMRES
Maximum Iteration Number (Maxit): 300
Linear solver convergence criteria (c_crit): 1.0e-3
AIM switching saturation criteria (AIM_sa): 0.01
/
# 3.10 FOFT/ monitor a list of elements.
This block output the information of the specified element, for which time-dependent data are to be written out for plotting.
Keywords:
**OBElements:** <char>, <char>.… [—]
This usually is the element name, separated by commas.
**ObsElemVar:** <int>, <int>, … [—]
The define of serial number of variables, separated by commas. Up to 10 numbers.
Example:
FOFT**
Element name (**OBElements**): AF1 1, A61 3
Serial number of variables (**ObsElemVar**): 1,2,3,4,5,6,7
/
# 3.11 COFT/ monitor a list of connections.
Specifies a list of connections for which time-dependent data are to be written out for plotting.
Keywords:
**OBConnections**: <char>, <char>… [—]
Connection name, separated by commas
**ObsConxVar:** <int>, <int>, [—]
Serial number of variables, separated by commas. Up to 10 numbers.
1 2 3
Example:
COFT**
Connection name (OBConnections): AC118AC119
Serial number of variables (ObsConxVar): 1,2,3,4,5,6,7,8,11
/
# 3.12 GOFT/ monitor a list of sinks/sources.
Specifies a list of sinks/sources for which time-dependent data are to be written out for plotting.
Keywords:
**OBSS:** <str>, <str>
Element name for the sinks/sources, separated by commas.
ObsVar:
Observe value
Example:
**GOFTS:**
**observation at wells (OBSS)**: 1L162,1L160,1F1H2, 1F1H0, 1P1I2, 1P1I0, 1E132,1E130, 1W1G2, 1W1G0,2A1A2, 2A1A0,260Y2, 260Y0,1K1F2, 1K1F0, 2A122, 2A120, 181D2, 181D0, 231C2, 231C0,1A1K2, 1A1K0
**observation variables including (ObsVar)**: 1, 2,3,4,5, 6, 7, 8, 9, 10, 11, 12, 13
# 3.13 MODEL/ Specifies related parameters for model design.
This block is used to specify the related parameters for the model design.
Keywords:
**HasAqu2:** <char>
The flow system has Aqu2. Ture or False. If it is true, the fluid properties must be specified in keyword FLUID and/or PVTDA.
**GasAqu2:** <char>
Including Gas simulation. True or False. If it is true, the gas properties must be specified in keyword FLUID and/or PVTDA.
**AquMixable:** <char>
Water and Aqu2 is mixable. True or False.
**ElemNameLength:** <int> (default is 5)
Length of element name.
RelaxedV:** <char> (default is false)
Using relaxed volume. Ture or False.
**AFDiffusion:** <char> [—] (default is false)
Accounting For Diffusion. Ture or False.
**PermUnit:** <int> [—]
Rock Permeability Unit (0 for m-2; 1 for D; 2 for mD).
**nIsoTherm:** <char> [—]
Non isothermal simulation. Ture or False.
**EosType****：**
Type of Equation State. "PR", "SRK" or "RK".
PR: Peng–Robinson equation of state
SRK: Soave modification of Redlich-Kwong
RK: Redlich-Kwong equation of state
This keyword is only needed if the gas phase is a non-condensible gas.
**RealGAcV:** <char> (default is true)
Real Gas Including Vapor. Ture or False.
GasCompo:
MassInMol:
DELVInLV:
CO2PModel:
Example:
MODEL**
The flow system has Aqu2 (HasAqu2): FALSE
Including Gas simulation (HasGas): TRUE
Water and Aqu2 is mixable (AquMixable): FALSE
Length of Element Name (ElemNameLength): 5
Using Relaxed Volume (RelaxedV): FALSE
Accounting For Diffusion (AFDiffusion): FALSE
Rock Permeability Unit (PermUnit): 0
Non Isothermal Simulation (nIsoTherm): TRUE
Type of Equation State (EoSType, can be one of "PR", "SRK", "RK"): PR
Real Gas Including Vapor (RealGAcV): TRUE
GasCompo:
/
# 3.14 RPCAP/ relative permeability function and capillary pressure function.
This block specifies default function type and related parameters for relative permeability function and capillary pressure function.
Keywords:
**DRPT:** <int>
Default relative permeability function type.
**DCPT:**<int>
Default capillary pressure function type.
**DRPParameters:** <float>, <float>, …
Default relative permeability function parameters.
**DCPParameters:** <float>, <float>, …
Default capillary pressure function type.
Example:
RPCAP**
Default Relative Permeability Function Type (DRPT): 7
Default Capillary Pressure Function Type (DCPT): 7
Default Relative Permeability Function Parameters (DRPParameters): 0.45000,9.6E-4,1.,,,,
Default Capillary Pressure Function Type (DCPParameters): 0.45000,1.0E-3,8.0E-05,5.E8,1.,,
# 3.15 FLUID/ Specifies related parameters about fluid properties.
This block mainly defines the related parameters about the fluid properties.
Keywords:
**FluidName:<char>**
Define of the model simulation fluid name.
**FluidDensity:<float>**
The Fluid density at standard condition. [kg/m3]
**FluidCom: <float>**
The factor of the Fluid compressibility. [1/Pa]
**FluidVis: <float>**
Define of the Fluid viscosity at standard condition. The unit is Pa.S.
**ThermalExans:<float>**
Thermal expansion coefficient for corresponding fluids.
**FluidBingG: <float>**
When fluid is considered as Bingham fluid, it is parameter G (PA / M).
**FluidPowN: <int>**
When fluid is considered to power-law fluid, it is Power-Law index N, otherwise it is set to 0.
**FluidPowH:<float>**
When fluid is considered as power law fluid, it is Power-Law parameter H, otherwise it is set to 0.
**FluidMinVis:<float>**
When power law fluid is considered, it is the minimum value of the power discharge fluid viscosity (Pa.s) (default is 1.0e-5 Pa.s); if the fluid is a non-power rhythm its value is 0.
**FluidType:<int>**
Define the type of fluid. Three types are included: Newtonian-0, PowerLaw-1, Bingham-2
**referT:<float>**
Fluid reference temperature[℃].
**referP:<float>**
Fluid reference pressure [Pa].
**ThermalCondc:<float>**
Fluid thermal conductivity. [W/m^2]
Example:
FLUID**
The model Simulation Fluid Name (**FluidName**): FLUID
Fluid density at standard condition (**FluidDensity**): 600
Fluid compressibility factor (**FluidCom**): 2.0
Fluid viscosity at standard condition (**FluidVis**): 0.5
Thermal expansion coefficient for the Fluid (**ThermalExansionC**): 2.0
Parameter G (**FluidBingG**):
Power-law index (**FluidPowN**):
Power-law parameter (**FluidPowH**):
Minimum viscosity for a power-law fluid and for a non-power-law fluid set to zero (**FluidMinVis**): 1.0E-5
FluidType (**FluidTpye,** Newtonian-0, PowerLaw-1, Bingham-2): 0
Parameter of minimum potential gradient (**FluidBingG**):
Fluid reference temperature (**referT**): 35
Fluid reference pressure (**referP**): 1.0e5
Fluid Thermal Conductivity (**ThermalCondc**): 2.4
# 3.16 PVTDA/ Specifies PVT data.
This block specifies PVT data.
Keywords:
**PVTatT:** **<float>, <float>, <float>**
PVT at Temperature.
**PB_I:<float>**
Initial Bubble Pressure for This Group.
**BW_I:<float>**
Water Formation Volume Factor.
**PVT_TABLE:** **<float>, <float>, <float>, <float>, <float>, <float>**
Pressure, Oil Formation Volume Factor, Gas Formation Volume Factor, Oil Gas Ratio, Oil viscosity, Gas Viscosity.
Example:
PVTDA**
PVT at Temperature (**PVTatT**):
Initial Bubble Pressure (**PB_I**):
Water Formation Volume Factor (**BW_I**):
Pressure, Oil Formation Volume Factor, Gas Formation Volume Factor, Oil Gas Ratio, Oil viscosity, Gas Viscosity (**PVT_TABLE**):
# 3.17 ENDCY/ End of input data.
It specifies the close the input file and start simulation.
# 3.18 ENDFI_/ data checking and marking.
Keywords that do not perform flow calculation simulation, which is often used to check input data or network Grid generation.
# 3.19 TIMBC_/The Gridblock name with time-dependent boundary condition.
This indicates the input time-dependent boundary condition.
Keywords:
Ele_Name:
Element name
N_PrimaryV:
Number of times
TimeUnit:
Time unit (1-minutes, 2-hours, 3-days, 4-years)
TimeBC:
Time of primary variable
Primary
Primary variable value.
Example:
Ele_Name: A1V29
N_PrimaryV: 1
TimeUnit: 3
timeBC: 0., 182.0, 365.0, 730.0
primaryVA: 0.8666267541935E+07, 0.8766267541935E+07, 0.8866267541935E+07, 0.8866267541935E+07
# 3.20 INDOM/ Specifies initial conditions for specific rocks.
This block specifies initial conditions for specific rocks.
Keywords:
**Rock_Name:** <char>
Define the rock name.
**IniCond:** <float>
This block introduces domain-specific initial conditions.
Example:
INDOM**
Rock name (**Rock_Name**): SAND1
Initial conditions (**IniCond**): 1.E5,0.0,10.2, 18.
# 3.21 DIFFU/ Specifies diffusion coefficients for different component at different phases.
This block illustrates the diffusion coefficients for different component at different phases.
Keywords:
**NoDif:<float>**
Dispersion is not considered.
**DiffC:<float>**
Diffusion coefficients for the first mass component at different phases.
Example:
DIFFU**
Diffusion Coefficients for first mass component at different phases (**DiffC**): 2.13e-5,0.e-8
Second mass component (**DiffC**): 2.13e-5,0.e-8
# 3.22 SVPAR/ Specifies the array input of element property.
Description：Array input for unit properties (porosity, absolute permeability, and so on).
# 3.23 TRCON/ Specifies the time-dependent connections information between elements.
Description：Specifies the time-dependent connections information between elements, such as opening or closing a connection at a certain time.
Keywords:
**CPair:<char>**
Connection name.
**t_openf:<float>**
The time the connection is opened
**t_unit:<char>**
Time unit.
# 3.24 OUTCT/ Specifies printout information.
The block specifies the printout information.
Keywords:
OutputFQ:** **<int>
Frequency (time step) of writing output file.
**SAVE_frequency:<int>**
Define the frequency of writing the SAVE files (default is 200).
**TimeSeriesPrint_frequency:<int>**
Output frequency for the time-dependent data.
**NumOutPutDataFile:<int>**
Number of output files
**NumOutPutSAVE:<int>**
Number of Save files
**ELEVA:<int>**
Serial number of variables which will be output for elements. Stop when blank space or serial number greater than 49 is read. By default (that is, ELEVA after unordered number data), select the first 6 variables. Variable selection code is: 1 pressure, 2 temperature, 3 Gas saturation, 4 water saturation, 5.gas mass fraction in Gas, 6. gas mass fraction in aqueous, 7. capillary pressure gas/water, 8. Porosity 9. X permeability, 10. Y permeability, 11. Z permeability, 12. Coordinates X, 13. Coordinate Y, 14. Coordinate Z, 15. Partition: domain decomposition region, 16. Oil-gas capillary pressure, 17. Gas viscosity, 18. Oil viscosity, 19. Gas Enthaply, 20. Gas density, 21. Water density. 22. Oil density. 23. Gas mass fraction in Oil. 24. Gas2 mass fraction in Gas, 25. Gas2 mass fraction in aqueous, 26. Gas2 mass fraction in oil, 27. Gas3 mass fraction in Gas 28. Gas3 mass fraction in aqueous, 29. Gas3 mass fraction in oil 30. Gas4 mass fraction in Gas 31. Gas4 mass fraction in aqueous 32. Gas4 mass fraction in oil 33. Water mass fraction in Gas 34. Water mass fraction in aqueous 35. Water mass fraction in oil 36. Oil saturation 37, water viscosity 38.water enthalpy 39. Oil ehthalpy 40. The relative permeability of gas 41. The relative permeability of water 42. The relative permeability of oil. 43. Solidphase saturation
\44. Solidphase enthalpy 45. Solidphase density 46. Gas pressure 47. Water pressure 48. Oil press
\49. Capiliary pressure of auq-oil.
**CONVA:<int>**
Serial number of variables which will be output for connections. Stop when blank space or serial number greater than 14 is read. Variable selection code is: 1. Gas FluxRate, 2.Water FluxRate, 3. Fluid Enthalpy, 4.Total Flux, 5.Heat Flux, 6. Por Vel (Gas), 7. Por Vel (Aqu), 8. PorVel(Aqu2), 9. Darcy Vel(Gas), 10. Darcy Vel(Aqu), 11.Darcy Vel (Oil), 12. X, 13. Y, 14. Z, 15. Oil FluxRate.
**ELETS:<int>**
Serial number of variables which will be output for time series data of elements. The corresponding keywords can be referred from ELEVA.
**CONTS:<int>**
Serial number of variables which will be output for time series data of connections. The corresponding keywords can be referred from CONVA.
**SASVA:<int>**
Serial number of variables which will be output for Sink and source term. The Variable selection code is the following: 1. "Gas Rate", 2."Water Rate", 3."Aqu2 Rate", 4."Heat rate", 5."Well Pressure", 6. "Cum. Gas", 7. "Cum. Water", 8."Cum. Oil", 9. Water saturation 10. Oil saturation 11. Gas capillary pressure 12. Oil capillary pressure 13.gas viscosity 14. Oil viscosity
The corresponding Sink Source_VariablesTerm_Units are the following: "kg/s", "kg/s", "kg/s", "J/s", "Pa", "kg", "kg", "kg", "kg/s", "kg/s", "kg/s", "kg", "kg", "kg", "kg".
**SASTS:<int>**
Serial number of variables which will be output for time series data of Sink and Source terms. The corresponding keywords can be referred from SASVA.
Example:
OUTCT_
output element variables max10(ELEVA): 1, 3, 5, 6, 24, 25, 27, 28, 20, 17
# 3.25 PREME/ Preprocessing of model mesh.
This part mainly introduces the key information involved in the grid preprocessing.
Keywords:
**MeshAction:<**char**>**
Mesh preprocessing Action. MeshAction=:
“ActiveGrid”: Exporting basic grid information statistics, remove inactive grids from MESH grids, link and save them to disk。
“EclMesh”: Converting the ECLIPSE corner point grids from Petrel and primary input contents to the “Karatsim” format.
“MINC”: Converting specific lithologies from current single-media MESH files to multiple-media meshes for multiple media.
“LinkMESH”：Connecting the elements generated by chunking in large-scale simulations into a whole mesh grid automatically
**KWFile:**<char>
Keyword and its filename (Global path of the file corresponding to the keyword). The keyword is necessary when MeshAction = ‘EclMesh’. Two strings, separated by commas. The keywords contained in the files exported by Petrel, containing: "COMPDAT","COORD","SPECGRID", "ZCORN", "PORO", "PERMX", and "OIL_WATER".
**RockMinC:** <char>
List the rock lithology for MINC processing. Up to 10 lithologies can be entered. The mesh with these lithologies will be converted to a multi-media mesh, and if the first lithology entered is "ALL", all lithologies in the mesh will be converted to multi-media.
**FTYPE**: <char>
A five-character word for selecting one of the three different conceptual models of fractures:
‘ONE-D’: a set of plane parallel infinite fractures with matrix block thickness between neighboring fractures equal to PAR(l)..
‘TWO-D’: two sets of plane parallel infinite fractures, with arbitrary angle between them. Matrix block thickness is PAR(l) for the first set, and PAR(2) for the second set. If PAR(2) is not specified explicitly, it will be set equal to PAR(l).
‘THRED’: three sets of plane parallel infinite fractures at right angles, with matrix block dimensions of PAR(l), PAR(2), and PAR(3), respectively. If PAR(2) and/or PAR(3) are not explicitly specified, they will be set equal to PAR(l) and/or PAR(2), respectively.
**DualFType:<**char**>**
The inter-media flow type. ‘MMVER’, ‘MMALL’ or blank.
‘blank’: (default) (default) global flow occurs only through the fracture continuum, while rock matrix and fractures interact locally by means of "interporosity" flow ("double-porosity" model).
‘MMVER’: global matrix-matrix flow is permitted only in the vertical; otherwise like the double-porosity model; for internal consistency this choice should only be made for flow systems with one or two predominantly vertical fracture sets.
“MMALL”：global matrix-matrix flow in all directions; for internal consistency only two continua, representing matrix and fractures, should be specified ("dual-permeability").
**NumInterM**: <int>
Total number of multi-correlated continuous media nesting (J).
**NVOL**:<int>
The total number of volume fractions (NVOL < J) is given explicitly. If NVOL < J, the total number of volume fractions for the volume fractions with subscripts NVOL+l, ... , the volume fractions of J are aliquoted to generate.
**LocationBG:<**float**>**
Specified volume fractions begin from fractures or media internal. ‘IN’ or ‘OUT’.
**PARAMS:**<float>
Parameters constraining the distribution of cracks. Up to 7 numbers can be entered, separated by commas: PRA(I) I=1,…,7
**VolumeFList:<**float>
Volume fraction of continuous media (from 0 to 1). When LocationBG = 'OUT ', specify the volume fraction of I. When LocationBG = 'IN', specify the volume fraction of J+ l - I. If LocationBG = 'OUT', I = 1 is the fracture continuum medium, I = 2 is the matrix medium closest to the fracture, I = 3 is the matrix medium adjacent to I = 2, and so on. The sum of all volume factors must not exceed 1.
Example:
PREME**
Mesh preprocessing Action (**MeshAction**): EclMesh
Keyword and its location (**KWFile**):
List rocks for MINC processing (**RockMinC**): ROCK2
The type of fracture media (**FTYPE**): ONE-D
The inter-media flow type (**DualFType**): MMVER
Number of multiple interacting continua (**NumInterM**): 5
Number of explicitly provided volume fractions (**NVOL**): 4
Specified volume fractions begin from fractures or media internal (**LocationBG**): 0.3
Fracture distribution parameters (**PARAMS**):
Volume fraction of continuum (**VolumeFList**): 0.2
# 3.26 COMPO_/Specifies some parameters of Composition model.
This block mainly defines some crucial parameters of the composition model.
Keywords:
CompoName**：<char>
Define the name of compositions. Three common gases are usually included here: hydrogen/nitrogen and methane, etc.
**CompoType**: <char>
Define the type of compositions. Different letters represent different input component types and the three common component types are as follows:
‘h’: hydrocarbon
‘g’: gas
‘s’: solvent
**Crit_TPZV**: <double>, <double>, <double>, <double>
Critical temperature, Critical pressure, Critical compression factor, and Critical volume.
Omegam :
Define some parameters of the component model.
Dipolm:
NBoilingP:
VPConstants:
MoleW:
Mole fraction.
HCConstants:
RefParams:
DiffuExp:
ViscoConstants:
SolubConstants:
CPartitionC:
FractionOC
DecayC
Parachor**：**
Define the parachor of different compoents.
| Component | Parachor |
|---|---|
| 41.0 | |
| 78.0 | |
| 80.1 | |
| 77.0 | |
| 108.0 | |
| 150.3 | |
| 181.5 | |
| 203.4 | |
| 225.0 | |
| 231.5 | |
| 271.0 | |
| 312.5 | |
| 351.5 | |
| 393.0 | |
| 446.2 |
PeneluxVC
BIJ
# 3.27 GMESH_/ the keywords including in the model construction software.
This block is about the pre-processing of the grid construction, the detail processes modeled using software GMESH.
Keywords:
**MeshType:** <str>
This variable indicating the coordinate system used in the geometric definition of the region. In may have one of two values:
‘Cylindrical’: this option indicates the region is defined geometrically in terms of Cartesian coordinates.
‘Cartesian’ : this option indicates the region is defined geometrically in terms of cylindrical coordiantes.
**ElemNameL**: <int>
The element name.
**LengthUnit**: <int>
Defines the length unit for the element.
**Dangle:** <int>
The angle between horizontal layer and Z axil
Origin_P** <int,int>
Coordinates of the origin.
Const_DX
Equal length division in X direction.
Const_DY:
Equal length division in Y direction.
Const_DZ
Equal length division in Z direction.
Vary_DX
Variable length division inX direction.
Vary_DY
Variable length division in Y direction.
Vary_DZ
Variable length division in Z direction.
Const_DR
Equal length division in R direction.
Vary_DR
Variable length division in R direction.
Log_DR:
Logarithmically divided along the R direction
RAD_R:
The value of a given radius.
Const_LYThick:
The sublayer is of equal thickness.
Vary_LYThick:
The sublayer is of variable thickness.
MeshSeq:
The Mesh sequence.
FilChar:
3-14 characters.
ZTopType:
It represents whether the top of the model is horizontal or undulating.
FormationInfo:
The information of the formation.
MRegion:
Define the region area.
RBound:
Define the boundary condition.
StartingR:
Define the initial radius value.
EXAMPLE:
ElemNameL:5
MeshType: xyz
Const_DY: 10,1.0
Const_DX: 2,1
Const_DZ: 2,1
meshSeq: 1
The angle between horizontal layer and Z axil (rotate around Y axil),
DAngle: -45 //positive for incline, negative for decline
FilChar: p
MRegion: SOIL1,cartesian, 0, 20, 0, 50, -1, 0
MRegion: SOIL2,Cartesian, 0, 20, 0, 50, -5, -1
RBound: I,BROCK,Cartesian, 0, 1, 0, 4, -1, 0
RBound: T,BROCK,Cartesian, 0, 1, 6, 10, -1, 0