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

2

Dear Haruki,

thank you for the detailed description of the problem.

Regarding the negative concentrations, unfortunately this is an issue that is difficult to completely avoid with the current FEM-based ComponentTransport formulation in OGS, particularly when sharp concentration gradients or strongly contrasting aqueous compositions are involved.

A more fundamental way to guarantee local conservation and better preserve positivity would be to use, for example, a finite-volume formulation or a mixed FEM formulation. We have been discussing such approaches internally, but at the moment they have not yet been implemented in OGS.

For your current model, I would first suggest trying the Flux Corrected Transport (FCT) stabilization scheme. It was implemented specifically to suppress numerical oscillations and negative concentrations in the ComponentTransport process.

There is an OGS benchmark showing how to activate and use FCT:

In the project file, it can be activated by adding

<numerical_stabilization>
FluxCorrectedTransport
</numerical_stabilization>

to the ComponentTransport process.

I would suggest trying this first with your simplified 1D model and checking whether it sufficiently suppresses the oscillations and negative concentrations you observed.

Best regards,
Haibing

3

Dear Haibing,

Thank you very much for your explanation and for suggesting Flux Corrected Transport (FCT). We agree that a finite-volume or mixed-FEM formulation would be a more fundamental approach to local conservation and positivity preservation.

Following your suggestion, we tested FCT using a simplified model derived directly from the OGS Calcite reactive-transport benchmark, calcite.prj. We selected this benchmark because it provides a much clearer demonstration of the issue than our full cement–groundwater model.

The main result is that, even for this simple benchmark-derived problem, FCT did not suppress the oscillation generated by the chemical-reaction update when advection was set to effectively zero. FCT produced a stable and monotonic solution when the same model was run without PHREEQC reactions. However, after Calcite/Dolomite equilibrium reactions were enabled, alternating positive and negative nodal concentrations appeared immediately after the PHREEQC reaction update, and the following PHREEQC calculation stopped because of negative concentrations.

To ensure that FCT was actually active, we used the Picard solution path rather than the Newton path. The model consisted of 100 linear elements with a uniform element length of 0.005 m. Advection was effectively zero, the dispersivities were zero, and the molecular diffusion coefficient of every transported component was 1.0E-10 m2/s. The aqueous composition at the fixed boundary differed from the initial porewater composition by approximately nine orders of magnitude for Mg and Cl.

We compared the following two cases under otherwise identical conditions:

  1. Component transport without PHREEQC reactions; and
  2. Component transport coupled to Calcite and Dolomite equilibrium reactions through the PhreeqcIO-based coupling.

For the transport-only version of the benchmark-derived model, FCT worked as expected. After 100 s, the Mg concentrations from the fixed boundary into the domain were:

x = 0.000 m: 1.000000E-3 mol/kgw
x = 0.005 m: 3.996813E-7 mol/kgw
x = 0.010 m: 1.607444E-10 mol/kgw
x = 0.015 m: 1.063847E-12 mol/kgw

The profile was non-negative and monotonic, and the transport-only calculation completed all 20 time steps up to 2000 s.

However, when the PHREEQC equilibrium reactions from the Calcite benchmark were enabled, an alternating nodal pattern appeared immediately after the first reaction update at 100 s:

x = 0.000 m:  1.000000E-3 mol/kgw
x = 0.005 m: -9.650113E-6 mol/kgw
x = 0.010 m:  2.693122E-6 mol/kgw
x = 0.015 m: -7.214489E-7 mol/kgw
x = 0.020 m:  1.934397E-7 mol/kgw
x = 0.025 m: -5.170395E-8 mol/kgw

At the next PHREEQC call, the calculation stopped with Negative concentration in system. The same major negative-concentration locations remained when the time step was reduced from 100 s to 1 s and then to 0.01 s.

Our current interpretation is that FCT successfully stabilizes the transport equation, but it is not applied to the separate chemical-reaction update used by the PhreeqcIO coupling in OGS 6.5.7/6.5.8.

From our reading of the implementation, the two stages are handled separately:

Transport stage
    -> assemble the transport M, K, and b
    -> apply FCT and the associated mass lumping

Chemical-reaction stage
    -> interpolate nodal concentrations to integration points
    -> perform local PHREEQC calculations
    -> assemble a new consistent reaction mass matrix and reaction RHS
    -> solve the global reaction-update equation

In particular, solveReactionEquation() assembles a separate reaction mass matrix through assembleReactionEquationConcrete() and solves the reaction-update system without calling the FCT procedure or applying mass lumping to that reaction mass matrix.

For constant porosity, the reaction update can be represented conceptually as

M_reac * (C_new - C_pre) / dt = b_reac

C_new = C_pre + dt * inverse(M_reac) * b_reac

Although interpolation with linear shape functions and the direct assembly of the reaction RHS use non-negative weights, inverse(M_reac) contains negative off-diagonal terms. A spatially localized PHREEQC reaction increment can therefore generate an alternating-sign nodal response. This is consistent with the numerical pattern above.

Consequently, the benchmark-derived test demonstrates that, in a model with zero or negligible advection and relatively small diffusion such as 1.0E-10 m2/s, FCT can remove the oscillation from the transport stage but does not remove the oscillation generated during the integration-point PHREEQC reaction and subsequent projection back to the nodes. In other words, applying FCT to calcite.prj does not resolve the chemical-update oscillation, even though it resolves the corresponding transport-only problem.

For this reason, it appears difficult to obtain a fundamental solution using the current integration-point, PhreeqcIO-based coupling in OGS 6.5.7/6.5.8 unless one of the following is added to the reaction-update stage:

  • mass lumping of the reaction mass matrix;
  • a conservative positivity-preserving limiter for the reaction increment;
  • a locally conservative finite-volume or mixed-FEM treatment; or
  • another conservative method that maps the reacted integration-point states back to the transport degrees of freedom without generating negative concentrations.

We have therefore been investigating the PhreeqcKernel approach, with one independent PHREEQC chemical system associated with each transport node. Our motivation is not that PhreeqcKernel automatically guarantees positivity, but that nodal chemistry can avoid the integration-point-to-node reaction projection and the separate consistent reaction-matrix solve that appear to generate the alternating pattern in the PhreeqcIO workflow. The aqueous totals transported by OGS can be passed directly to the nodal PHREEQC systems, and the reacted nodal totals can then be returned to the same transport degrees of freedom.

This nodal approach still requires careful verification of water-mass scaling, mineral amounts, charge and redox-state reconstruction, and global elemental conservation. Nevertheless, for our low-advection, diffusion-dominated cement–groundwater problem, it currently appears to be the more practical route than relying only on transport-side FCT with the existing PhreeqcIO reaction projection.

Could you please comment on the following points?

  1. Is our understanding correct that FCT in the current ComponentTransport implementation is applied only to the transport assembly and not to the separate equation assembled in solveReactionEquation()?
  2. Is there an existing option in OGS 6.5.7/6.5.8 to mass-lump or otherwise positivity-limit the chemical-reaction update matrix, independently of the transport matrix?
  3. Would a mass-lumped reaction update be considered a reasonable short-term modification, provided that elemental conservation is verified?
  4. From the OGS development perspective, is nodal chemistry based on PhreeqcKernel considered an appropriate alternative for avoiding the integration-point-to-node reaction projection, or are there known limitations that would make another approach preferable?

We would be happy to provide the simplified Calcite project files and the concentration profiles if they would be useful for reproducing the issue.

Best regards,
Haruki

4

Dear Haruki,

Thank you very much for your detailed report. Sorry for the delayed reply. I just came back from my 1-week vacation.

  1. Is our understanding correct that FCT in the current ComponentTransport implementation is applied only to the transport assembly and not to the separate equation assembled in solveReactionEquation()?

From your report, it looks like that FCT is not correctly activated for PhreeqC coupling, or not fully suppressing the oscillation. In my opinion this is not be the right behavior. I will need some deeper debugging to confirm this. Just a question: Have you tried OGS 6.5.8 with the same benchmark you tested? Does it deliver the save behavior?

  1. Is there an existing option in OGS 6.5.7/6.5.8 to mass-lump or otherwise positivity-limit the chemical-reaction update matrix, independently of the transport matrix?
  2. Would a mass-lumped reaction update be considered a reasonable short-term modification, provided that elemental conservation is verified?

Mass lumping feature can also be switched on. What it does is to artificially increase diffusion. But I guess it will not help you too much. My experience is that FCT delivers much stronger suppressing effect. Actually, FCT activates mass lumping internally. So I will more focus on the above issue and figure out why it is not working properly with the calcite problem.

  1. From the OGS development perspective, is nodal chemistry based on PhreeqcKernel considered an appropriate alternative for avoiding the integration-point-to-node reaction projection, or are there known limitations that would make another approach preferable?

What PhreeqcKernel do is to directly access the Phreeqc internal memory for data exchange, and avoid using the PhreeqcIO interface that is file stream based. This feature is not fully completed yet. I just started testing this feature. Also, it is intended for a faster data transfer, it will have no effect on the negative concentration problem that we are facing.

I will get back to you within 1~2 weeks after some further investigation on the problem.

Best,
Haibing

5

Dear Haibing,

Welcome back from your vacation, and thank you very much for taking the time to review our report and follow up!
Regarding your question about **OGS 6.5.8**:
Yes, we tested the exact same benchmark on both **OGS 6.5.7** and **OGS 6.5.8**, and confirmed that **OGS 6.5.8 delivers the identical behavior** (the spatial projection in the reaction assembly step still produces the local negative concentration).
We look forward to your further investigation and findings in the coming weeks. Please let us know if you need our benchmark input files or any additional test data.
Thank you again for your great support!

Best regards,
Haruki