Mass balance issue when using MINERAL_SCALE_FACTOR and AFFINITY_POWER

[Steve Benbow]

The attached case appears to show an issue with mass balance when MINERAL_SCALE_FACTOR and AFFINITY_POWER are used to set the mineral reaction rate.  Apologies in advance if it is a case of user error.  I'm using PFLOTRAN v6.0.1.

The case has five reactive minerals: Ca-montmorillonite, quartz, gypsum, dolomite and kaolinite.  The Ca-montmorillonite rate is set using PREFACTORs for neutral, acidic and alkaline conditions; the quartz and gypsum rates use PREFACTORs for neutral conditions; dolomite and kaolinite use simple rate constants.  The case is intended to simulate a single cell (no fluxes in/out), so in particular the total Ca content should be conserved.  The simulation period is 1e4y and the case runs in a few seconds.

The attached figure shows results for two cases: one where MINERAL_SCALE_FACTOR and AFFINITY_POWER are not set for the Ca-montmorillonite reaction, and one where they are set (both +ve values) by uncommenting lines 126-127 in the input file.  The plots show evolution of mineral vol fractions (top) and the corresponding evolution of the total element molalities in the pore fluid (bottom) for both cases.

I find that when MINERAL_SCALE_FACTOR and AFFINITY_POWER are not set for the Ca-montmorillonite rate, total Ca (solid + aqueous) is approximately conserved.  The conservation is not perfect - total mol of Ca in the box starts at ~906 mol and ends at ~873 mol - I would probably have expected tighter conservation of mass.  (I try to control it with VOLUME_FRACTION_CHANGE_GOVERNOR.)

When MINERAL_SCALE_FACTOR and AFFINITY_POWER are set for the Ca-montmorillonite rate (with the rest of the input identical) the final mass is ~3585 mol (starting mass is ~906 mol as before).  A lot of Ca seems to be gained in the system in the form of gypsum, with seemingly no corresponding Ca-montmorillonite dissolution and no change in aqueous Ca to compensate - this appears non-physical.  (See orange curve in the top-right plot attached)

I was unable to find a regression test involving MINERAL_SCALE_FACTOR and AFFINITY_POWER to help identify any mistakes that I might have made in the input, so apologies again if this is a user error.

Steve

[Hammond, Glenn E]

Steve,

I will investigate this. Can you send me your Python script for plotting? I assume it is Matplotlib….

Glenn

From: 'Steve Benbow' via pflotran-users <pflotra...@googlegroups.com>
Date: Monday, August 18, 2025 at 8:32 AM
To: pflotran-users <pflotra...@googlegroups.com>
Subject: [pflotran-users: 8531] Mass balance issue when using MINERAL_SCALE_FACTOR and AFFINITY_POWER

Check twice before you click! This email originated from outside PNNL.

--
You received this message because you are subscribed to the Google Groups "pflotran-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pflotran-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/pflotran-users/f4cf28c6-6812-4854-a0fb-ab33a13f129bn%40googlegroups.com.

[Hammond, Glenn E]

Steve,

I fear that our definition of mass balance may be misleading. The mass is not summed across phases. In other words, Ca++ in minerals is not summed. The mass of each mineral is summed.

Glenn

To view this discussion visit https://groups.google.com/d/msgid/pflotran-users/PH0PR09MB7836DEFD4ED868C6D583EC619A32A%40PH0PR09MB7836.namprd09.prod.outlook.com.

[Steve Benbow]

Hi Glenn,

Thanks for taking the time to investigate.  My Python script is attached.  It's called as follows:

    python3 plot_onecell.py <filename.h5> <cell index> <log or linear> <pattern> [<log or linear> <secondary_item>]

The script produces the plots and also exports corresponding plotted csv data to file.

The earlier plots of mineral VF and Ca speciation, with pH on secondary axis, can be produced with:

    python3 plot_onecell.py onecell2.h5 0 "log" "*VF*" "linear" "pH"
    python3 plot_onecell.py onecell2.h5 0 "log" "*Ca*[[Mm]]" "linear" "pH"
   
I've attached a Paraview plot of the mineral data for comparison with the earlier plots (for the case when MINERAL_SCALE_FACTOR and AFFINITY_POWER are set).  It looks similar so I think the Python is okay.

The attached xls uses the csv data exported by the Python script to calculate total Ca in the system.  Since it's one cell with no fluxes in/out I'd expect this to be constant, but the calculated total Ca increases from 1.94E+2 mol at the start of the simulation to 2.87E+3 mol at the end.  These numbers are different from the ones I quoted in the first post (PEBCAK! - they are now the same values as reported in the mas.dat file), but I think the conclusion that there appears to be a Ca mass imbalance is still the same.  This seems clear from the mineral plot, where there doesn't seem to be any Ca-montmorillonite dissolution to compensate for the gypsum.

Once again, probably a user error, but the input is quite minimal and I can't see it.  I appreciate you taking the time to look at it.

Steve

[Peter Lichtner]

Attached is a derivation showing that total mass is conserved.

Peter

[Hammond, Glenn E]

Thanks Peter. I should refactor to sum total mass of each primary species across all phases….

Glenn

From: Peter Lichtner <peter.l...@gmail.com>
Date: Friday, August 22, 2025 at 10:09 AM
To: Hammond, Glenn E <glenn....@pnnl.gov>
Cc: pflotran-users <pflotra...@googlegroups.com>
Subject: Re: [pflotran-users: 8535] Mass balance issue when using MINERAL_SCALE_FACTOR and AFFINITY_POWER

Check twice before you click! This email originated from outside PNNL.

Attached is a derivation showing that total mass is conserved.

Peter

To view this discussion visit https://groups.google.com/d/msgid/pflotran-users/PH0PR09MB78365DAFDCC32116D0F925F59A32A%40PH0PR09MB7836.namprd09.prod.outlook.com.

[Peter Lichtner]

Yes

...Peter <iPhone>

On Aug 22, 2025, at 11:20, Hammond, Glenn E <glenn....@pnnl.gov> wrote:



[Steven Benbow]

Hi Glenn,

I wasn't sure where you were proposing to refactor below, but just to be clear, the total mass that I calculate in the xls is the mass across all phases, which doesn't seem to be conserved - but this is the same as what is reported in mas.dat.  So I think that mas.dat is already summing over all phases.  The attached plot confirms this, I think.

The quantity reported as Total_Ca++ in the hdf5 (and in the python plots) is the total Ca molality in solution, which I think is what is most useful to report there.

So in terms of reported quantities I'm not sure that anything needs to be refactored.  The current summing seems correct to me.  But I do seem to have a mass balance issue, which I guess could be caused by summing in the mass balance calculations (but seems unlikely?  my suspicion is still user error...).

Steve

On 22/08/2025 18:20, 'Hammond, Glenn E' via pflotran-users wrote:

Caution: This sender is from outside of Quintessa. Do not click links or open attachments unless you recognise the sender and know the content is safe.

You received this message because you are subscribed to a topic in the Google Groups "pflotran-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/pflotran-users/naBLNykFOhs/unsubscribe.
To unsubscribe from this group and all its topics, send an email to pflotran-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/pflotran-users/PH0PR09MB7836454D110793F7FA0DE8F09A3DA%40PH0PR09MB7836.namprd09.prod.outlook.com.

-- 
Dr Steven J Benbow BSc MSc PhD CMath MIMA
Quintessa Ltd, The Hub, 14 Station Road, Henley-on-Thames, Oxfordshire, RG9 1AY, UK
Tel: 01491 636246  DD: 01491 630051  Web: http://www.quintessa.org

Quintessa Limited is an employee-owned company registered in England, Number 3716623.
Registered office: Quintessa Ltd, The Hub, 14 Station Road, Henley-on-Thames, Oxfordshire, RG9 1AY, UK

If you have received this e-mail in error, please notify pri...@quintessa.org and delete it from your systems.  Our privacy and other policies are available from https://www.quintessa.org/about-us/policies

[Hammond, Glenn E]

Steve,

Yes, I also confirmed that the calculation of total mass is correct. I could tell from closer inspection of the code, but I also hand calculated just to be sure. The code follows in case anyone cares:

https://bitbucket.org/pflotran/pflotran/src/e390943cf878b055687e77ae24dc319c4ce995f7/src/pflotran/reactive_transport.F90#lines-663

BTW, it is easier to hand calculate mass using molarity instead of molality as you do not need liquid density in the calculation.

The mass balance issue involves the explicit update of mineral volume fractions at the end of the time step. If any mineral exists (e.g., less than a single molecule), the mineral can still dissolve beyond the amount (mass) that exists. For extremely fast rates, this is problematic.

I note that the time step governor based on change in mineral volume fraction is broken (you will see that dvfmx and dvf/dt are both zero). When I fix the governor and set it to a very low value, the discrepancy gets much smaller.

I confirmed that it is not an issue with using MINERAL_SCALE_FACTOR or AFFINITY_POWER. When I set both to 1, the results are nearly identical to the run that does not employ them. And as I ramp up both to the values that you use, the discrepancy increases (faster rates). When I set the governor to 1.d-7 (in the attached) in the fixed code,

https://bitbucket.org/pflotran/pflotran/branch/glenn/fix-vf-governor

the results are much closer. This tells me that it is an issue with the faster rate  (due to scale factor and affinity power), large time step and broken vol frac governor. 

Regarding the governor, I am working on a fix that will not break the regression tests.

Glenn

To view this discussion visit https://groups.google.com/d/msgid/pflotran-users/7c4583c3-e0a8-41c1-af48-9b0a020a10e4%40quintessa.org.

[Steven Benbow]

Hi Glenn,

It hadn't occurred to me that the vf governor might have been the cause
of the issue.  Thanks again for looking into it and working on the fix. 
Very much appreciated.

Regarding your point "The mass balance issue involves the explicit

update of mineral volume fractions at the end of the time step. If any
mineral exists (e.g., less than a single molecule), the mineral can
still dissolve beyond the amount (mass) that exists. For extremely fast

rates, this is problematic.", have you considered adding an option to
scale the rate for dissolving minerals so that the small amount of
mineral would take the whole timestep to dissolve to zero in order to
prevent dissolving beyond the mass that exists?  I.e. you could
calculate DT' = min(amount/rate) across all dissolving minerals, and
then scale the rates for the dissolving minerals by DT'/DT in the cells
where it is necessary.  Obviously this would introduce some inaccuracy,
and wouldn't be as accurate as taking smaller timesteps to better
capture the exact time at which mass tends to zero, but might be
sufficiently accurate when trading off with the ability to use longer
timesteps.  Just a thought.  I wouldn't be surprised if there are other
complications that make this impractical.

Steve