Hydro-Thermal Process (HT)

Introduction

The Hydro-Thermal (HT) process models coupled groundwater flow and heat transport in porous media. Both subprocesses are governed by coupled parabolic partial differential equations, linked via temperature dependent fluid properties and advective heat transport.

Key features

• Monolithic and staggered coupling schemes.
• Fluid compressibility and optional solid thermal expansion coupling via the Biot coefficient.
• Hydrodynamic thermal dispersion (longitudinal and transversal dispersivities).
• Numerical stabilisation for advection-dominated problems.
• Fracture flow support via aperture size parameter.
• Surface flux calculation.

Physical variables

• Primary variables: pressure $p$ and temperature $T$.
• Secondary variable: Darcy velocity $\mathbf{q}$.

Theoretical background

Both the flow and heat transport processes are derived from integral conservation laws. At the moment there is no coupling by source or sink terms, i.e., the coupling is implemented only through density changes due to temperature changes in the buoyancy term of the groundwater flow. The coupling scheme is referred to as the Boussinesq approximation.

Mass balance equation

The Darcy velocity is given by:

where $\boldsymbol{\kappa}$ is the permeability, $p$ is the pore pressure, $\varrho_f$ is the liquid density, and $\mathbf{g}$ is the gravitational force.

The mass balance equation reads:

where $Q_H$ is the sink or source term.

Assuming the solid matrix is compressible (i.e., $\frac{\partial \phi }{\partial t}\neq0$), the mass balance equation can be expanded as

where $S_s$ is the specific storage of the solid phase, $\alpha_B$ is the Biot coefficient, $\alpha_T^s$ is the linear solid thermal expansivity, and $T$ is the temperature.

The value of $S_s$ can be computed as $(\alpha_B - \phi)(1-\alpha_B) / K$, where $K$ is the drained bulk modulus. Equivalently, $S_s = (\alpha_B - \phi) / K_s$ with $K_s = K / (1 - \alpha_B)$ the intrinsic bulk modulus of the solid phase.

In certain scenarios, such as far-field simulations, the fluid density is often assumed constant. Consequently, the volume balance equation can be derived by dividing the mass balance equation by the constant fluid density:

Note: In OGS, the fluid part of the thermal expansion term is always present; its solid part $3(\alpha_B-\phi)\alpha_T^s$ is optional.

Heat transport equation

where:

Weak formulation

The weak form is obtained by multiplying the strong form by a test function $v \in H_0^1(\Omega)$ and integrating over the domain.

For the pressure field, the weak form is:

for the mass balance, and

for the volume balance, where $\mathbf{n}$ is the outward unit normal to the domain boundary $\Gamma$.

Note: The term $\frac{\phi}{\varrho_f} \frac{\partial \varrho_f}{\partial p}$ represents the fluid compressibility contribution to the storage coefficient.

For the temperature field, the weak form is:

Finite element discretization

Both equations are discretized into the standard form $\mathbf{M} \dot{\mathbf{u}} + \mathbf{K} \mathbf{u} = \mathbf{f}$ by substituting the Galerkin discretization $p \approx \sum_j N_j a_j$ and choosing $v = N_i$.

Pressure equation

For the mass balance equation, the resulting discretized matrices and vectors are:

In OGS, the thermal expansion term is always assembled; only its solid part $3(\alpha_B - \phi)\alpha_T^s$ is optional and drops out when the solid thermal expansivity is undefined. It is added to the right hand vector for the staggered scheme as

and it is added to the mass matrix for the monolithic scheme as

The corresponding terms for the volume balance equation are obtained by dividing all the above integrals by the fluid density.

Temperature equation

The advection term is part of $\mathbf{K}^T$ in both schemes; they differ only in where the Darcy velocity $\mathbf{q}$ comes from, namely from the pressure solution of the same equation system in the monolithic scheme and from the previous solution of the pressure equation in the staggered scheme.

The right hand side vector is given by

Definition in the project file

The HT process is declared in the <processes> block of the project file.

Monolithic scheme (default)

<process>
    <name>HydroThermal</name>
    <type>HT</type>
    <integration_order>2</integration_order>
    <process_variables>
        <temperature>T</temperature>
        <pressure>p</pressure>
    </process_variables>
    <specific_body_force>0 -9.81</specific_body_force>
    <secondary_variables>
        <secondary_variable name="darcy_velocity"/>
    </secondary_variables>
</process>

Balance equation type: mass or volume

Similar to other hydraulics-related processes, HT supports an optional <equation_balance_type> tag with values mass or volume to select the balance equation type for the hydraulics process. For example:

<process>
    <name>HydroThermal</name>
    <type>HT</type>
    <integration_order>2</integration_order>
    <equation_balance_type>mass</equation_balance_type>
    ...
</process>

By default, the value is volume. The volume balance requires the liquid phase density property to be of type Constant; with any other property type, for instance Linear, Function or WaterDensityIAPWSIF97Region1, OGS issues a fatal error prompting the user to add the tag. In such cases, the input values for Neumann boundary conditions and source/sink terms must be adjusted accordingly, for example:

• Neumann condition: from volume rate per area (SI unit: $\text{m}\cdot\text{s}^{-1}$) to mass rate per area (SI unit: $\text{kg}\cdot\text{m}^{-2}\text{s}^{-1}$).
• source/sink: from volume rate (SI unit: $\text{m}^{3}\text{s}^{-1}$) to mass rate (SI unit: $\text{kg}\cdot\text{s}^{-1}$).

The choice of the balance equation type is orthogonal to the material property combinations listed in Material property restrictions below: the mass balance scales $\mathbf{M}^p$, $\mathbf{K}^p$, $\mathbf{f}^p$ and the thermal expansion term by $\varrho_f$ uniformly, and it does not change which combination applies.

Thermal expansion in the mass/volume balance equation

The fluid part of the thermal expansion term, $-\frac{\phi}{\varrho_f}\frac{\partial \varrho_f}{\partial T}$, is always computed. The solid part, $3(\alpha_B - \phi)\alpha_T^s$, is added only if the linear solid thermal expansivity (property name: thermal_expansivity) and the Biot’s coefficient (property name: biot_coefficient) are defined in the input project file. If linear solid thermal expansivity is anisotropic, the average of its components can be used.

Staggered scheme

To use the staggered coupling scheme, add:

<coupling_scheme>staggered</coupling_scheme>

In the staggered scheme, the heat transport equation (process ID 0) and the hydraulic equation (process ID 1) are solved sequentially.

Process variables

The HT process requires two process variables: temperature and pressure. Each should have 1 component. For more details, see Process variables .

Media properties

The HT process requires properties for the porous medium, the liquid phase, and the solid phase.

Medium properties

Material property restrictions for Biot’s coefficient, porosity, storage, and thermal expansivity combinations

The Biot coefficient and the solid thermal expansivity are optional material parameters used to account for the solid phase contribution to thermal expansion. Once these parameters are set in the project file, they are associated with the specific storage of the solid phase $S_s$, which is input via storage, and the effective thermal expansivity $\beta_{\mathrm{eff}}=\left(3(\alpha_B-\phi)\alpha_T^s-\frac{\phi}{\varrho_f}\frac{\partial \varrho_f}{\partial T}\right)$. This can be treated as a parameter combination when both the Biot coefficient and the solid thermal expansivity are defined. The rules are:

1. $\alpha_B$ is read only when $\alpha_T^s$ is defined. Therefore:
   • $\alpha_T^s$ defined, $\alpha_B$ undefined → fatal error (undefined property).
   • $\alpha_B$ defined, $\alpha_T^s$ undefined → $\alpha_B$ is silently ignored (never read).
2. The specific storage $S_s$ given by the solid-phase storage property is always read.
3. The combination $\alpha_B = 1$ and $S_s \neq 0$ leads to a fatal error because $S_s=(\alpha_B-\phi)(1-\alpha_B)/K$. The guard is only reached when $\alpha_T^s$ is defined, since that is when $\alpha_B$ is read at all. The evaluated values are compared once per element, at its integration points, during the initialisation of the local assemblers, at $t=0$ and $\Delta t=0$. That is the literal time zero, not the initial time of the time loop, which may be a different value. This covers every property that depends on neither the primary variables nor the time, whatever its type. A property depending on the time is compared at $t=0$ only, which the simulation need not pass through at all, so a violation occurring at any later time is not caught here. A property that does depend on the primary variables evaluates to NaN at that point and is not rejected: a Biot coefficient of NaN makes the comparison $\alpha_B = 1$ false, and a specific storage of NaN is excluded explicitly. Both cases are decided by the same comparison, which is repeated at every integration point during the assembly, where the properties have values. With $\alpha_T^s$ undefined the combination is accepted silently, because $\alpha_B$ is then never read.

The ranges of these four parameters are listed in the following table:

This means that if these parameters are left unchecked and fall outside the stated ranges, only the combination $\alpha_B=1,\,S_s\neq0$ is detected, and only when $\alpha_T^s$ is defined.

The relation $S_s=(\alpha_B-\phi)(1-\alpha_B)/K$ also constrains the reverse direction: $S_s=0$ with $\alpha_B<1$ implies $\alpha_B=\phi$ or $K=\infty$. That combination is not guarded, so row 4 below runs without a warning even though only $\alpha_B=1$ makes $S_s=0$ consistent with that relation. Row 1 carries no such inconsistency, since $\alpha_B$ is undefined there and never read; its risk is the vanishing storage term alone.

The following table summarises the behaviour of the different combinations of the four parameters (Legend: ✓ defined, ✗ undefined):

Liquid phase properties

$\large^{\star}$ Functional dependencies (e.g. on temperature) can be specified. See the OGS User Guide .

Solid phase properties

Features

Specific body force

The gravity vector is specified as:

<specific_body_force>0 -9.81</specific_body_force>

Aperture size

For lower-dimensional fracture elements, an aperture size parameter can be specified:

<aperture_size>
    <parameter>fracture_aperture</parameter>
</aperture_size>

Surface flux calculation

Surface flux output can be configured:

<calculatesurfaceflux>
    ...
</calculatesurfaceflux>

Numerical stabilisation

The HT process supports numerical stabilisation methods for advection-dominated transport.

Benchmarks

See the OGS benchmark gallery for Hydro-Thermal examples .