Questions about discretization artefacts in OGS 6.5.7–Phreeqc coupling

Help & Support
1

Dear OGS developers and users,

We are investigating reactive transport with the officially distributed OGS 6.5.7 binary. Our project uses ComponentTransport and selects the chemistry interface as follows:<chemical_system chemical_solver=“Phreeqc”>

According to the OGS source code and runtime messages, this selects the PhreeqcIO implementation. The official binary rejects chemical_solver="PhreeqcKernel" as not currently accessible.

Our model represents cementitious materials, bentonite-based material, an EDZ, and a groundwater mixing-cell boundary. The aqueous compositions can differ sharply across material and boundary interfaces. In simplified one-dimensional tests, we encountered two main numerical and state-handling issues, together with a related limitation concerning material-specific reaction sets. We would appreciate advice on whether these are known limitations, configuration issues, or avoidable by another supported approach.

1. Discretization artefacts in the transport–chemistry coupling

Our present understanding of the reaction step is as follows:

  1. Each element integration point is assigned a local chemical system.
  2. Nodal component concentrations are interpolated to the integration points.
  3. PhreeqcIO performs a chemical calculation for each local chemical system.
  4. For each transported component, the difference between the post-reaction and interpolated pre-reaction concentration is integrated with the shape functions and assembled into a nodal reaction right-hand side.
  5. A global equation containing the consistent mass matrix is solved to obtain the updated nodal concentrations.

For a linear two-node line element, the chemical calculation at either Gauss point samples an interpolated finite-element field containing contributions from both end nodes. Therefore, when the nodal compositions differ sharply, both integration-point chemical systems receive contributions from both sides of the element. For higher-order elements, some shape-function values at integration points can be negative, so the interpolated value can also fall outside the range of the nodal values.

In addition, a spatially local reaction contribution is distributed to multiple element nodes and then acted on by the inverse of the assembled consistent mass matrix. In our tests, this produces alternating positive and negative nodal concentration changes away from the reacting boundary.

We tested a strongly simplified one-dimensional case containing only the boundary-side shotcrete region and the mixing-cell boundary. Reducing the molecular diffusion coefficient from 1E-12 to 1E-30 m2/s did not materially change the early spatial oscillation of Na. At t = 3.835E-3 s, alternating changes were already visible over approximately 0.05 m; both calculations subsequently produced invalid input to Phreeqc at about t = 0.221147 s. This suggests that the observed early propagation is mainly associated with the discrete reaction mapping and global reaction solve rather than physical diffusion.

Could you please advise:

  • Is this behaviour a known limitation of the integration-point chemistry coupling used by ComponentTransport with PhreeqcIO in OGS 6.5.7?
  • Is there an existing option for mass lumping, a positivity-preserving treatment, a cell-centred formulation, or another recommended discretization for strongly contrasting aqueous compositions?
  • Would mesh refinement or a different element order merely reduce the artefact, or is there an established criterion or benchmark for selecting a stable mesh and time step in this situation?
  • Is a nodal or strictly cell-local chemistry coupling available or planned, in which a local aqueous state is reacted and returned without the integration-point-to-node reaction projection and consistent-mass inversion?
  • Is there another OGS-supported workflow that avoids these oscillations while retaining feedback of chemically induced porosity and transport-property changes?

2. Different reaction sets for different materials

Our material regions require different mineral sets, kinetic RATES definitions, and possible secondary phases. In our current PhreeqcIO project, these definitions appear to be common to all local chemical systems, while material dependence is represented mainly through zero or nonzero initial reactant amounts and material-dependent parameters.

Setting an initial amount to zero is not necessarily equivalent to removing a phase or reaction from a material-specific chemical system, particularly when precipitation of a secondary phase is possible. Is there a supported way in OGS 6.5.7, or in a later release, to assign genuinely different equilibrium phases, kinetic reactants, and RATES definitions to different material regions?

At present, the alternatives we can identify are:

  1. Use PHREEQC TRANSPORT with cell-based transport and chemistry, accepting that feedback of porosity, permeability, and diffusion changes to OGS would require separate handling.
  2. Modify and rebuild OGS to use a nodal or cell-local chemistry state, a positivity-preserving reaction projection or mass lumping, explicit transport of the required oxidation-state variables, and material-specific reaction sets.

Before adopting a custom implementation, we would like to ask whether there is a better existing method, an example project, a development branch, or recommended settings that we have overlooked.

We can provide a minimal reproducible project file, mesh, thermodynamic database, generated Phreeqc input, selected output, logs, and concentration profiles if useful.

Thank you very much for your advice.

Best regards,

Haruki Ishii
Consulting Department
Geosphere Environmental Technology Corp. (GETC)
URL: http://www.getc.co.jp