poet/bench/dolo

Description of dolo benchmark

Quick start

mpirun -np 4 ./poet dolo_diffu_inner.R dolo_diffu_inner_res
mpirun -np 4 ./poet --dht --interp dolo_interp_long.R dolo_interp_long_res

List of Files

• dolo_interp.R: POET input script for a 400x200 simulation grid
• dolo_diffu_inner_large.R: POET input script for a 400x200 simulation grid
• phreeqc_kin.dat: PHREEQC database containing the kinetic expressions for dolomite and celestite, stripped down from phreeqc.dat
• dol.pqi: PHREEQC input script for the chemical system

Chemical system

This model describes a simplified version of dolomitization of calcite, itself a complex and not yet fully understood natural process which is observed naturally at higher temperatures (Möller and De Lucia, 2020). Variants of such model have been widely used in many instances especially for the purpose of benchmarking numerical codes (De Lucia et al., 2021 and references therein).

We consider an isotherm porous system at 25 °C in which pure water is initially at equilibrium with calcite (calcium carbonate; brute formula: CaCO₃). An MgCl₂ solution enters the system causing calcite dissolution. The released excess concentration of dissolved calcium Ca²⁺ and carbonate (CO₃^2-) induces after a while supersaturation and hence precipitation of dolomite (calcium-magnesium carbonate; brute formula: CaMg(CO₃)₂). The overall dolomitization reaction can be written:

Mg²⁺ + 2 ⋅ calcite → dolomite + 2 ⋅ Ca²⁺

The precipitation of dolomite continues until enough Ca²⁺ is present in solution. Further injection of MgCl₂ changes its saturation state causing its dissolution too. After enough time, the whole system has depleted all minerals and the injected MgCl₂ solution fills up the domain.

Both calcite dissolution and dolomite precipitation/dissolution follow a kinetics rate law based on transition state theory (Palandri and Karhaka, 2004; De Lucia et al., 2021).

rate = -S_m k_m (1-SR_m)

where the reaction rate has units mol/s, S_m (m²) is the reactive surface area, k_m (mol/m²/s) is the kinetic coefficient, and SR is the saturation ratio, i.e., the ratio of the ion activity product of the reacting species and the solubility constant, calculated internally by PHREEQC from the speciated solution.

For dolomite, the kinetic coefficient results from the sum of two mechanisms, r_/acid/ and r_/neutral/:

rate_dolomite = S_dolomite (k_/acid/ + k_/neutral/) * (1 - SR_dolomite)

where:

k_/acid/ = 10^-3.19 e^{-36100 / R} ⋅ act(H⁺)^0.5

k_/neutral/ = 10^-7.53 e^{-52200 / R}

R (8.314462 J K^-1 mol^-1) is the gas constant.

Similarly, the kinetic law for calcite reads:

k_/acid/ = 10^-0.3 e^{-14400 / R} ⋅ act(H⁺)^0.5

k_/neutral/ = 10^-5.81 e^{-23500 / R}

The kinetic laws as implemented in the phreeqc_kin.dat file accepts one parameter which represents reactive surface area in m². For the benchmarks the surface areas are set to

• S_dolomite: 0.005 m²
• S_calcite: 0.05 m²

The initial content of calcite in the domain is of 0.0002 mol per kg of water. A constant partial pressure of 10^-1675 atm of O_2(g) is maintained at any time in the domain in order to fix the redox potential of the solution to an oxidizing state (pe around 9).

Note that Cl is unreactive in this system and only effects the computation of the activities in solution.

POET simulations

Several benchmarks based on the same chemical system are defined here with different grid sizes, resolution and boundary conditions. The transported elemental concentrations are 7: C(4), Ca, Cl, Mg and the implicit total H, total O and Charge as required by PHREEQC_RM.

dolo_diffu_inner.R

• Grid discretization: square domain of 1 ⋅ 1 m² discretized in 100x100 cells
• Boundary conditions: All sides of the domain are closed. Fixed concentration of 0.001 molal of MgCl₂ is defined in the domain cell (20, 20) and of 0.002 molal MgCl₂ at cells (60, 60) and (80, 80)
• Diffusion coefficients: isotropic homogeneous α = 1E-06
• Time steps & iterations: 10 iterations with Δ t of 200 s

• DHT parameters:

  H	O	Charge	C(4)	Ca	Cl	Mg	Calcite	Dolomite
  10	10	3	5	5	5	5	5	5

• Hooks: the following hooks are defined:

  1. dht_fill:
  2. dht_fuzz:
  3. interp_pre_func:
  4. interp_post_func:

dolo_interp_long.R

• Grid discretization: rectangular domain of 5 ⋅ 2.5 m² discretized in 400 × 200 cells
• Boundary conditions: Fixed concentrations equal to the initial state are imposed at all four sides of the domain. Fixed concentration of 0.001 molal of MgCl₂ is defined in the domain center at cell (100, 50)
• Diffusion coefficients: isotropic homogeneous α = 1E-06
• Time steps & iterations: 20000 iterations with Δ t of 200 s

• DHT parameters:

  H	O	Charge	C(4)	Ca	Cl	Mg	Calcite	Dolomite
  10	10	3	5	5	5	5	5	5

• Hooks: the following hooks are defined:

  1. dht_fill:
  2. dht_fuzz:
  3. interp_pre_func:
  4. interp_post_func:

References

• De Lucia, Kühn, Lindemann, Lübke, Schnor: POET (v0.1): speedup of many-core parallel reactive transport simulations with fast DHT lookups, Geosci. Model Dev., 14, 7391–7409, 2021. https://doi.org/10.5194/gmd-14-7391-2021
• Möller, Marco De Lucia: The impact of Mg²⁺ ions on equilibration of Mg-Ca carbonates in groundwater and brines, Geochemistry 80, 2020. https://doi.org/10.1016/j.chemer.2020.125611
• Palandri, Kharaka: A Compilation of Rate Parameters of Water-Mineral Interaction Kinetics for Application to Geochemical Modeling, Report 2004-1068, USGS, 2004.