Questions about discretization artefacts and preservation of redox state 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. Preservation of redox state when total-element components are transported

Our current project transports Fe, S, and C as total-element components, rather than transporting oxidation-state-specific components such as Fe(2), Fe(3), S(-2), S(4), S(6), C(-4), and C(4) separately.

From our reading of the OGS 6.5.7 PhreeqcIO implementation, the configured component amounts, pH, and pe are stored for each local chemical system. After a Phreeqc calculation, selected output is read back into these fields. Before the next reaction calculation, the transported component values are interpolated to the integration points and used to construct the next Phreeqc input.

Our concern is therefore not that all Phreeqc state is necessarily discarded immediately after a chemistry calculation. Rather, when only total-element components are used as the OGS transport variables, the reacted solution is projected onto totals such as total Fe, total S, and total C for transport. The oxidation-state-specific aqueous inventories are then not independent transported degrees of freedom.

This may be acceptable when instantaneous redox equilibrium is intended and the oxidation-state distribution is uniquely reconstructed from the imposed equilibrium conditions. However, it may be important for kinetic reactions. For example, two solutions can have the same total Fe but very different Fe(2) and Fe(3) totals. A RATES expression depending on TOT("Fe(2)"), a redox-sensitive activity, pe, or a saturation index could then give a very different rate after the transported total has been used to reconstruct the next solution.

We have observed Phreeqc convergence or mass-balance errors involving C, S, and Fe in the full cement–groundwater reaction system. We do not yet know whether these errors are caused by loss or inconsistent reconstruction of redox-state information, by the negative/non-finite concentrations generated by the discretization issue above, or by a combination of both.

Could you please clarify:

  • In the official OGS 6.5.7 PhreeqcIO implementation, what is the intended state passed from one reaction step to the next when the configured transport components are total Fe, S, and C?
  • How is the pe value read from one Phreeqc calculation used when the transported totals are mixed or changed before the next calculation? Is pe a retained local chemical-state variable, an imposed value, or only an output unless fixing_pe is enabled?
  • Does PhreeqcIO officially support transporting oxidation-state-specific components such as Fe(2), Fe(3), S(-2), S(4), S(6), C(-4), and C(4) as separate OGS process variables and returning their respective valence-state totals consistently?
  • If oxidation-state-specific components are used, should the corresponding total-element component, for example Fe, be omitted to avoid redundant or inconsistent constraints?
  • What component basis is recommended for redox-sensitive kinetic transport: oxidation-state-specific totals, total elements plus pe, alkalinity, an electron-equivalent balance, or another formulation?
  • Are there examples, known fixes, or later developments that demonstrate robust coupling for redox-sensitive RATES expressions?

3. 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