Compositional multiphase fluid model

Overview

The compositional fluid models implemented in this framework provide a robust and flexible approach for modelling the thermodynamic behaviour of multicomponent, multiphase fluid mixtures. These models are intended to support complex subsurface flow simulations, including primary hydrocarbon recovery, miscible gas injection, carbon sequestration, and geothermal energy extraction.

Fluid composition and phase partitioning

In this model, the fluid is described by $N_c$ components, with $z_c$ being the total mole fraction of component $c$ . The fluid can partition into up to three phases: a liquid phase (denoted $\ell$ ), a vapour phase (denoted $v$ ), and potentially an aqueous phase (denoted $a$ ).

Therefore, by taking into account the molar phase component fractions (which is the fraction of the molar mass of phase $p$ represented by component $c$ ), the following partition matrix establishes the component distribution within a two-phase (liquid-vapour) system:

$\begin{bmatrix} x_1 & x_2 & x_3 & \cdots & x_{N_c} \\ y_1 & y_2 & y_3 & \cdots & y_{N_c} \end{bmatrix}$

where $x_c$ is the mole fraction of component $c$ in the liquid phase and $y_c$ is the mole fraction of component $c$ in the vapour phase.

If a third aqueous phase is also present, with component mole fractions denoted by $w_c$ , this partition matrix naturally extends to three rows to fully describe the component distribution across all three phases:

$\begin{bmatrix} x_1 & x_2 & x_3 & \cdots & x_{N_c} \\ y_1 & y_2 & y_3 & \cdots & y_{N_c} \\ w_1 & w_2 & w_3 & \cdots & w_{N_c} \end{bmatrix}$

Structure of the fluid property calculations

The evaluation of fluid properties follows a sequential thermodynamic structure designed to ensure consistency across all phase properties. At a given pressure, temperature, and overall composition, the framework executes the following overarching evaluations:

State Conversion: If the system is formulated using mass variables, the overall mass fractions are converted into overall mole fractions.
Phase Split (Flash Calculation): The thermodynamic stability of the mixture is evaluated. If unstable, a flash calculation splits the mixture into separate phases (e.g., liquid, vapour, and potentially an aqueous phase), determining the phase fractions and the composition of each individual phase.
Phase Properties: Using the computed phase compositions, the model evaluates properties for each specific phase:
- Density: Evaluated typically via a cubic Equation of State (EoS) combined with volume shift corrections.
- Viscosity: Evaluated using established compositional correlations or empirical models.
- Thermal Properties: If the simulation is non-isothermal, the phase enthalpies and internal energies are computed using ideal gas heat capacities and EoS departure functions.
Total Fluid Properties: The individual phase properties are homogenised to yield total fluid properties (e.g., total density). Derivatives with respect to pressure, temperature, and composition are rigorously propagated using the chain rule.

Structure of this Documentation

This documentation is divided into the following sections:

Phase Split: Details the thermodynamic stability and flash calculations (Negative Flash and K-Value Flash), including the underlying Equations of State.
Density: Describes the models used to compute molar and mass phase densities.
Viscosity: Details the models used to evaluate the resistance to flow for each phase.
Enthalpy: Explains the calculation of thermal properties.
Miscellaneous: Outlines internal unit conversions and mass-to-molar variable handling.
References: Provides the bibliography for the scientific correlations implemented in the code.

Phase split

The primary purpose of the phase split (or flash) calculation is to determine the number of fluid phases present at equilibrium, the fraction of the total mixture occupying each phase, and the detailed molar composition of these phases. This is a fundamental requirement for compositional modelling, as the physical properties (density, viscosity) of a fluid depend heavily on its phase state.

Phase split equations

The phase split relies on two fundamental physical principles: the conservation of mass (material balance) and thermodynamic equilibrium.

Material Balance

For a two-phase system (liquid and vapour), the overall mole fraction of a component, $z_i$ , must equal the sum of its molar contributions in the liquid phase, $x_i$ , and the vapour phase, $y_i$ . Introducing the vapour phase fraction, $V$ , and the equilibrium ratio (K-value), $K_i = y_i / x_i$ , the material balance is expressed by the Rachford-Rice equation (Rachford and Rice, 1952):

$\sum_{i=1}^{N_c} \frac{z_i (K_i - 1)}{1 + V (K_i - 1)} = 0$

The root of this equation provides the vapour fraction $V$ . The phase compositions are subsequently derived via:

$x_i = \frac{z_i}{1 + V(K_i - 1)}$

$y_i = K_i x_i$

Thermodynamic Equilibrium

At thermodynamic equilibrium, the chemical potential, or equivalently the fugacity, of each component must be identical across all existing phases:

$f_i^L(P, T, x) = f_i^V(P, T, y)$

Using the fugacity coefficient, $\phi_i$ , where $f_i = x_i \phi_i P$ , the equilibrium condition dictates the K-values:

$K_i = \frac{\phi_i^L}{\phi_i^V}$

The phase split models in this framework either rigorously solve for these fugacity coefficients using an Equation of State or use simplified, pre-tabulated K-values.

Equations of state

An Equation of State (EoS) mathematically relates the pressure, volume, temperature, and composition of a fluid. In compositional modelling, the EoS serves a dual purpose: 1. It calculates the fugacity coefficients, $\phi_i$ , required to establish thermodynamic phase equilibrium. 2. It calculates the compressibility factor, $Z$ , required to determine the phase density.

The fugacity coefficient of a component is derived from the exact thermodynamic relationship involving the integration of the EoS volume departure with respect to pressure.

Cubic equation of state

The Peng-Robinson (PR) (Peng and Robinson, 1976; Robinson and Peng, 1978) and Soave-Redlich-Kwong (SRK) (Soave, 1972) equations of state are widely used for hydrocarbon systems encompassing natural gas, gas condensates, and volatile oils. The fundamental pressure-volume-temperature (PVT) relationships for these equations are defined as:

Peng-Robinson (PR):

$P = \frac{RT}{v - b} - \frac{a(T)}{v(v + b) + b(v - b)}$
Soave-Redlich-Kwong (SRK):

$P = \frac{RT}{v - b} - \frac{a(T)}{v(v + b)}$

In these PVT relationships, several intermediate parameters hold specific physical meanings:

$v$ represents the molar volume of the fluid phase.
$R$ is the universal gas constant, serving as the fundamental scaling factor relating macroscopic thermodynamic state variables to the molar energy.
$a(T)$ is the temperature-dependent attractive parameter, representing the intermolecular attractive forces between the fluid molecules.
$b$ is the repulsive co-volume parameter, representing the effective physical volume occupied by the fluid molecules themselves.

To efficiently solve these equations in a multicomponent mixture, they are rewritten into a generalized cubic equation for the compressibility factor, $Z$ :

$Z^3 + E_2 Z^2 + E_1 Z + E_0 = 0$

where the coefficients are defined by the dimensionless mixture parameters $A$ and $B$ :

$E_2 = (\delta_1 + \delta_2 - 1)B - 1$

$E_1 = A + \delta_1 \delta_2 B^2 - (\delta_1 + \delta_2)B(B + 1)$

$E_0 = -[AB + \delta_1 \delta_2 B^2(B + 1)]$

The constants $\delta_1$ and $\delta_2$ determine the specific EoS used:

Peng-Robinson (PR): $\delta_1 = 1 + \sqrt{2}$ , $\delta_2 = 1 - \sqrt{2}$
Soave-Redlich-Kwong (SRK): $\delta_1 = 0$ , $\delta_2 = 1$

The mixture parameters $A$ and $B$ are calculated using standard Van der Waals mixing rules over the phase mole fractions, $x_i$ :

$A = \sum_i \sum_j x_i x_j (1 - k_{ij}) \sqrt{A_i A_j}$

$B = \sum_i x_i B_i$

where $k_{ij}$ are the binary interaction coefficients (BICs). The pure component parameters $A_i$ and $B_i$ are defined as:

$A_i = \Omega_a \frac{P_{r,i}}{T_{r,i}^2} \alpha_i$

$B_i = \Omega_b \frac{P_{r,i}}{T_{r,i}}$

where $P_{r,i} = P/P_{c,i}$ and $T_{r,i} = T/T_{c,i}$ are the reduced pressure and temperature of component $i$ .

The dimensionless constants $\Omega_a$ and $\Omega_b$ are derived mathematically by applying the critical point constraints (where the first and second pressure-volume derivatives vanish) to dictate the weight of the attractive and repulsive forces. The exact values implemented in the code are:

Peng-Robinson: $\Omega_a = 0.457235529$ , $\Omega_b = 0.077796074$
Soave-Redlich-Kwong: $\Omega_a = 0.42748$ , $\Omega_b = 0.08664$

The alpha function, $\alpha_i$ , accounts for the temperature dependence of the attractive parameter and relies on the acentric factor, $\omega_i$ , of each component. For both equations, it takes the form:

$\alpha_i = \left( 1 + m_i \left( 1 - \sqrt{T_{r,i}} \right) \right)^2$

The parameter $m_i$ is computed differently depending on the chosen equation of state:

Peng-Robinson: If $\omega_i < 0.49$ :

$m_i = 0.37464 + 1.54226 \omega_i - 0.26992 \omega_i^2$

If $\omega_i \ge 0.49$ :

$m_i = 0.3796 + 1.485 \omega_i - 0.164423 \omega_i^2 + 0.016666 \omega_i^3$
Soave-Redlich-Kwong:

$m_i = 0.480 + 1.574 \omega_i - 0.176 \omega_i^2$

Finally, the fugacity coefficient, $\phi_i$ , for each component is calculated using the derived roots for $Z$ :

$\ln \phi_i = (Z - 1) \frac{B_i}{B} - F - G \left( 2 \sum_j x_j (1 - k_{ij}) \sqrt{A_i A_j} - A \frac{B_i}{B} \right) E$

where the intermediate terms are defined as:

$E = \ln \left( \frac{Z + \delta_1 B}{Z + \delta_2 B} \right), \quad F = \ln(Z - B), \quad G = \frac{1}{(\delta_1 - \delta_2)B}$

Soreide-Whitson equation of state

The Soreide-Whitson model represents a specialized modification of the Peng-Robinson EoS tailored specifically for aqueous-hydrocarbon mixtures, such as those encountered in carbon sequestration (CO2-brine systems) or water-flooding in oil reservoirs. It adapts the standard cubic EoS by capturing the “salting-out” effect, where increased brine salinity reduces the solubility of hydrocarbons and CO2 in the aqueous phase. It is highly recommended whenever accurate modelling of gas dissolution in saline water is critical.

This model introduces two primary modifications to the standard Peng-Robinson formulation, applying special treatment exclusively to the water component and its interactions within the aqueous phase and follows closely the description given by Soreide and Whitson (1992):

First, the standard alpha function, $\alpha_i$ , is replaced specifically for the water component by a salinity-dependent formulation. If $c_{sw}$ represents the brine salinity, the water alpha function is evaluated as:

$\alpha_w = \left[ 1 + 0.4530 \left( 1 - T_{r,w} \left( 1 - 0.0103 c_{sw}^{1.1} \right) \right) + 0.0034 \left( T_{r,w}^{-3} - 1 \right) \right]^2$

where $T_{r,w}$ is the reduced temperature of water.

Second, unlike standard equations of state that use constant binary interaction coefficients (BICs), this model calculates dynamic aqueous phase BICs, $k_{ij}^{AQ}$ , exclusively for pairs where one component is water and the other is a hydrocarbon or specific non-hydrocarbon.

For water-hydrocarbon pairs, the interaction coefficient is modeled as a polynomial function of the hydrocarbon’s reduced temperature, $T_{r,i}$ , its acentric factor, $\omega_i$ , and the brine salinity:

$k_{ij}^{AQ} = A_0 (1 + \alpha_0 c_{sw}) + A_1 T_{r,i} (1 + \alpha_1 c_{sw}) + A_2 T_{r,i}^2 (1 + \alpha_2 c_{sw})$

The constants for this polynomial depend on the acentric factor. A special treatment is applied for methane, identified by an acentric factor $\omega_i \le 0.02$ (methane’s acentric factor is approximately 0.0108). To achieve higher accuracy for methane-brine systems, a specifically tuned set of empirical constants is used instead of the generalized hydrocarbon constants.

For standard hydrocarbons ( $\omega_i > 0.02$ ):

$A_0 = 1.1120 - 1.7369 \omega_i^{-0.1}$
$A_1 = 1.1001 + 0.8360 \omega_i$
$A_2 = -0.15742 - 1.0988 \omega_i$
$\alpha_0 = 0.017407$
$\alpha_1 = 0.033516$
$\alpha_2 = 0.011478$

For methane ( $\omega_i \le 0.02$ ):

$A_0 = 1.616705 - 1.884534 \omega_i^{-0.1}$
$A_1 = 0.8157014 + 0.8723632 \omega_i$
$A_2 = -0.0887821 - 0.8767864 \omega_i$
$\alpha_0 = 0.09988448$
$\alpha_1 = 0.1485516$
$\alpha_2 = 0.2111324$

For specific non-hydrocarbon components interacting with water, the dynamic BICs are calculated using explicitly tuned correlations based on the component type:

For carbon dioxide (CO2):

$k_{ij}^{AQ} = A_0 a_0 + A_1 a_1 T_{r,i} + A_2 a_2$

where $A_0 = -0.31092$ , $A_1 = 0.23580$ , $A_2 = -21.2566$ , $a_0 = 1 + 0.15587 c_{sw}^{0.7505}$ , $a_1 = 1 + 0.17837 c_{sw}^{0.979}$ , and $a_2 = \exp(-6.7222 T_{r,i} - c_{sw})$ .

For nitrogen (N2):

$k_{ij}^{AQ} = A_0 a_0 + A_1 a_1 T_{r,i}$

where $A_0 = -1.70235$ , $A_1 = 0.44338$ , $a_0 = 1 + 0.025587 c_{sw}^{0.75}$ , and $a_1 = 1 + 0.08126 c_{sw}^{0.75}$ .

For hydrogen sulfide (H2S):

$k_{ij}^{AQ} = A_0 + A_1 T_{r,i}$

where $A_0 = -0.20441$ and $A_1 = 0.23426$ .

For hydrogen (H2) the correlation due to Chabab et al. (2024) is used:

$k_{ij}^{AQ} = D_0 (1 + a_0 c_{sw}) + D_1 T_{r,i} (1 + a_1 c_{sw}) + D_2 \exp(D_3 T_{r,i})$

where $D_0 = -2.11917$ , $D_1 = 0.14888$ , $D_2 = -13.01835$ , $D_3 = -0.43946$ , $a_0 = -0.0226322$ , and $a_1 = -0.0044736$ .

Component naming requirements

Because the Soreide-Whitson EoS applies highly specific correlations based on the exact type of component interacting with water, it relies on the component’s assigned name to identify its type. This name matching is case-insensitive.

To trigger the correct dynamic binary interaction coefficients, you must use the following exact names for the special components:

Water: H2O
Carbon dioxide: CO2
Nitrogen: N2
Hydrogen sulfide: H2S
Hydrogen: H2

If a component’s name does not match any of these predefined strings (for example, if you name it “CarbonDioxide”, “C1”, or “Methane”), the model will automatically default to categorising it as a generic hydrocarbon (hc). While this is the intended and correct behavior for actual hydrocarbons (where methane is safely identified via its acentric factor), using an unrecognized name for a special non-hydrocarbon (like CO2) will cause the simulator to mistakenly apply the generic hydrocarbon-water correlation instead of its highly tuned specific correlation.

Model parameters

Equations of state are assigned per-phase within the specific fluid model XML block. The equationsOfState attribute takes a list specifying the EoS type for each phase corresponding to the order defined in the phaseNames list. Standard component properties necessary for the EoS calculations must also be provided in matching component order.

equationsOfState: List specifying the EoS type (e.g., PengRobinson, SoaveRedlichKwong, SoreideWhitson) corresponding to each phase.
componentCriticalPressure: Critical pressure of each component. Unit: [Pa].
componentCriticalTemperature: Critical temperature of each component. Unit: [K].
componentAcentricFactor: Acentric factor of each component. Unit: [dimensionless].
componentVolumeShift: Volume shift parameter for each component (optional). Unit: [dimensionless].
componentBinaryCoeff: Flattened $N \times N$ matrix of binary interaction coefficients (optional). Unit: [dimensionless].
salinity: Brine salinity utilized by the Soreide-Whitson aqueous model (optional, defaults to 0.0). Unit: [mol/kg].

Negative two-phase flash

The negative flash is a robust, rigorous approach to solving the phase equilibrium problem using an Equation of State (EoS).

Stability analysis

Before executing a flash calculation, the model must ascertain if the homogeneous single-phase mixture is thermodynamically unstable at the specified pressure $P$ and temperature $T$ . This is performed via a tangent plane distance (TPD) analysis developed by Michelsen (Michelsen, 1982a; Michelsen, 1982b). The TPD assesses whether the Gibbs free energy of the mixture can be minimised by splitting into multiple phases. A phase with feed composition $z$ is stable if and only if:

$g(\mathbf{y}) = \sum_{i=1}^{N_c} y_i \left( \ln y_i + \ln \phi_i(\mathbf{y}) - \ln z_i - \ln \phi_i(\mathbf{z}) \right) \ge 0$

for any permissible trial composition $\mathbf{y}$ , where $\phi_i$ denotes the fugacity coefficient of component $i$ .

To determine stability, testing is initiated from basic starting points using Wilson’s K-values (Wilson, 1969; Whitson and Torp, 1981) to generate a lighter trial mixture ( $y_i = z_i / K_i$ ) and a heavier trial mixture ( $y_i = z_i K_i$ ). Wilson’s K-values are defined as:

$K_i = \frac{P_{ci}}{P} \exp \left( 5.37 (1 + \omega_i) \left( 1 - \frac{T_{ci}}{T} \right) \right)$

where $P_{ci}$ , $T_{ci}$ , and $\omega_i$ are the critical pressure, critical temperature, and acentric factor of component $i$ , respectively.

The stationarity condition, $\ln Y_i + \ln \phi_i(\mathbf{y}) - h_i = 0$ (where $h_i = \ln z_i + \ln \phi_i(\mathbf{z})$ and $Y_i$ are unnormalized trial phase moles), is solved using successive substitution. If both initial states converge to a solution where $g(\mathbf{y}) \ge 0$ , the mixture is stable; otherwise, it is unstable.

Phase labeling

If the stability test confirms a single stable phase, the two-phase flash calculation is bypassed. To appropriately label the single phase as liquid or vapour for relative permeability evaluations, the model utilises a “Li-temperature” correlation. This pseudo-critical temperature is weighted by the critical volumes ( $V_{ci}$ ):

$T_{cp} = \frac{\sum_{i=1}^{N_c} T_{ci} V_{ci} z_i}{\sum_{i=1}^{N_c} V_{ci} z_i}$

If $T_{cp} < T$ , the mixture is logically categorised as a vapour; otherwise, it is labelled as a liquid.

Flash calculation

If the mixture is deemed unstable, the system must determine the phase split by ensuring thermodynamic equilibrium, meaning the fugacities in the liquid and vapour phases must be equal ( $\phi_{iL} = \phi_{iV}$ ).

Based on the material balance and the equilibrium constants ( $K_i = y_i/x_i$ ), the mole fractions in the liquid ( $x_i$ ) and vapour ( $y_i$ ) phases are given by:

$x_i = \frac{z_i}{1 + (K_i - 1)V}, \quad y_i = \frac{K_i z_i}{1 + (K_i - 1)V}$

The vapour molar fraction, $V$ , is determined by solving the Rachford-Rice equation alongside the EoS:

$F(V) = \sum_{i=1}^{N_c} (x_i - y_i) = \sum_{i=1}^{N_c} \frac{z_i(1 - K_i)}{1 + (K_i - 1)V} = 0$

The flash algorithm:

Initialize: An initial set of K-values is chosen using Wilson’s formula.
Solve for V: The Rachford-Rice equation is solved for the vapour fraction $V$ (using successive substitution, followed by Newton iterations).
Compute compositions: The corresponding liquid ( $x_i$ ) and vapour ( $y_i$ ) mole fractions are computed.
Evaluate EoS: These compositions are used to calculate the component fugacities $\phi_{iL}$ and $\phi_{iV}$ via the Equation of State.
Check convergence: Convergence is reached when $\sum_{i=1}^{N_c} (\phi_{iL} - \phi_{iV})^2 < \epsilon$ .
Update K-values: If not converged, the algorithm employs successive substitution iterations, constantly updating the K-values using fugacity coefficients derived from the EoS:

$K_i^{(new)} = K_i^{(old)} \frac{\phi_{iL}}{\phi_{iV}} = K_i^{(old)} \exp \left( \ln \phi_i^L - \ln \phi_i^V \right)$

The term “negative” flash arises because the algorithm temporarily permits the vapour fraction, $V$ , to converge to values slightly outside the physically meaningful bounds of $[0, 1]$ . This mathematical relaxation prevents the solver from getting artificially trapped at phase boundaries, vastly improving convergence robustness near the critical point or saturation envelopes. Upon convergence, if $V$ is negative or greater than unity, the system truncates to a single-phase solution.

Recommendations

The negative two-phase flash is the recommended standard for rigorous compositional modelling of miscible gas injection, primary depletion of volatile oils, and scenarios where fluid compositions change drastically over time. An extended three-phase variant is available for immiscible water systems.

Parameters

stabilityThreshold: Tangent plane distance below which a mixture is unstable (default: -1.0e-8). Unit: [dimensionless].
stabilityTolerance: Tolerance for stationarity in the stability test. Unit: [dimensionless].
stabilityMaxIterations: Maximum successive substitution steps for stability analysis. Unit: [dimensionless].
flashTolerance: Convergence tolerance for the fugacity ratio error. Unit: [dimensionless].
flashMaxIterations: Maximum successive substitution steps for the flash solve. Unit: [dimensionless].

K-value flash

The K-value flash model offers a computationally simplified alternative to the rigorous equation of state flash calculation. Rather than iteratively evaluating complex equations of state to find fugacity coefficients, this model determines phase splits using pre-tabulated K-values.

The model supports both two-phase and three-phase variants. The two-phase flash calculates equilibrium partitioning between a liquid and a vapour phase. The three-phase flash introduces an additional aqueous phase, allowing for the predefined partitioning of components into water alongside the hydrocarbon phases.

Methodology

In this formulation, the equilibrium ratio for each component in the secondary phases relative to the primary liquid phase, $K_i = y_i / x_i$ , is treated strictly as a predefined function of pressure and temperature.

Users input these functions as either tables or symbolic functions. To optimize performance during the simulation, the model discretizes these functions during initialization to construct a multi-dimensional hypercube of K-values spanning the phase, component, pressure, and temperature dimensions.

The pressure and temperature coordinates for this hypercube grid can be explicitly chosen and specified by the user. If they are omitted, the model automatically inspects the provided table functions, extracts all unique coordinate points from the underlying data, and constructs the grid dynamically. If only a single point is found, an artificial offset is added to permit interpolation.

Equilibrium calculation

During the simulation step, the model performs bilinear interpolation on the pre-built hypercube to evaluate the K-values and their partial derivatives at the current pressure and temperature.

Because the K-values are explicitly known, the iterative fugacity updates are bypassed. Instead, the vapour fraction, $V$ , is found by directly solving the Rachford-Rice equation (Rachford and Rice, 1952) using successive substitution and Newton iterations:

$\sum_{i=1}^{N_c} \frac{z_i (K_i - 1)}{1 + V(K_i - 1)} = 0$

Once the vapour fraction is determined, the phase compositions are directly evaluated:

$x_i = \frac{z_i}{1 + V(K_i - 1)}, \quad y_i = K_i x_i$

This direct solution approach intrinsically supports a “negative” flash. The mathematical solver allows the vapour fraction, $V$ , to temporarily converge to values outside the physically meaningful bounds of $[0, 1]$ . If the resulting $V \le 0$ or $V \ge 1$ , the system cleanly truncates the compositions to represent a single-phase liquid or single-phase vapour solution, respectively.

Recommendations and pitfalls

The K-value flash is highly computationally efficient and is recommended for thermal recovery processes (like steam injection) or relatively stable black-oil systems where fluid compositions remain within a narrow, predictable envelope.

Pitfalls: In reality, thermodynamic K-values are heavily dependent on the overall mixture composition, not just pressure and temperature. Relying on tabulated K-values implies assuming a fixed compositional path. Care should be taken when using this model for highly variable compositional processes such as miscible gas injection (e.g., CO2 or hydrocarbon gas flooding), or strong depletion of gas condensates, where composition dependency is critical to accurate phase behaviour.

Model parameters

The K-value flash fluid models are assigned using specific XML blocks that pair the flash model with viscosity and density models.

kValueTables: List of function names linking to the pre-tabulated K-values. For an $N$ -phase system, this requires $(N-1) \times N_c$ function names.
pressureCoordinates: List of pressure values for interpolation grid generation (optional). Unit: [Pa].
temperatureCoordinates: List of temperature values for interpolation grid generation (optional). Unit: [K].

Immiscible water flash

The immiscible water flash is an efficient three-phase equilibrium model designed for systems containing an aqueous phase alongside hydrocarbon liquid and vapour phases. Instead of solving a fully coupled three-phase equation of state problem, this model vastly simplifies the phase split by assuming water forms a strictly pure aqueous phase and is completely immiscible in the hydrocarbon phases, while hydrocarbons are entirely insoluble in the aqueous phase.

Note

When using the immiscibile water flash, the water component still needs to be explicitly specified as one of the components with the name H2O.

Methodology

In this formulation, the mixture composition is explicitly separated into a water component and a hydrocarbon mixture before any thermodynamic equilibrium calculations occur. The water component is identified using its assigned component name (H2O).

Given a total feed mole fraction $z_i$ for each component, the mole fraction of the aqueous phase, $V_{aq}$ , is determined to be exactly equal to the total feed mole fraction of the water component:

$V_{aq} = z_{w}$

The composition of the aqueous phase is fixed at 100% water. The remaining feed defines the total hydrocarbon mole fraction, $z_{hc}$ :

$z_{hc} = 1 - z_{w}$

The feed composition of the remaining hydrocarbon mixture is then normalised:

$z_{i,hc} = \frac{z_i}{z_{hc}}$

for all non-water components, while the normalised hydrocarbon feed fraction for water is set to zero.

Equilibrium calculation

During the simulation step, if the total hydrocarbon fraction $z_{hc}$ is negligibly small, the model bypasses the flash and assigns the system as a single-phase aqueous fluid.

If hydrocarbons are present, the model performs a rigorous two-phase negative flash exclusively on the normalised hydrocarbon mixture ( $z_{i,hc}$ ) using the chosen equations of state for the liquid and vapour phases.

This embedded flash calculation determines the intermediate hydrocarbon vapour fraction, $V_{V,hc}$ , and the corresponding compositions of the hydrocarbon liquid ( $x_i$ ) and hydrocarbon vapour ( $y_i$ ) phases using the standard Rachford-Rice equation and iterative fugacity updates.

Once the two-phase hydrocarbon flash converges, the overall phase fractions for the entire three-phase system are scaled back using the total hydrocarbon mole fraction:

$V_{V} = z_{hc} V_{V,hc}$

$V_{L} = z_{hc} (1 - V_{V,hc})$

The derivatives of the phase fractions and compositions with respect to pressure, temperature, and total composition are analytically transformed using the chain rule to account for this scaling and normalisation.

Recommendations and pitfalls

The immiscible water flash is highly computationally efficient and is recommended for multi-phase reservoir simulations where water solubility in the hydrocarbon phases (and vice versa) is negligible. By avoiding a full three-phase fugacity solver, it drastically reduces simulation time and improves non-linear solver stability.

Pitfalls: Because this model strictly enforces total immiscibility, it cannot accurately model scenarios where mutual solubility is physically significant, such as high-pressure CO2 injection into saline aquifers (where CO2 dissolves heavily into brine, and water vaporises into the gas stream). For such complex mutual solubility problems, using a two-phase flash paired with the Soreide-Whitson equation of state might be required.

Model parameters

The immiscible water flash fluid models are assigned using specific XML blocks that pair the three-phase flash model with distinct viscosity and density models for the hydrocarbon and aqueous phases.

equationsOfState: List specifying the EoS type for each phase (liquid, vapour, aqueous). Note that the aqueous EoS entry is largely ignored by the pure-water density model but must be provided for list alignment.

Density

The accurate determination of phase density is critical for modelling fluid flow, as it directly impacts volumetric calculations, gravity drainage, and phase mobilities.

Compositional density model

The compositional density model first calculates the uncorrected molar volume, $v_{EOS}$ , of a phase using the compressibility factor, $Z$ , obtained from the cubic Equation of State:

$v_{EOS} = \frac{Z R T}{P}$

To improve the accuracy of liquid phase densities, a Peneloux volume shift s described by Péneloux et al. (1982). correction is applied. The corrected molar volume, $v$ , is calculated by subtracting a composition-weighted sum of component-specific dimensional volume shifts, $c_i$ , from the EOS molar volume:

$v = v_{EOS} - \sum_{i=1}^{N_c} x_i c_i$

where $x_i$ is the mole fraction of component $i$ in the phase.

The dimensional volume shifts, $c_i$ , are computed from the user-provided non-dimensional volume shift parameters, $s_i$ (specified as componentVolumeShift), along with the critical temperature, $T_{c,i}$ , and critical pressure, $P_{c,i}$ , of each component:

$c_i = s_i \Omega_B R \frac{T_{c,i}}{P_{c,i}}$

The EOS-specific constant $\Omega_B$ takes the value $0.077796074$ for the Peng-Robinson equation of state and $0.08664$ for the Soave-Redlich-Kwong equation of state.

The molar density of the phase, $\rho_{molar}$ , is the inverse of the corrected molar volume:

$\rho_{molar} = \frac{1}{v}$

Finally, the mass density of the phase, $\rho_{mass}$ , is obtained by multiplying the molar density by the apparent molecular weight of the mixture:

$\rho_{mass} = \rho_{molar} \sum_{i=1}^{N_c} x_i MW_i$

where $MW_i$ is the molecular weight of component $i$ . This model is recommended for all hydrocarbon phases.

Parameters

componentVolumeShift: Component-specific volume shift parameters. Unit: [dimensionless].

Phillips brine density model

For aqueous phases containing dissolved salts, cubic equations of state often struggle to capture complex ionic interactions. To address this, the Phillips brine density model utilizes the empirical correlation developed by Phillips et al. (1981) to pre-compute a 2D lookup table of volume shifts specifically for the water component as a function of pressure and temperature.

First, the apparent molar weight of the brine, $MW_{brine}$ , is calculated from the salinity, $S$ (molality), the salt molar weight, $MW_{salt}$ , and the water molar weight, $MW_{H2O}$ :

$MW_{brine} = \frac{MW_{H2O}}{1 + S\left(MW_{H2O} - MW_{salt}\right)}$

The target mass density of the brine, $\rho_{brine}$ , is then evaluated. The pressure must be in Pascal and must be less than $5 \times 10^7$ Pascal. The temperature must be in Kelvin and must be between $283.15$ and $623.15$ Kelvin. Using these parameters, a preprocessing step constructs a two-dimensional table storing the brine density for the specified salinity as a function of pressure and temperature. For salinities above a threshold of $0.25$ mol/kg, the density is calculated directly using the Phillips correlation, $\rho_{Phillips}(P, T, S)$ , using the expression:

$\rho_{Phillips}(P, T, S) = A + Bx + Cx^2 + Dx^3$

$x = c_1 \exp(a_1 S) + c_2 \exp(a_2 T) + c_3 \exp(a_3 P)$

We refer the reader to Phillips et al. (1981), equations (4) and (5), pages 14 and 15 for the definition of the coefficients involved in the previous equation. For low salinities ( $S < 0.25$ ), the model linearly interpolates between the pure water density, $\rho_{pure}(P, T)$ , and the Phillips correlation to maintain physical continuity towards the pure water limit:

$\rho_{brine}(P, T) = \begin{cases} \rho_{Phillips}(P, T, S) & \text{if } S \ge 0.25 \\ \frac{S}{0.25} \rho_{Phillips}(P, T, S) + \left(1 - \frac{S}{0.25}\right) \rho_{pure}(P, T) & \text{if } S < 0.25 \end{cases}$

The 2D table of dimensional volume shifts for water, $c_{H2O}(P, T)$ , is generated by taking the difference between the target molar volume of the brine and the molar volume of pure water predicted by the chosen Equation of State, $v_{EOS,H2O}$ :

$c_{H2O}(P, T) = \frac{MW_{brine}}{\rho_{brine}(P, T)} - v_{EOS,H2O}(P, T)$

During the simulation, this pre-computed volume shift couples the empirical correlation to the overall phase EOS density calculation. The corrected molar volume of the aqueous phase, $v$ , is obtained by adding the water component’s volume shift weighted by its mole fraction, $x_{H2O}$ , to the uncorrected phase EOS molar volume, $v_{EOS}$ :

$v = v_{EOS} + x_{H2O} c_{H2O}(P, T)$

The phase molar density, $\rho_{molar}$ , is calculated as the inverse of the corrected molar volume:

$\rho_{molar} = \frac{1}{v}$

The phase mass density, $\rho_{mass}$ , is then computed by multiplying the phase molar density by the apparent molar weight of the phase mixture:

$\rho_{mass} = \rho_{molar} \left( x_{H2O} MW_{brine} + \sum_{i \neq H2O} x_i MW_i \right)$

Note

While the volume shift tightly couples the Phillips correlation to the EOS density, the brine viscosity is not derived from this density. Instead, the phase viscosity is modelled independently using a distinct temperature- and salinity-dependent multiplier applied to the pure water viscosity.

Parameters

salinity: Brine salinity. Unit: [mol/kg].
saltMolarWeight: The molar weight of the salt component. Unit: [kg/mol].

Immiscible water density

A simplified exponential density model is available for pure, immiscible water phases. It evaluates the mass density, $\rho$ , based on isothermal compressibility and thermal expansion from a reference state:

$\rho = \rho_{ref} \exp(c_w (P - P_{ref})) \exp(-\alpha_w (T - T_{ref}))$

The phase molar density, $\rho_{molar}$ , is then calculated by dividing the mass density by the molecular weight of water, $MW_{H2O}$ :

$\rho_{molar} = \frac{\rho}{MW_{H2O}}$

Parameters

waterReferencePressure: Reference pressure for the water density calculation. Unit: [Pa].
waterReferenceTemperature: Reference temperature for the water density calculation. Unit: [K].
waterDensity: Water mass density at the reference pressure and temperature. Unit: [kg/m^3].
waterCompressibility: Isothermal compressibility of water. Unit: [1/Pa].
waterExpansionCoefficient: Volumetric coefficient of thermal expansion of water. Unit: [1/K].

Viscosity

Viscosity dictates the fluid’s resistance to shear flow and fundamentally controls phase mobilities in porous media. The framework provides multiple approaches to evaluate viscosity.

Constant viscosity model

The simplest approach applies a static, user-defined viscosity. While the viscosity does not change with state variables such as pressure, temperature, or composition, it is possible to assign different constant viscosities to the different fluid phases. This is recommended solely for idealised simulations, synthetic testing, or when fluid flow is dominated entirely by pressure gradients rather than mobility contrasts.

Parameters

constantPhaseViscosity: The constant viscosity value applied to each phase. Unit: [Pa.s].

Lohrenz-Bray-Clark (LBC) viscosity model

The Lohrenz, Bray, and Clark correlation (Lohrenz et al., 1964) is the industry-standard methodology for calculating the viscosity of compositional hydrocarbon mixtures.

The model evaluates viscosity by combining a low-pressure dilute gas contribution with a dense fluid residual term.

Dilute gas viscosity

The low-pressure, dilute gas viscosity $\mu_{i}^*$ of each pure component $i$ is estimated as a function of temperature using the correlations of Stiel and Thodos (1961). This relies on the viscosity-reducing parameter $\xi_i$ , defined as:

$\xi_i = \frac{T_{c,i}^{1/6}}{M_i^{1/2} P_{c,i}^{2/3}}$

where $T_{c,i}$ is the critical temperature [K], $P_{c,i}$ is the critical pressure [atm], and $M_i$ is the molar weight [g/mol].

For standard non-polar components, the dilute viscosity (in centipoise) depends on the reduced temperature $T_{r,i} = T / T_{c,i}$ :

$\mu_{i}^* = \frac{34.0 \times 10^{-5} T_{r,i}^{0.94}}{\xi_i} \quad \text{for } T_{r,i} \le 1.5$

$\mu_{i}^* = \frac{17.78 \times 10^{-5} (4.58 T_{r,i} - 1.67)^{0.625}}{\xi_i} \quad \text{for } T_{r,i} > 1.5$

For hydrogen gas ( $M_i < 2.1 \times 10^{-3}$ kg/mol), the specific correlation is:

$\mu_{H_2}^* = 90.71 \times 10^{-5} (0.1375 T - 1.67)^{0.625}$

These pure viscosities are then blended to find the dilute mixture viscosity, $\mu^*$ . The user may select between the Herning-Zipperer, Wilke, or Brokaw mixing rules to aggregate the individual components with mole fractions $y_i$ .

The Herning-Zipperer mixing rule (Herning and Zipperer, 1936) evaluates the mixture viscosity as:

$\mu^* = \frac{\sum_i y_i \mu_i^* \sqrt{M_i}}{\sum_i y_i \sqrt{M_i}}$

The Wilke mixing rule (Wilke, 1950) evaluates the mixture viscosity as:

$\mu^* = \sum_i \frac{y_i \mu_i^*}{\sum_j y_j \phi_{ij}}$

where the interaction parameter $\phi_{ij}$ is defined as:

$\phi_{ij} = \frac{\left[1 + (\mu_i^*/\mu_j^*)^{1/2} (M_i/M_j)^{-1/4}\right]^2}{\sqrt{8(1 + M_i/M_j)}}$

The Brokaw mixing rule (Brokaw, 1968) evaluates the mixture viscosity as:

$\mu^* = \sum_i \frac{y_i \mu_i^*}{\sum_j y_j \phi_{ij} \sqrt{\mu_i^*/\mu_j^*}}$

where the interaction parameter $\phi_{ij}$ depends on the molar weight ratios $A_{ij} = M_i / M_j$ and $B_{ij} = [4 M_i M_j / (M_i + M_j)^2]^{0.25}$ :

$\phi_{ij} = \frac{B_{ij}}{\sqrt{A_{ij}}} \left[ 1 + \frac{A_{ij} - A_{ij}^{0.45}}{2 + 2A_{ij} + B_{ij}(1+A_{ij}^{0.45})/(1+B_{ij})} \right]$

Dense fluid contribution

A residual viscosity, representing the effects of elevated pressure and density, is calculated using the Lohrenz et al. (1964) correlation. This evaluates a fourth-order polynomial dependent on the reduced density of the phase.

First, the mixture critical parameters are calculated using Kay’s mixing rule:

$T_c = \sum_i y_i T_{c,i}, \quad P_c = \sum_i y_i P_{c,i}, \quad V_c = \sum_i y_i V_{c,i}, \quad M = \sum_i y_i M_i$

The mixture viscosity-reducing parameter $\xi$ and reduced density $\rho_r$ are defined using the phase molar density $\rho_m$ :

$\xi = \frac{T_c^{1/6}}{M^{1/2} P_c^{2/3}}, \quad \rho_r = \rho_m V_c$

The dense fluid polynomial $f(\rho_r)$ is then evaluated:

$f(\rho_r) = 0.1023 + 0.023364 \rho_r + 0.058533 \rho_r^2 - 0.040758 \rho_r^3 + 0.0093324 \rho_r^4$

This residual is added to the dilute gas viscosity to obtain the final phase viscosity $\mu$ :

$\mu = \mu^* + \frac{f(\rho_r)^4 - 10^{-4}}{\xi}$

This model is strongly recommended for all miscible gas, volatile oil, and general compositional hydrocarbon simulations.

The model relies on critical component volumes; when these are not provided, they are estimated internally using a correlation due to Ihmels (2010)

$V_c = \frac{2.215 \times 10^{-6} T_c}{0.025 + 10^{-6} P_c}$

Parameters

viscosityMixingRule: Defines the dilute gas mixing rule. Valid options are HerningZipperer, Wilke, or Brokaw.
componentCriticalVolume: The critical volumes of each component, required for the reduced density calculation. Unit: [m^3/mol].

Phillips brine viscosity model

For aqueous phases, the Phillips brine viscosity model uses empirical correlations from Phillips et al. (1981) to calculate viscosity as a function of temperature and salinity. In this specific implementation, the brine viscosity is evaluated independent of pressure. This model is recommended for saline aquifers and CO2 sequestration.

The model determines the final brine viscosity $\mu$ by applying a salinity- and temperature-dependent multiplier to the pure water viscosity $\mu_w$ :

$\mu = \mu_w(T) \left[ c_0(S) + c_1(S) T \right]$

where $T$ is the temperature in degrees Celsius and $S$ is the salinity in mol/kg. The pure water viscosity $\mu_w(T)$ is not evaluated from a closed-form analytical equation; rather, it is interpolated from a pre-computed tabulated function of pure water saturation viscosities corresponding to the given temperature.

The coefficients $c_0(S)$ and $c_1(S)$ dictate the scaling effect of the dissolved salts and are evaluated as follows:

$c_0(S) = 1 + a S + b S^2 + c S^3$

$c_1(S) = d (1 - \exp(k S))$

The correlation uses the specific empirical constants $a = 0.0816$ , $b = 0.0122$ , $c = 0.000128$ , $d = 0.000629$ , and $k = -0.7$ .

Parameters

salinity: The salinity of the brine. Unit: [mol/kg].
temperatureCoordinates: List of temperature values for the interpolation of tabulated pure water properties. Unit: [K].

Immiscible water viscosity

A simplified exponential viscosity model is available for pure, immiscible water phases. It evaluates the viscosity, $\mu$ , based on pressure and temperature variations from a reference state:

$\mu = \mu_{ref} \exp(c_{\mu} (P - P_{ref})) \exp(-\alpha_{\mu} (T - T_{ref}))$

where $\mu_{ref}$ is the reference viscosity, $c_{\mu}$ is the viscosity compressibility, and $\alpha_{\mu}$ is the viscosity expansion coefficient.

Parameters

waterReferencePressure: Reference pressure for the water viscosity calculation. Unit: [Pa].
waterReferenceTemperature: Reference temperature for the water viscosity calculation. Unit: [K].
waterViscosity: Water viscosity at the reference pressure and temperature. Unit: [Pa.s].
waterViscosityCompressibility: The compressibility (normalized derivative with respect to pressure) of the water viscosity. Unit: [1/Pa].
waterViscosityExpansionCoefficient: The coefficient of thermal expansion (normalized derivative with respect to temperature) of water viscosity. Unit: [1/K].

Enthalpy

For non-isothermal simulations, such as thermal recovery processes or geothermal operations, the accurate calculation of phase enthalpies and internal energies is required to resolve the energy balance equation.

Compositional enthalpy model

The compositional enthalpy model for a multiphase fluid system computes the total enthalpy of the phase as the sum of an ideal (polynomial) enthalpy and a real fluid correction (or departure enthalpy), which is derived from the equation of state.

The total enthalpy of a fluid phase is defined as:

$H_{total} = H_{ideal} + H_{dep}$

Ideal gas contribution

The ideal enthalpy represents the enthalpy of the mixture assuming it behaves as an ideal gas, relying on the heat capacity of the individual components. The specific heat capacity $C_{p,i}$ for each component $i$ is modeled as a 4th-order Poling polynomial (Poling et al., 2001) based on the temperature difference from a reference state, $\Delta T = T - T_{ref}$ :

$C_{p,i}(T) = a_{i,0} + a_{i,1}\Delta T + a_{i,2}\Delta T^2 + a_{i,3}\Delta T^3 + a_{i,4}\Delta T^4$

To find the ideal enthalpy for a component, this heat capacity is integrated with respect to temperature from the reference temperature to the current temperature, and then added to a baseline reference enthalpy $H_{ref,i}$ :

$H_{ideal,i}(T) = H_{ref,i} + \Delta T \left( a_{i,0} + \frac{a_{i,1}}{2}\Delta T + \frac{a_{i,2}}{3}\Delta T^2 + \frac{a_{i,3}}{4}\Delta T^3 + \frac{a_{i,4}}{5}\Delta T^4 \right)$

The total ideal enthalpy for the phase is simply the mole-fraction weighted sum of the individual component ideal enthalpies:

$H_{ideal} = \sum_i x_i H_{ideal,i}(T)$

Departure function

The dimensionless enthalpy departure is given by:

$\frac{H_{dep}}{RT} = Z - 1 + \frac{A + T\frac{\partial A}{\partial T}}{B(\delta_1 - \delta_2)} \ln \left( \frac{Z + \delta_1 B}{Z + \delta_2 B} \right)$

where $Z$ is the compressibility factor of the mixture.

This model is recommended whenever non-isothermal compositional physics are required.

Mass formulation conversion

When the simulator is configured to use the mass formulation, the underlying inputs and thermodynamic calculations remain entirely in the molar formulation. The final phase molar enthalpy $H_{total}$ [J/mol] is converted to a mass-based phase enthalpy $H_{mass}$ [J/kg] by dividing by the phase molar weight $M_{phase}$ :

$H_{mass} = \frac{H_{total}}{M_{phase}}$

where the phase molar weight is the mole-fraction weighted sum of the component molar weights $M_i$ :

$M_{phase} = \sum_i x_i M_i$

Because all inputs must be provided in molar units, if a user’s specific heat capacity coefficients or reference enthalpies are originally defined in a mass formulation (e.g., $\hat{a}_{i,k}$ in [J/(kg K^n)]), they must be scaled by the respective component’s molar weight $M_i$ [kg/mol] prior to input to yield the molar coefficients $a_{i,k}$ :

$a_{i,k} = \hat{a}_{i,k} M_i$

Implementing the Michaelides Model

The the Michaelides (1981) enthalpy model which is used for the CO2-brine model (CO2-brine model) can be implemented using the polynomial framework by setting specific values for the water component.

The rational function $Y(m)$ below serves as a universal algebraic template to compute the required simulator inputs as a function of the salt molality $m$ [mol/kg]. Because the underlying thermodynamic derivations rely on fixed mixing rules and molecular weights, the calculation for every target parameter condenses into this single functional form:

$Y(m) = \frac{p_0 + p_1 m + p_2 m^2 + p_3 m^3}{m + 55.5247}$

To determine the input value for a specific simulator parameter, locate that parameter in the table below and substitute its corresponding row of coefficients ( $p_0$ through $p_3$ ) along with your desired molality $m$ into the $Y(m)$ equation. The evaluated result provides the direct numerical input required for the value.

The values in the table are hard-coded to a reference temperature of 293.15K and molar weights of 18.01 g/mol for water and 58.44 g/mol for salt.

Parameter	$p_0$	$p_1$	$p_2$	$p_3$
$H_{ref,w}$	6.49471e+04	-2.97460e+04	-1.19193e+04	8.46747e+02
$a_{w,0}$	4.64595e+03	4.69199e+02	2.29413e+02	-1.56721e+01
$a_{w,1}$	-7.53304e+00	-5.72452e+00	-2.69521e+00	1.84321e-01
$a_{w,2}$	3.73590e-02	1.20531e-02	9.02790e-03	-6.15000e-04
$a_{w,3}$	0.00000e+00	0.00000e+00	0.00000e+00	0.00000e+00
$a_{w,4}$	0.00000e+00	0.00000e+00	0.00000e+00	0.00000e+00

Note

Because the original Michaelides enthalpy model relies on a 3rd-order polynomial fit with respect to temperature (Michaelides, 1981), the higher-order simulator heat capacity coefficients ( $a_{w,3}$ and $a_{w,4}$ ) evaluate exactly to zero.

This formulation is bounded by specific thermodynamic assumptions: it achieves about ±3% accuracy relative to standard geothermal fluid tables, remains valid for saturation temperatures up to 320 °C (with possible divergence below 0 °C), and applies to salinity levels up to roughly 5-6 mol/kg using equivalent NaCl content. It further assumes fixed molar masses of 18.01 g/mol for water and 58.44 g/mol for NaCl, meaning any deviation without polynomial refitting will introduce scaling errors. Finally, the model is limited to liquid-phase (brine) enthalpy, and phase transitions such as flashing to vapor must be handled separately using appropriate gas mixing rules.

Parameters

enthalpyReferenceTemperature: (Optional) A scalar defining the reference temperature $T_{ref}$ with a default of 298.15 K (25 C). Unit: [K].
referenceEnthalpy: (Optional) A 2D array mapping the reference enthalpy for each component within each phase at the reference temperature. The dimensions must match [number of phases, number of components]. Defaults to zero. Unit: [J/mol].
componentHeatCapacityCoefficients: (Required) A 2D array containing the 5 polynomial coefficients ( $a_0$ to $a_4$ ) for the specific heat capacity of each component. The dimensions must strictly be [number of components, 5]. To maintain dimensional consistency with a heat capacity in [J/(mol K)] and temperature difference in [K], the units for the specific coefficients are:
- $a_{i,0}$ : [J/(mol K)]
- $a_{i,1}$ : [J/(mol K^2)]
- $a_{i,2}$ : [J/(mol K^3)]
- $a_{i,3}$ : [J/(mol K^4)]
- $a_{i,4}$ : [J/(mol K^5)]

Miscellaneous

Unit and Variable Conversions

Thermodynamic models, including Equations of State and Flash calculations, are fundamentally derived from molar balances and molar properties. Consequently, the internal computations within the compositional fluid framework are strictly executed using molar variables (e.g., molar fractions, molar densities).

Mass Formulation Support

While the core thermodynamics rely on molar formulations, many reservoir simulators frame their governing conservation equations in terms of mass to ensure strict mass conservation across numerical schemes.

To bridge this gap seamlessly, the fluid model supports an automated conversion layer:

Entry Conversion: Upon entering the fluid property update routine, if the overall fluid state is provided in mass fractions, the framework dynamically converts these inputs into molar fractions using the designated component molecular weights.
Internal Computation: All stability, flash, density, and viscosity calculations proceed purely in molar space.
Exit Transformation: Before the update routine concludes, the resulting phase compositions, phase fractions, and molar densities are transformed back into mass variables. Additionally, rigorous chain-rule transformations are applied to all derivatives (with respect to pressure, temperature, and composition) to ensure the Jacobian matrices populated by the simulator remain mathematically consistent with the mass-based primary variables.

Note on Component Limits

When compositional fluid models are evaluated in standalone mode (for example, using the PVTDriver task (PVT Driver)), the maximum number of allowed components is 9. However, when these fluid models are actively coupled with a flow solver during a full simulation, the maximum number of components is restricted to 5.

List of available compositional fluid models

The following table outlines the available compositional fluid models (specified via their XML elements), detailing the number of phases they support, their underlying thermodynamic phase split (flash) methods, the phase density and viscosity models used, and their recommended applications:

Compositional fluid models

CompositionalTwoPhaseFluid XML Element: CompositionalTwoPhaseFluid

Number of phases: 2.
Possible applications: Idealised simulations, synthetic testing, or when fluid flow is dominated entirely by pressure gradients.
Phase split method: Negative Two-Phase Flash.
Phase density models: Compositional density.
Phase viscosity models: Constant viscosity.

CompositionalTwoPhaseFluidLohrenzBrayClark XML Element: CompositionalTwoPhaseFluidLohrenzBrayClark

Number of phases: 2.
Possible applications: Miscible gas injection, primary depletion of volatile oils, and scenarios where fluid compositions change drastically.
Phase split method: Negative Two-Phase Flash.
Phase density models: Compositional density.
Phase viscosity models: Lohrenz-Bray-Clark viscosity.

CompositionalTwoPhaseFluidPhillipsBrine XML Element: CompositionalTwoPhaseFluidPhillipsBrine

Number of phases: 2.
Possible applications: Carbon sequestration (CO2-brine systems) or water-flooding in oil reservoirs.
Phase split method: Negative Two-Phase Flash.
Phase density models: Phillips brine density model for the aqueous phase; Compositional density for the gas phase.
Phase viscosity models: Phillips brine viscosity model for the aqueous phase; Lohrenz-Bray-Clark viscosity for the gas phase.

CompositionalThreePhaseFluidLohrenzBrayClark XML Element: CompositionalThreePhaseFluidLohrenzBrayClark

Number of phases: 3.
Possible applications: Multi-phase reservoir simulations where water solubility in the hydrocarbon phases is negligible.
Phase split method: Immiscible Water Flash.
Phase density models: Immiscible water density for the aqueous phase; Compositional density for the hydrocarbon phases.
Phase viscosity models: Immiscible water viscosity for the aqueous phase; Lohrenz-Bray-Clark viscosity for the hydrocarbon phases.

CompositionalTwoPhaseKValueFluidLohrenzBrayClark XML Element: CompositionalTwoPhaseKValueFluidLohrenzBrayClark

Number of phases: 2.
Possible applications: Thermal recovery processes or relatively stable black-oil systems.
Phase split method: K-Value Flash.
Phase density models: Compositional density.
Phase viscosity models: Lohrenz-Bray-Clark viscosity.

CompositionalTwoPhaseKValueFluidPhillipsBrine XML Element: CompositionalTwoPhaseKValueFluidPhillipsBrine

Number of phases: 2.
Possible applications: Saline aquifers and CO2 sequestration with stable compositions where pre-tabulated K-values are sufficient.
Phase split method: K-Value Flash.
Phase density models: Phillips brine density model; Compositional density.
Phase viscosity models: Phillips brine viscosity model; Lohrenz-Bray-Clark viscosity.

References

Brokaw, R. S. (1968). Viscosity of gas mixtures. Washington, D.C.: National Aeronautics and Space Administration.
Chabab, S., Kerkache, H., Bouchkira, I., Poulain, M., Baudouin, O., Moine, É., Ducousso, M., Hoang, H., Galliéro, G., and Cézac, P. (2024). Solubility of H2 in water and NaCl brine under subsurface storage conditions: Measurements and thermodynamic modeling. International Journal of Hydrogen Energy, 50, 648-658. https://doi.org/10.1016/j.ijhydene.2023.10.290
Herning, F., and Zipperer, L. (1936). Beitrag zur Berechnung der Zähigkeit Technischer Gasgemische aus den Zähigkeitswerten der Einzelbestandteile; Calculation of the Viscosity of Technical Gas Mixtures from Viscosity of the individual Gases. Das Gas- und Wasserfach, 79, 49-54 and 69-73.
Ihmels, E. C. (2010). The Critical Surface. Journal of Chemical & Engineering Data, 55 (9), 3474-3480. https://doi.org/10.1021/je100167w
Lohrenz, J., Bray, B. G., and Clark, C. R. (1964). Calculating Viscosities of Reservoir Fluids From Their Compositions. Journal of Petroleum Technology, 16 (10), 1171-1176. https://doi.org/10.2118/915-pa
Michelsen, M. L. (1982a). The isothermal flash problem. Part I. Stability. Fluid Phase Equilibria, 9 (1), 1-19. https://doi.org/10.1016/0378-3812(82)85001-2
Michelsen, M. L. (1982b). The isothermal flash problem. Part II. Phase-split calculation. Fluid Phase Equilibria, 9 (1), 21-40. https://doi.org/10.1016/0378-3812(82)85002-4
Péneloux, A., Rauzy, E., and Fréze, R. (1982). A consistent correction for Redlich-Kwong-Soave volumes. Fluid Phase Equilibria, 8 (1), 7-23. https://doi.org/10.1016/0378-3812(82)80002-2
Peng, D.-Y., and Robinson, D. B. (1976). A New Two-Constant Equation of State. Industrial & Engineering Chemistry Fundamentals, 15 (1), 59-64. https://doi.org/10.1021/i160057a011
Phillips, S. L., Igbene, A., Fair, J. A., Ozbek, H., and Tavana, M. (1981). A technical databook for geothermal energy utilization. Lawrence Berkeley National Laboratory.
Poling, B. E., Prausnitz, J. M., and O’Connell, J. P. (2001). The Properties of Gases and Liquids. McGraw-Hill.
Rachford, H. H., Jr., and Rice, J. D. (1952). Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium. Journal of Petroleum Technology, 4 (10), 19-3. https://doi.org/10.2118/952327-g
Robinson, D. B., and Peng, D. (1978). The Characterization of the Heptanes and Heavier Fractions for the GPA Peng-Robinson Programs. United States: Gas Processors Association.
Soave, G. (1972). Equilibrium constants from a modified Redlich-Kwong equation of state. Chemical Engineering Science, 27 (6), 1197-1203. https://doi.org/10.1016/0009-2509(72)80096-4
Søreide, I., and Whitson, C. H. (1992). Peng-Robinson predictions for hydrocarbons, CO2, N2, and H2 S with pure water and NaCI brine. Fluid Phase Equilibria, 77, 217-240. https://doi.org/10.1016/0378-3812(92)85105-h
Stiel, L. I., and Thodos, G. (1961). The viscosity of nonpolar gases at normal pressures. AIChE Journal, 7 (4), 611-615. https://doi.org/10.1002/aic.690070416
Whitson, C. H., and Torp, S. B. (1981). Evaluating Constant Volume Depletion Data. SPE Annual Technical Conference and Exhibition. https://doi.org/10.2118/10067-ms
Wilke, C. R. (1950). A Viscosity Equation for Gas Mixtures. The Journal of Chemical Physics, 18 (4), 517-519. https://doi.org/10.1063/1.1747673
Wilson, G. M. (1969). A modified Redlich-Kwong equation of state, application to general physical data calculations. 65th National AIChE Meeting, Cleveland, OH, 15.