Aquifer Boundary Condition

Overview

Aquifer boundary conditions allow simulating flow between the computational domain (the reservoir) and one or multiple aquifers. In GEOS, we use a Carter-Tracy aquifer model parameterized in Aquifer tags of the FieldSpecifications XML input file blocks.

Aquifer model

An aquifer $A$ is a source of volumetric flow rate $q^A_f$ , where $f$ is the index of a face connecting the aquifer and the reservoir. We use a Carter-Tracy model in GEOS to compute this volumetric flow rate.

Once $q^A_f$ is computed, the aquifer mass contribution $F^A_f$ is assembled and added to the mass conservation equations of the reservoir cell $K$ connected to face $f$ . The computation of $F^A_f$ depends on the sign of the volumetric flow rate $q^A_f$ .

The upwinding procedure is done as follows: if the sign of $q^A_f$ indicates that flow goes from the aquifer to the reservoir, the aquifer contribution to the conservation equation of component $c$ is:

$F^A_{f,c} = \rho^A_w y^A_{w,c} q^A_f$

where $\rho^A_w$ is the aquifer mass/molar water phase density and $y^A_{w,c}$ is the aquifer mass/molar fraction of component $c$ in the water phase. We assume that the aquifer is fully saturated with the water phase.

If the sign of $q^A_f$ indicates that flow goes from the reservoir into the aquifer, the aquifer contribution to the mass/molar conservation equation of component $c$ is computed as:

$F^A_{f,c} = \sum_{\ell = 1}^{n_p} ( \rho_{\ell} S_{\ell} y_{\ell,c} )_K q^A_f$

where $n_p$ is the number of fluid phases, $(\rho_{\ell})_K$ is the reservoir cell mass/molar density of phase $\ell$ , $(S_{\ell})_K$ is the reservoir cell saturation of phase $\ell$ , and $(y_{\ell,c})_K$ is the reservoir cell mass/molar fraction of component $c$ in phase $\ell$ .

In the next section, we review the computation of the aquifer volumetric flow rate $q^A_f$ .

Carter-Tracy analytical aquifer

The Carter-Tracy aquifer model is a simplified approximation to a fully transient model (see R. D. Carter and G. W. Tracy, An improved method for calculating water influx, Transactions of the AIME, 1960).

Although the theory was developed for a radially symmetric reservoir surrounded by an annular aquifer, this method applies to any geometry where the dimensionless pressure can be expressed as a function of a dimensionless time.

The two main parameters that govern the behavior of the aquifer are the time constant and the influx constant. These two parameters are precomputed at the beginning of the simulation and are later used to compute the aquifer volumetric flow rate.

Time constant

The time constant, $T_c$ , has the dimension of time (in seconds).

It is computed as:

$T_c = \frac{\mu^A_w \phi^A c_t^A (r^A_0)^2}{k^A}$

where $\mu^A_w$ is the aquifer water phase viscosity, $\phi^A$ is the aquifer porosity, $c_t^A$ is the aquifer total compressibility (fluid and rock), $r^A_0$ is the inner radius of the aquifer, and $k^A$ is the aquifer permeability.

The time constant is used to convert time ( $t$ , in seconds) into dimensionless time, $t_D$ using the following expression:

$t_D = \frac{t}{T_c}$

Influx constant

The influx constant, $\beta$ , has the dimension of $m^3.Pa^{-1}$ .

It is computed as:

$\beta = 6.283 h^A \theta^A \phi^A c^A_t (r^A_0)^2$

where $h^A$ is the aquifer thickness, $\theta^A$ is the aquifer angle, $\phi^A$ is the aquifer porosity, $c_t^A$ is the aquifer total compressibility (fluid and rock), and $r^A_0$ is the inner radius of the aquifer.

Aquifer volumetric flow rate

Let us consider a reservoir cell $K$ connected to aquifer $A$ through face $f$ , and the corresponding aquifer volumetric flow rate $q^A_f$ over time interval $[ t^n, t^{n+1}]$ .

The computation of $q^A_f$ proceeds as follows:

$q^A_f = \alpha^A_f ( a - b ( p_K( t^{n+1} ) - p_K( t^n ) ) )$

where $\alpha^A_f$ is the area fraction of face f, and $p_K( t^{n+1} )$ and $p_K( t^n )$ are the pressures in cell $K$ at time $t^{n+1}$ and time $t^n$ , respectively.

The area fraction of face $f$ with area $|f|$ is computed as:

$\alpha^A_f = \frac{ |f| }{ \sum_{ f_i \in A } |f_i| }$

The coefficient $a$ is computed as:

$a = \frac{1}{T_c} \frac{ \beta \Delta \Phi^A_K( t^n_D ) - W^A( t^n_D ) P_D^{\prime}( t^{n+1}_D ) }{ P_D ( t^{n+1}_D ) - t^{n+1}_D P_D^{\prime} ( t^{n+1}_D ) }$

and the coefficient $b$ is given by the formula:

$b = \frac{1}{T_c} \frac{\beta}{ P_D( t^{n+1}_D ) - t^{n+1}_D P^{\prime}_D( t^{n+1}_D ) }$

where $\Delta \Phi^A_K( t^n_D ) := p^A - p_K( t^n ) - \rho^A_w g ( z_K - z^A )$ is the potential difference between the reservoir cell and the aquifer at time $t^n$ , $P_D( t_D )$ is the dimensionless pressure evaluated at dimensionless time $t_D$ , $P^{\prime}_D( t_D )$ is the derivative of the dimensionless pressure with respect to dimensionless time, evaluated at dimensionless time $t_D$ .

The functional relationship of dimensionless pressure, $P_D$ , as a function of dimensionless time is provided by the user. A default table is also available, as shown below. The cumulative aquifer flow rate, $W^A( t^n_D )$ , is an explicit quantity evaluated at $t^n_D$ and updated at the end of each converged time step using the formula:

$W^A( t^{n+1}_D ) = W^A( t^{n}_D ) + ( t^{n+1} - t^{n} ) \sum_{f \in A} q^A_f$

with $W^A( 0 ) := 0$ .

Parameters

The main Carter-Tracy parameters and the expected units are listed below:

aquiferPorosity: the aquifer porosity $\phi^A$ .
aquiferPermeability: the aquifer permeability $k^A$ (in m2).
aquiferInitialPressure: the aquifer initial pressure $p^A$ (in Pa), used to compute $\Delta \Phi^A_K$ .
aquiferWaterViscosity: the aquifer water viscosity $\mu^A_w$ (in Pa.s).
aquiferWaterDensity: the aquifer water mass/molar density $\rho^A_w$ (in kg/m3 or mole/m3).
aquiferWaterPhaseComponentNames: the name of the components in the water phase. These names must match the component names listed in the fluid model of the Constitutive block. This parameter is ignored in single-phase flow simulations.
aquiferWaterPhaseComponentFraction: the aquifer component fractions in the water phase, $y^A_{w,c}$ . The components must be listed in the order of the components in aquiferWaterPhaseComponentNames. This parameter is ignored in single-phase flow simulations.
aquiferTotalCompressibility: the aquifer total compressibility (for the fluid and the solid) $c^A_t$ (in 1/Pa).
aquiferElevation: the elevation of the aquifer (in m).
aquiferThickness: the thickness of the aquifer (in m).
aquiferInnerRadius: the aquifer inner radius (in m).
aquiferAngle: the angle subtended by the aquifer boundary from the center of reservoir (in degrees, must be between 0 and 360).
allowAllPhasesIntoAquifer: flag controlling the behavior of the aquifer when there is flow from the reservoir to the aquifer. If the flag is equal to 1, all phases can flow into the aquifer. If the flag is equal to 0, only the water phase can flow into the aquifer. The default value of this optional parameter is 0.
pressureInfluenceFunctionName: the name of the table providing the dimensionless pressure as a function of dimensionless time. This table must be defined as a TableFunction in the Functions block of the XML file. If this optional parameter is omitted, a default pressure influence table is used.
setNames: the names of the face sets on which the aquifer boundary condition is applied.

Note

Following the GEOS convention, the z-coordinate is increasing upward. This convention must be taken into account when providing the aquiferElevation. In other words, the z-value is not a depth.

The full list of parameters is provided below:

XML Element: Aquifer

Name	Type	Default	Description
allowAllPhasesIntoAquifer	integer	0	Flag to allow all phases to flow into the aquifer. This flag only matters for the configuration in which flow is from reservoir to aquifer. - If the flag is equal to 1, then all phases, including non-aqueous phases, are allowed to flow into the aquifer. - If the flag is equal to 0, then only the water phase is allowed to flow into the aquifer. If you are in a configuration in which flow is from reservoir to aquifer and you expect non-aqueous phases to saturate the reservoir cells next to the aquifer, set this flag to 1. This keyword is ignored for single-phase flow simulations
aquiferAngle	real64	required	Angle subtended by the aquifer boundary from the center of the reservoir [degress]
aquiferElevation	real64	required	Aquifer elevation (positive going upward) [m]
aquiferInitialPressure	real64	required	Aquifer initial pressure [Pa]
aquiferInnerRadius	real64	required	Aquifer inner radius [m]
aquiferPermeability	real64	required	Aquifer permeability [m^2]
aquiferPorosity	real64	required	Aquifer porosity
aquiferThickness	real64	required	Aquifer thickness [m]
aquiferTotalCompressibility	real64	required	Aquifer total compressibility (rock and fluid) [Pa^-1]
aquiferWaterDensity	real64	required	Aquifer water density [kg.m^-3]
aquiferWaterPhaseComponentFraction	real64_array	{0}	Aquifer water phase component fraction. This keyword is ignored for single-phase flow simulations.
aquiferWaterPhaseComponentNames	string_array	{}	Aquifer water phase component names. This keyword is ignored for single-phase flow simulations.
aquiferWaterViscosity	real64	required	Aquifer water viscosity [Pa.s]
bcApplicationTableName	groupNameRef		Name of table that specifies the on/off application of the boundary condition.
beginTime	real64	-1e+99	Time at which the boundary condition will start being applied.
direction	R1Tensor	{0,0,0}	Direction to apply boundary condition to.
endTime	real64	1e+99	Time at which the boundary condition will stop being applied.
functionName	groupNameRef		Name of function that specifies variation of the boundary condition.
initialCondition	integer	0	Boundary condition is applied as an initial condition.
logLevel	integer	0	Log level
name	groupName	required	A name is required for any non-unique nodes
pressureInfluenceFunctionName	groupNameRef		Name of the table describing the pressure influence function . If not provided, we use a default pressure influence function
scale	real64	0	Scale factor for value of the boundary condition.
setNames	groupNameRef_array	required	Name of sets that boundary condition is applied to.

Examples

Setting up the Aquifer boundary condition requires two additional pieces of information in the XML input file: a set of faces to specify where the aquifer boundary conditions will apply, and an aquifer tag that specifies the physical characteristics of the aquifer and determines how the boundary condition is applied.

To specify a set of faces: on simple grids, in the Geometry block of the XML file, we can define a Box that selects and assigns a name to a set of faces. To be included in a set, the faces must be fully enclosed in the Box (all vertices of a face must be inside the box for the face to be included to the set). The name of this box is a user-defined string, and it will be used in the aquifer tag to locate the face set. Here is an example of XML code to create such a face set from a box:

<Geometry>
  ...
  <Box
    name="aquifer"
    xMin="{ 999.99, 199.99, 3.99 }"
    xMax="{ 1010.01, 201.01, 6.01 }"/>
  ...
</Geometry>

Note

This step captures faces, not cells. For now, the user must ensure that the box actually contains faces (GEOS will proceed even if the face set is empty).

For more complex meshes, sch as those imported using the VTKMesh, using a Box to perform a face selection is challenging. We recommend using a surface in the vtk file instead, which will be used to locate the face set.

To specify the aquifer characteristics: in the FieldSpecifications block of the XML file, we include an Aquifer tag. For single-phase flow, the aquifer definition looks like:

<FieldSpecifications>
  ...
  <Aquifer
    name="aquiferBC"
    aquiferPorosity="2e-1"
    aquiferPermeability="3e-13"
    aquiferInitialPressure="9e6"
    aquiferWaterViscosity="0.00089"
    aquiferWaterDensity="962.81"
    aquiferTotalCompressibility="1e-10"
    aquiferElevation="4"
    aquiferThickness="18"
    aquiferInnerRadius="2000"
    aquiferAngle="20"
    setNames="{ aquifer }"/>
  ...
</FieldSpecifications>

For compositional multiphase flow, the user must include additional parameters to specify the water composition. We have additional influx controls over the aquifer with allowAllPhasesIntoAquifer. This is illustrated below for the CO2-brine fluid model:

<FieldSpecifications>
  ...
  <Aquifer
    name="aquiferBC"
    aquiferPorosity="2e-1"
    aquiferPermeability="3e-13"
    aquiferInitialPressure="9e6"
    aquiferWaterViscosity="0.00089"
    aquiferWaterDensity="962.81"
    aquiferWaterPhaseComponentFraction="{ 0.0, 1.0 }"
    aquiferWaterPhaseComponentNames="{ co2, water }"
    aquiferTotalCompressibility="1e-10"
    aquiferElevation="4"
    aquiferThickness="18"
    aquiferInnerRadius="2000"
    aquiferAngle="20"
    allowAllPhasesIntoAquifer="1"
    setNames="{ aquifer }"/>
  ...
</FieldSpecifications>

Finally, for both single-phase and multiphase flow, if a pressureInfluenceFunctionName attribute is specified in the Aquifer tag, a TableFunction must be included in the Functions block of the XML file as follows:

<Functions>
  ...
  <TableFunction
    name="pressureInfluenceFunction"
    coordinates="{ 0.01, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5,
                   2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0,
                   50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 200.0, 800.0, 1600.0, 3200.0, 6400.0, 12800.0 }"
    values="{ 0.112, 0.229, 0.315, 0.376, 0.424, 0.469, 0.503, 0.564, 0.616, 0.659, 0.702, 0.735, 0.772, 0.802,
              0.927, 1.02, 1.101, 1.169, 1.275, 1.362, 1.436, 1.5, 1.556, 1.604, 1.651, 1.829, 1.96, 2.067, 2.147,
              2.282, 2.388, 2.476, 2.55, 2.615, 2.672, 2.723, 3.0537, 3.7468, 4.0934, 4.44, 4.7866, 5.1331 }"/>
  ...
</Functions>

Note

The values provided in the table above are the default values used internally in GEOS when the user does not specify the pressure influence function in the XML file.