In Bourgeat et al. (2013) a numerical experiment is presented which investigates the transport behaviour of hydrogen in clay rock. This test example was investigated by several institutions using different simulation tools and the results were then compared in detail.
To verify the phase transition, the diffusive mass transport and the twophase behaviour, this numerical experiment was recalculated with the TH2M process class. The material parameters and boundary conditions used essentially correspond to those of the experiment and are summarised as follows:
Parameter  Symbol  Value  Unit 

binary diffusion coefficient  $D$  3.0e9  Pa 
viscosity liquid  $\mu_\text{LR}$  1.0e3  Pa s 
viscosity gas  $\mu_\text{GR}$  9.0e6  Pa s 
Henrycoefficient  $H$  7.65e6  mol Pa$^{1}$m$^{3}$ 
molar mass hydrogen  $M_\mathrm{H_2}$  2.0e3  kg mol$^{1}$ 
molar mass water  $M_\mathrm{H_2O}$  1.0e2  kg mol$^{1}$ 
density liquid  $\rho_\text{LR}$  eq. (1)  kg m$^{3}$ 
intrinsic permeability  $\mathbf{k}$  5.0e20  m$^2$ 
porosity  $\phi$  0.15  1 
The van Genuchten model was used as the saturation relation with $m=$0.329 and $p_b$=2.0e6 Pa, $s_\mathrm{L}^\mathrm{res}=$0.4 and $s_\mathrm{G}^\mathrm{res}=$0.0. The relative permeabilities $k^\mathrm{rel}\mathrm{L}$ and $k^\mathrm{rel}\mathrm{G}$ were determined using the van Genuchten and van GenuchtenMualem models, respectively.
Unlike in the cited article, the density of the liquid phase was not considered constant, but was represented by a linear equation of state of the form
$$ \rho_\text{LR}=\rho^\text{ref}\text{LR}\left(1+\chi{c,\text{L}}c^\text{C}_\text{L}\right) ;;;(1) $$
Here, $\rho^\text{ref}\text{LR}$=1000 kg$^{3}$ is a reference density and $\chi{c,\text{L}}=M_{H_2}{\rho^\text{ref}_\text{LR}}^{1}$=2.0e6 m$^{3}$mol$^{1}$ is the slope of the function with respect to the hydrogen concentration.
This deviation from the specifications in Bourgeat et al. (2013) is necessary because the phase transition model in TH2M evaluates the ratio of the densities of solvent and solution to calculate the dissolved gas fraction. Therefore, the equation of state must take into account the concentration of dissolved gases, even if the actual influence on the density seems negligible.
The experiment considers a 200 m long, porous quasi1D region, at the left edge of which a hydrogen mass flow is injected. This injection is kept constant for 500,000 years and then stopped. The chosen mass flow is low enough that all the injected hydrogen is initially dissolved in the liquid phase and is transported there primarily by diffusion.
The validity of the test is carried out by comparing the variables gas saturation, gas pressure and water pressure with the results of various codes from Bourgeat et al. (2013).
Bourgeat, A. P., Granet, S., & Smaï, F. (2013). Compositional twophase flow in saturated–unsaturated porous media: benchmarks for phase appearance/disappearance. Simulation of Flow in Porous Media, 1–29. https://doi.org/10.1515/9783110282245.81
import os
out_dir = os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out")
if not os.path.exists(out_dir):
os.makedirs(out_dir)
import numpy as np
from ogs6py.ogs import OGS
model = OGS(INPUT_FILE="bourgeat.prj", PROJECT_FILE=f"{out_dir}/modified.prj")
# This Jupyter notebook version of this test runs not as far as its cTestcounterpart,
# it'll stop after approx. 800 ka
model.replace_text(2.5e13, xpath="./time_loop/processes/process/time_stepping/t_end")
# The cTest version shows only a few outputtimesteps while this version is supposed to
# output much higher resolution time steps to be able to compare the results. Thus, every
# timestep will be written and the maximum timestepsize of the adaptive time stepping
# method will be reduced drastically
time_end = 1e11
model.replace_text(
time_end, xpath="./time_loop/processes/process/time_stepping/maximum_dt"
)
# The following for loop generates a text with output times,which is then replaced by
# the ogs6py API in the project file.
new_line = "\u000A"
timesteps = str(0.4 * 1e6 * 86400 * 365.25) + new_line
for t in np.arange(0.6, 1.1, 0.1):
timesteps += str(t * 1e6 * 86400 * 365.25) + new_line
model.replace_text(1, xpath="./time_loop/output/timesteps/pair/each_steps")
model.replace_text(timesteps, xpath="./time_loop/output/fixed_output_times")
model.write_input()
# Run OGS
model.run_model(logfile=f"{out_dir}/out.txt", args=f"o {out_dir} m .")
OGS finished with project file /var/lib/gitlabrunner/builds/e3EQ9HiK/1/ogs/build/releaseall/Tests/Data/TH2M/H2/dissolution_diffusion/phase_appearance/modified.prj.
Execution took 45.34620499610901 s
Project file written to output.
# Colors
cls = ["#e6191d", "#337fb8", "#4eae4c", "#984ea3", "#984ea3", "#feff32"]
import vtuIO
# Read PVDoutput
pvdfile = vtuIO.PVDIO(f"{out_dir}/result_bourgeat.pvd", dim=2)
point = {"A": (0.0, 0.0, 0.0)}
time = pvdfile.timesteps
saturation = pvdfile.read_time_series("saturation", point)
gas_pressure = pvdfile.read_time_series("gas_pressure", point)
liquid_pressure = pvdfile.read_time_series("liquid_pressure_interpolated", point)
num_results = [1.0  saturation["A"], gas_pressure["A"], liquid_pressure["A"]]
time_years = time / 365.2425 / 86400
WARNING: Default interpolation backend changed to VTK. This might result in
slight changes of interpolated values if defaults are/were used.
import pandas as pd
# Read the reference data from CSV files
refs = [
pd.read_csv("references/bourgeat_sG.csv"),
pd.read_csv("references/bourgeat_pGR.csv"),
pd.read_csv("references/bourgeat_pLR.csv"),
]
header = list(refs[0].keys())
indices = {"Gas saturation": 0, "Gas pressure": 1, "Liquid pressure": 2}
labels = ["$s_{G}$", "$p_{GR}$", "$p_{LR}$"]
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (12, 4)
# Loop over gas_saturation, gas_pressure, and liquid_pressure
for i in indices:
index = indices[i]
ref_time = refs[index]["time"]
fig1, (ax1, ax2) = plt.subplots(1, 2, sharey=True)
fig1.suptitle(i + r" vs. time at $x=0$m")
ax1.set_xscale("log")
ax1.set_xlabel("time / a", fontsize=12)
ax1.set_ylabel(labels[index], fontsize=12)
for r in range(1, len(refs[index].columns)):
ax1.plot(
ref_time,
refs[index][refs[index].keys()[r]],
linewidth=1,
linestyle="",
label=refs[index].keys()[r],
)
ax1.plot(
time_years,
num_results[index],
color="black",
linewidth=2,
linestyle="",
label="OGSTH$^2$M",
)
ax2.set_xlabel("time / a", fontsize=12)
for r in range(1, len(refs[index].columns)):
ax2.plot(
ref_time,
refs[index][refs[index].keys()[r]],
linewidth=1,
linestyle="",
label=refs[index].keys()[r],
)
ax2.plot(
time_years,
num_results[index],
color="black",
linewidth=2,
linestyle="",
label="OGSTH$^2$M",
)
ax1.legend()
fig1.savefig("results_sG_pGR_pLR.pdf")
After about 10,000 years, the dissolution capacity of the liquid phase is exhausted and a separate gas phase is formed. From this point on, an increase in water pressure can be observed. This is shown in the liquid pressure plot by comparison with the results shown in Bourgeat et al. (2013).
Regarding the accuracy (especially during times t<10,000 years), it should be noted that the reference graphs from the original paper were digitised manually; in the paper, the graphs were only linearly plotted (like those shown here in the right column), so no real statement can be made about the accuracy of the results for times t<10,000 years.
The comparison of the time evolution of the water pressure, gas pressure und gas saturation shows that the numerical solution calculated by OGS6TH2M is very close to the results shown in Bourgeat et al. (2013). Both the time of formation of the gas phase (approximately at t=15,000 a and the shape and magnitude of the resulting pressure rise and pressure drop due to switching off the source term are within the differences of the reference solutions.
It can be concluded that both the transport and the phase transition behaviour of the TH2M process have been successfully derived and implemented. The selfcontrolled appearance and disappearance of a free gas phase should be emphasised.
This article was written by Norbert Grunwald. If you are missing something or you find an error please let us know.
Generated with Hugo 0.117.0
in CI job 389152

Last revision: October 19, 2022