Relative Permeability Functions#

TOUGH3 provides several relative permeability functions for two-phase or three-phase flow problems. The original TOUGH2 two-phase functions have been retained (from IRP = 1 to IRP = 8). The modified Brooks-Corey and van Genuchten models and two versions of the hysteresis model (regular and simple) are implemented (IRP = 10, 11, 12, and 13). The three-phase functions from TMVOC and ECO2M start at IRP = 31. For TMVOC, if one of the two-phase functions is chosen, the NAPL relative permeability is assumed to be zero. For ECO2M, the two-phase functions are only applicable when an aqueous phase and a single CO₂-rich phase are present. If one of the two-phase functions is chosen, the relative permeability of the CO₂-rich phase will be the same function of saturation, regardless whether the CO₂-rich phase is liquid or gas. The notation used below is:

$k_{rl}$: aqueous phase relative permeability
$k_{rg}$: gas (for all two-phase EOSs and TMVOC) or gaseous CO₂ (for ECO2M) phase relative permeability
$k_{rn}$: NAPL (for TMVOC) or liquid CO₂ (for ECO2M) phase relative permeability
$S_l$, $S_g$, and $S_n$ are the saturations of aqueous, gas (or gaseous CO₂) and NAPL (or liquid CO₂) phases, respectively

Linear Functions (IRP=1)#

\(k_{rl}\) increases linearly from 0 to 1 in range \(RP(1) \le S_l \le RP(3)\)
\(k_{rg}\) increases linearly from 0 to 1 in range \(RP(2) \le S_l \le RP(4)\)
Restrictions: \(RP(3) \gt RP(1)\), \(RP(4) \gt RP(2)\)

If $RP(5) \gt 0$, $k_{rn}$ increases linearly from 0 to 1 in the range $RP(1) \le S_n \le RP(3)$

If $RP(5) \gt 0$ and $RP(6) \gt 0$, $k_{rn}$ increases linearly from 0 to 1 in the range $RP(5) \le S_n \le RP(6)$

Restrictions: $RP(6) \lt RP(5)$

IRP=2#

\[k_{rl} = S_l^{RP(1)}\]

\[k_{rg} = 1\]

Corey’s Curves (IRP=3)#

After Corey, 1954.

\[k_{rl} = \hat{S}^4\]

\[k_{rg} = \left( 1 - \hat{S} \right)^2 \left( 1 - \hat{S}^2 \right)^2\]

where $\hat{S} = \frac{S_l - S_{lr}}{1 - S_{lr} - S_{gr}}$.

Parameters:

$RP(1) = S_{lr}$
$RP(2) = S_{gr}$

Restrictions: $RP(1) + RP(2) \lt 1$

Grant’s Curves (IRP=4)#

After Grant, 1977.

\[k_{rl} = \hat{S}^4\]

\[k_{rg} = 1 - k_{rl}\]

Parameters:

$RP(1) = S_{lr}$
$RP(2) = S_{gr}$

Restrictions: $RP(1) + RP(2) \lt 1$

All Phases Perfectly Mobile (IRP=5)#

$k_{rg} = k_{rg} = 1$ for all saturations.

No parameters.

Functions of Fatt and Klikoff (IRP=6)#

After Fatt & Klikoff, 1959.

\[k_{rl} = \left( S^\ast \right)^3\]

\[k_{rg} = \left( 1 - S^\ast \right)^3\]

where $S^\ast = \frac{S_l - S_{lr}}{1 - S_{lr}}$.

Parameters:

$RP(1) = S_{lr}$

Restriction: $RP(1) \lt 1$

van Genuchten-Mualem Model (IRP=7)#

After Mualem, 1976 and van VanGenuchten, 1980.

\[\begin{split}k_{rl} = \begin{cases} \sqrt{S^\ast} \left( 1 - \left( 1 - \left( S^\ast \right)^{\frac{1}{\lambda}} \right)^{\lambda} \right)^2 & \text{if $S_l \lt S_{ls}$} \\ 1 & \text{if $S_l \ge S_{ls}$} \\ \end{cases}\end{split}\]

Gas relative permeability can be chosen as one of the following two forms, the second of which is due to Corey, 1954.

\[\begin{split}k_{rg} = \begin{cases} 1 - k_{rl} & \text{if $S_{gr} = 0$} \\ \left( 1 - \hat{S} \right)^2 \left( 1 - \hat{S}^2 \right) & \text{if $S_{gr} \gt 0$} \\ \end{cases}\end{split}\]

subject to the restriction $0 \le k_{rl}$, $k_{rg} \le 1$.

Here,

\[S^\ast = \frac{S_l - S_{lr}}{S_{ls} - S_{lr}}\]

\[\hat{S} = \frac{S_l - S_{lr}}{1 - S_{lr} - S_{gr}}\]

Parameters:

$RP(1) = \lambda$
$RP(2) = S_{lr}$
$RP(3) = S_{ls}$
$RP(4) = S_{gr}$

Note

Parameter $\lambda$ is $m$ in van Genuchten’s notation, with $m = 1 - \frac{1}{n}$; parameter $n$ is often written as $\beta$.

Functions of Verma et al. (IRP=8)#

After Verma, 1985.

\[k_{rl} = \hat{S}^3\]

\[k_{rg} = A + B \hat{S} + C \hat{S}^2\]

where $\hat{S} = \frac{S_l - S_{lr}}{S_{ls} - S_{lr}}$

Parameters as measured by Verma et al. (1985) for steam-water flow in an unconsolidated sand:

Parameters:

$RP(1) = S_{lr} = 0.2$
$RP(2) = S_{ls} = 0.895$
$RP(3) = A = 1.259$
$RP(4) = B = -1.7615$
$RP(5) = C = 0.5089$

Modified Brooks-Corey Model (IRP=10)#

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 + \lambda}{\lambda}}\]

\[\begin{split}k_{rg} = \begin{cases} \left( 1 - S_{ek} \right)^2 \left( 1 - S_{ek}^{\frac{2 + \lambda}{\lambda}} \right) & \text{if $RP(3) = 0$} \\ 1 - k_{rl} & \text{if $RP(3) \ne 0$} \\ \end{cases}\end{split}\]

where

\[S_{ek} = \frac{S_l - S_{lrk}}{1 - S_{lrk} - S_{gr}}\]

Parameters:

$RP(1) = S_{lrk}$
$RP(2) = S_{gr}$
$RP(3) =$ flag to indicate which equation is used for $k_{rg}$

Modified van Genuchten Model (IRP=11)#

A modified version of the van Genuchten 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 van Genuchten model is invoked by setting both IRP and ICP to 11.

\[k_{rl} = S_{ekl}^{\gamma} S_{ekl}^{\left( 1 - \gamma \right) \eta} \left( 1 - \left( 1 - S_{ekl}^{\frac{1 - \gamma}{m}} \right)^m \right)^2\]

\[\begin{split}k_{rg} = \begin{cases} \left( 1 - S_{ekg} \right)^{\zeta} \left( 1 - S_{ekg}^{\frac{1}{m}} \right)^{2m} & \text{if $RP(3) = 0$} \\ 1 - k_{rl} & \text{if $RP(3) \ne 0$} \\ \end{cases}\end{split}\]

where

\[S_{ekl} = \frac{S_l - S_{lrk}}{1 - S_{lrk}}\]

\[S_{ekg} = \frac{S_l}{1 - S_{gr}}\]

Parameters:

$RP(1) = S_{lrk}$, if negative, $S_{lrk} = 0$ for calculating $k_{rg}$, and absolute value is used for calculating $k_{rl}$
$RP(2) = S_{gr}$, if negative, $S_{gr} = 0$ for calculating $k_{rl}$, and absolute value is used for calculating $k_{rg}$
$RP(3) =$ flag to indicate which equation is used for $k_{rg}$
$RP(4) = \eta$ (default = 1/2)
$RP(5) = \varepsilon_k$, use linear function between $k_{rl}$ ($S_e = 1 - \varepsilon_k$) and 1.0
$RP(6) = a_{fm}$, constant fracture-matrix interaction reduction factor, in combination with Active Fracture Model
$RP(7) = \zeta$ (default = 1/3)

Regular Hysteresis (IRP=12)#

The hysteretic form of the van Genuchten model (Lenhard & Parker, 1987, Parker & Lenhard, 1987) has been implemented. Details of the implementation are described in Doughty, 2013. The regular hysteresis model is invoked by setting both IRP and ICP to 12.

\[k_{rl} = \sqrt{\bar{S}_l} \left( 1 - \left( 1 - \frac{\bar{S}_{gt}}{1 - \bar{S}_l^{\Delta}} \right) \left( 1 - \left( \bar{S}_l + \bar{S}_{gt} \right)^{\frac{1}{m}}\right)^m - \frac{\bar{S}_{gt}}{1 - \bar{S}_l^{\Delta}} \left( 1 - \left( \bar{S}_l^{\Delta} \right)^{\frac{1}{m}} \right)^m \right)^2\]

\[k_{rg} = k_{rgmax} \left( 1 - \left( \bar{S}_l + \bar{S}_{gt} \right) \right)^{\gamma} \left( 1 - \left( \bar{S}_l + \bar{S}_{gt} \right)^{\frac{1}{m}} \right)^{2m}\]

where

\[\bar{S}_l = \frac{S_l - S_{lr}}{1 - S_{lr}}\]

\[\bar{S}_l^{\Delta} = \frac{S_l^{\Delta} - S_{lr}}{1 - S_{lr}}\]

\[\bar{S}_{gt} = \frac{S_{gr}^{\Delta} \left( S_l - S_l^{\Delta} \right)}{\left( 1 - S_{lr} \right) \left( 1 - S_l^{\Delta} - S_{gr}^{\Delta} \right)}\]

\[S_{gr}^{\Delta} = \frac{1}{\frac{1}{1 - S_l^{\Delta}} + \frac{1}{S_{gr, max}} - \frac{1}{1 - S_{lr}}}\]

$S_l^{\Delta}$ is the turning-point saturation, and $S_{gr}^{\Delta}$ is the residual gas saturation.

Parameters:

$RP(1) = m$, van Genuchten $m$ for liquid relative permeability (need not equal $CP(1)$ or $CP(6)); :math:`k_{rl}$ uses the same $m$ for drainage and imbibition.
$RP(2) = S_{lr}$, $k_{rl} (S_{lr}) = 0$, $k_{rg} (S_{lr}) = k_{rgmax}$. Must have $S_{lr} \gt S_{lmin}$ in capillary pressure ($CP(2)). :math:`S_{lr}$ is minimum saturation for transition to imbibition branch. For $S_l \lt S_{lr}$, curve stays on primary drainage branch even if $S_l$ increases.
$RP(3) = S_{grmax}$, maximum possible value of $S_{gr}^{\Delta}$. Note that the present version of the code requires that $S_{lr} + S_{grmax} \lt 1$, otherwise there will be saturations for which neither fluid phase is mobile, which the code cannot handle. Setting $S_{grmax} = 0$ effectively turns off hysteresis. As a special option, a constant, non-zero value of Sgr may be employed by setting $CP(10) \gt 1$ and making $RP(3)$ negative. The code will set $S_{gr}^{\Delta}$ = $-RP(3)$ for all grid blocks at all times.
$RP(4) = \gamma$, typical values 0.33 - 0.50
$RP(5) = k_{rgmax}$
$RP(6) =$ fitting parameter for $k_{rg}$ extension for $S_l \lt S_{lr}$ (only used when $k_{rgmax} \lt 1$); determines type of function for extension and slope of $k_{rg}$ at $S_l = 0$:
- ≤0: use cubic spline for $0 \lt S_l \lt S_{lr}$, with slope at $S_l = 0$ of $RP(6)$
- >0: use linear segment for $0 \lt S_l \lt RP(8) S_{lr}$ and cubic spline for $RP(8) S_{lr} \lt S_l \lt S_{lr}$, with slope at $S_l = 0$ of $-RP(6)$.
$RP(7) =$ numerical factor used for $k_{rl}$ extension to $S_l \gt S_l^\ast$. $RP(7)$ is the fraction of $S_l^\ast$ at which $k_{rl}$ curve departs from the original van Genuchten function. Recommended range of values: 0.95-0.97. For $RP(7) = 0$, $k_{rl} = 1$ for $S_l \gt S_l^\ast$ (not recommended).
$RP(8) =$ numerical factor used for linear $k_{rg}$ extension to $S_l \lt S_{lr}$ (only used when $k_{rgmax} \lt 1$). $RP(8)$ is the fraction of $S_{lr}$ at which the linear and cubic parts of the extensions are joined.
$RP(9) =$ flag to turn off hysteresis for $k_{rl}$ (no effect on $P_c$ and $k_{rg}$; to turn off hysteresis entirely, set $S_{grmax} = 0$ in $RP(3)$).
- 0: hysteresis is on for $k_{rl}$
- 1: hysteresis is off for $k_{rl}$ (force $k_{rl}$ to stay on primary drainage branch ($k_{rl}^d$) at all times)
$RP(10) = m_{gas}$, van Genuchten m for gas relative permeability (need not equal $CP(1)$ or $CP(6)$); $k_{rg}$ uses same $m_{gas}$ for drainage and imbibition. If zero or blank, use $RP(1)$ so that $m_{gas} = m$.

Simple Hysteresis (IRP=13)#

The regular hysteresis option (IRP = ICP = 12) provides a rigorous representation of hysteretic relative permeability and capillary pressure curves. However, it can significantly slow down TOUGH3 simulations, because small time steps are often required at turning points, when a grid block switches between drainage and imbibition, because the slopes of the characteristic curves are discontinuous. Moreover, several control parameters are needed, which generally must be determined by trial and error, for the code to run smoothly. An alternative means of capturing the essence of hysteresis, while maintaining continuous slopes and requiring no additional control parameters, is the simple hysteresis algorithm of Patterson & Falta, 2012, which is invoked with IRP = ICP = 13. Presently this option is only available when ECO2N is being used.

The Mualem, 1976 relative permeability model is used for the non-wetting phase:

\[k_{rn} = \sqrt{1 - \bar{S}_{wn}} \left( 1 - \bar{S}_{wn}^{\frac{1}{m}} \right)^{2m}\]

where

\[\bar{S}_{wn} = \frac{S_w - S_{wr}}{1 - S_{wr} - S_{nr}}\]

and $S_{wr}$ and $S_{nr}$ are residual saturations of the wetting and non-wetting phases, respectively. Hysteresis is implemented by considering $S_{nr}$ to be a variable, which is calculated from the maximum historical non-wetting phase saturation in a grid block, $S_{nmax}$. The user has the option to specify $S_{nr}$ as a linear function of the historical $S_{nmax}$:

\[S_{nr} = f_{snr} S_{nmax}\]

or $S_{nr}$ can be calculated using a modified form of the Land, 1968 relationship

\[S_{nr} = \frac{S_{nmax}}{1 + C S_{nmax}}\]

with

\[C = \frac{1}{S_{nrmax}} - \frac{1}{1 - S_{wr}}\]

where $f_{snr}$ and $S_{nrmax}$, the maximum residual non-wetting phase saturation, are user-specified material properties. $S_{nr}$ is calculated during every Newton-Raphson iteration. If $S_n$ drops below $S_{nr}$ by dissolution or compression, $S_{nmax}$ is recalculated as

\[S_{nmax} = \frac{S_n}{f_{snr}} \text{ or } S_{nmax} = \frac{S_n}{1 - C S_n}\]

Wetting-phase relative permeability (non-hysteretic) is from van Genuchten (1980)

\[k_{rw} = \sqrt{\bar{S}_w} \left( 1 - \left( 1 - \bar{S}_w^{\frac{1}{m}}\right)^m \right)^2\]

where

\[\bar{S}_w = \frac{S_w - S_{wr}}{S_{ws} - S_{wr}}\]

Parameters:

$RP(1) = m$ to use in $k_{rw}$
$RP(2) = S_{wr}$
$RP(3) = S_{ws}$ (recommend 1)
$RP(4)$
- <0: $= -f_{snr}$ in linear trapping model
- >0: $S_{nrmax}$ in Land trapping model
$RP(5) = m_{gas}$, $m$ to use in $k_{rn}$: if zero or blank, use $RP(1)$
$RP(6) =$ power to use in first term in $k_{rn}$ (default 1/2)
$RP(7)$
- =0: use $\left( 1 - \bar{S}_{wn} \right)$ in first term in $k_{rn}$ (Mualem, 1976)
- >0: use $S_g$ in first term in $k_{rn}$ (Charbeneau, 2007), so that $k_{rn}$ does not go to 1 when immobile liquid phase is present

All Phases Perfectly Mobile (IRP=31)#

$k_{rg} = k_{rl} = k_{rn} = 1$

No parameters.

Modified Version of Stone’s First Three-Phase Method (IRP=32)#

After Stone, 1970.

\[k_{rg} = \left( \frac{S_g - S_{gr}}{1 - S_{gr}} \right)^n\]

\[k_{rl} = \left( \frac{S_l - S_{lr}}{1 - S_{lr}} \right)^n\]

\[k_{rn} = \left( \frac{1 - S_g - S_l - S_{nr}}{1 - S_g - S_{lr} - S_{nr}} \right) \left( \frac{1 - S_{lr} - S_{nr}}{1 - S_l - S_{nr}} \right) \left( \frac{\left( 1 - S_g - S_{lr} - S_{nr} \right)\left( 1 - S_l \right)}{1 - S_{nr}} \right)^n\]

When $S_n = 1 - S_l - S_g - S_s$ is near irreducible liquid saturation, $S_{nr} \le S_n \le S_{nr} + 0.005$, liquid relative permeability is taken to be

\[k_{rn}^{'} = k_{rn} \frac{S_n - S_{nr}}{0.005}\]

Parameters:

$RP(1) = S_{lr}$
$RP(2) = S_{nr}$
$RP(3) = S_{gr}$
$RP(4) = n$

Three-Phase Functions of Parker et al. (IRP=33)#

After Parker et al., 1987.

\[k_{rg} = \sqrt{\bar{S}_g} \left( 1 - \left( \bar{S}_n \right)^{\frac{1}{m}} \right)^{2m}\]

\[k_{rl} = \sqrt{\bar{S}_l} \left( 1 - \left( 1 - \left( \bar{S}_l \right)^{\frac{1}{m}} \right)^m \right)^2\]

\[k_{rn} = \sqrt{\bar{S}_n - \bar{S}_l} \left( \left( 1 - \left( \bar{S}_l \right)^{\frac{1}{m}} \right)^m - \left( 1 - \left( \bar{S}_n \right)^{\frac{1}{m}} \right)^m \right)^2\]

with

\[m = 1 - \frac{1}{n}\]

\[\bar{S}_g = \frac{S_g}{1 - S_m}\]

\[\bar{S}_l = \frac{S_l - S_m}{1 - S_m}\]

\[\bar{S}_n = \frac{S_l + S_n - S_m}{1 - S_m}\]

where $k_{rg}$, $k_{rl}$, and $k_{rn}$ are limited to values between 0 and 1.

Parameters:

$RP(1) = S_m$
$RP(2) = n$

IRP=34#

Same as Modified Version of Stone’s First Three-Phase Method (IRP=32), except that

\[k_{rg} = 1 - \left( \frac{S_n + S_l - S_{lr}}{1 - S_{lr}} \right)^n\]

Power Law (IRP=35)#

Phases $\beta = l, n, g$.

\[k_{r \beta} = \left( \frac{S_{\beta} - S_{\beta r}}{1 - S_{\beta r}} \right)^n\]

Parameters:

$RP(1) = S_{lr}$
$RP(2) = S_{nr}$
$RP(3) = S_{gr}$
$RP(4) = n$

Functions Used by Faust (1985) for Two-Phase Buckley-Leverett Problem#

After Faust, 1985.

\[k_{rl} = \frac{\left( S_l - 0.16 \right)^2}{0.64}\]

\[k_{rg} = 0\]

\[k_{rn} = \frac{\left( 0.8 - S_l \right)^2}{0.64}\]

where $k_{rl}$ and $k_{rn}$ are limited to values between 0 and 1.

No parameters.

IRP=37#

Same are Modified Version of Stone’s First Three-Phase Method (IRP=32), except a correction factor is applied to $k_{rn}$ such as to make $k_{rn}$ equal to $k_{rg}$ for two-phase conditions with the same aqueous phase saturation.

Custom#

If the user wishes to employ other relative permeability relationships, these need to be programmed into subroutine RELP in module Utility.f90. The routine has the following structure:

SUBROUTINE RELP(SATU,RELPERM,NNPH,NMAT,USRX)
    ...
    RELP_FUNCTION: SELECT CASE (IRP(NMAT))
            CASE (1)
                    CALL RELP_LINEAR(...)
            CASE (2)
                    ...
            ...
            ...
            END SELECT RELP_FUNCTION

END

To code an additional relative permeability function, the user needs to insert a code segment analogous to that shown above, beginning with a CASE option which would be identical to IRP and calls a subroutine for the additional relative permeability function.