Chemistry Surface PK Time-Step Cut Issue in ATS-PFLOTRAN Reactive Transport

[Xueyuan Kang]

Hello all,

I’m trying to use ATS–PFLOTRAN to simulate watershed-scale mineral weathering. The reactive transport configuration (chemistry PKs + transport PKs) is taken largely from Xu et al. (2022, WRR ),  while my flow model is applied to a different watershed near Copper Creek. The model starts normally, but after ~0.5 simulation days it consistently fails during the chemistry surface PK step.

 

Below is the error message:

I enabled debug cell for cell 2520, and the logs show:

• Extremely low surface saturation (sat_new ≈ 7.46e-11)
• Zero initial concentrations (tcc_old = 0)
• Negative advective contribution for some species (cons(adv) < 0)

A snippet of the debug output:

Initially I suspected the issue was caused by cells drying out due to evapotranspiration, which could cause PFLOTRAN chemistry to fail (as mentioned in an earlier discussion thread). However:

1. Removing ET from surface-water_source does not fix the problem.
2. Adding a constant precipitation to the climate forcing does allow the model to run smoothly.

This suggests the failure is triggered by extremely low surface saturation combined with the reactive transport configuration from Xu et al. (2022)? but I’m not sure whether the problem lies in:

• the surface transport formulation under near-dry conditions,
• PFLOTRAN’s reaction time-stepping when concentrations are extremely low,
• or a mismatch between the chemistry PK setup and my hydrologic conditions.

I also tested a simpler example from the ATS demos (https://github.com/amanzi/ats-demos/blob/ats-demos-1.5/13_integrated_hydro_reactive_transport/integrated_hydro_reactive_transport-superslab.xml). That is, I used the chemistry configuration from the demo while keeping my own flow model. However, the simulation failed with the same issue. My ATS version is 1.5.

Has anyone encountered similar behavior, especially when using PFLOTRAN chemistry on nearly dry surface cells? Any suggestions on stabilizing the chemistry PK under low-saturation conditions would be greatly appreciated.

Thank you in advance!

Best,

Xueyuan

[Coon, Ethan]

What version of the code are you using?  Hopefully 1.6?  If you are using 1.5, please update to 1.6 first.  Please attach your input xml file.

Realize that “surface saturation” here is actually ponded depth, in units of m.  So that is order 1e-11 meters of water — frankly I’m surprised that the flow solution thinks that this is non-zero ponded depth.

Can you look at the same debug cell in the flow solve — is the pressure on the surface > atmospheric?

It may be as simple as setting a “threshold” on the water, so that we don’t call chemistry unless surface water content is greater than some threshold (rather than greater than 0).  I’m not entirely sure how that works without looking at the code, but I want to make sure you’re on 1.6 first before we start looking harder.

Ethan

From: ats-...@googlegroups.com <ats-...@googlegroups.com> on behalf of Xueyuan Kang <xueyuan...@gmail.com>
Date: Tuesday, December 9, 2025 at 10:08 AM
To: ats-...@googlegroups.com <ats-...@googlegroups.com>
Subject: [EXTERNAL] Chemistry Surface PK Time-Step Cut Issue in ATS-PFLOTRAN Reactive Transport

This Message Is From an External Sender

This email was sent from a non-ORNL address. If suspicious, use the Report Phish button in Outlook.

 

--
You received this message because you are subscribed to the Google Groups "Amanzi-ATS Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ats-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/ats-users/CAHHf1KSAJVWCXO1voSmNe811xD-UJwaQZyWgZWkXU7ZZ9q3_3w%40mail.gmail.com.

[Sergi Molins Rafa]

Yes, there is a saturation threshold in the chemistry pk but the default value works mostly fine. That would be a second step. 

I would try to update to 1.6, and run non reactive tracers before doing reactions. 

Sergi

To view this discussion visit https://groups.google.com/d/msgid/ats-users/SA0PR09MB6953865EC4BB4DEB7311D9C5A9A3A%40SA0PR09MB6953.namprd09.prod.outlook.com.

[Sergi Molins Rafa]

---------- Forwarded message ---------
From: Sergi Molins Rafa <smo...@lbl.gov>
Date: Tue, Dec 9, 2025 at 10:07 AM
Subject: Re: Chemistry Surface PK Time-Step Cut Issue in ATS-PFLOTRAN Reactive Transport
To: Xueyuan Kang <xueyuan...@gmail.com>

Hi Xueyuan

Did you run the same problem with one or two non-reactive tracers (still using PFLOTRAN) as explained in the Exercise 2a of the 2025 ATS Short Course (https://github.com/amanzi/ats-short-course/tree/ats-short-course-20250908/04_reactive_transport)? 

Please do so before running reactive transport. You could then try with a single tracer but using CrunchFlow as engine using the files available in the short course materials to rule out a PFLOTRAN-specific issue. Should not take much time if you're there already. 

The probable cause is the low saturation but we first need to figure out if this is too hard of a geochemical problem or just a numerical issue.

Thanks

Sergi

[Xueyuan Kang]

Hi Sergi and Ethan,

Thank you both for your replies. I am currently using ATS 1.5.
I checked the corresponding cell in the overland flow PK, and its surface pressure is only a little bit greater than atmospheric:

Regarding the idea of introducing a saturation (or pond depth) “threshold” to prevent chemistry from being called in essentially dry cells—would this be something that can be enabled through an XML parameter? Or would it require modifying and adding functionality in the source code?

I did attempt to update to 1.6 previously, but ran into issues during the geochemical engine installation. I will try again, but if there is a straightforward XML-level parameter controlling this, it would help me continue running the model in the meantime.

Thanks again.

Best,
Xueyuan

[Xueyuan Kang]

Hello,

Attached is my input XML file (v1.5). I am wondering whether there is an XML parameter that can control the saturation threshold (or ponded-depth threshold) for triggering surface chemistry.

In the meantime, I will also work on updating to version 1.6.

Thanks,
Xueyuan

[Coon, Ethan]

Please don’t do anything until you update to 1.6.  ATS 1.5 to 1.6 changed a lot of reactive transport, and fixed many bugs.

Ethan

From: ats-...@googlegroups.com <ats-...@googlegroups.com> on behalf of Xueyuan Kang <xueyuan...@gmail.com>
Date: Tuesday, December 9, 2025 at 5:53 PM
To: Amanzi-ATS Users <ats-...@googlegroups.com>
Subject: [EXTERNAL] Re: Chemistry Surface PK Time-Step Cut Issue in ATS-PFLOTRAN Reactive Transport

This Message Is From an External Sender

This email was sent from a non-ORNL address. If suspicious, use the Report Phish button in Outlook.

 

Hello,

To view this discussion visit https://groups.google.com/d/msgid/ats-users/7a40e1ee-dbd7-4487-a3a1-7fbbb8ba1b2cn%40googlegroups.com.

[Sam Shaheen]

This issue (which I also encountered in 1.5) seems fixed in 1.6

Setting "saturation tolerance" in the surface chemistry PK is the XML parameter Sergi is referring to. I found I basically had to turn surface chemistry off with this for things to converge in 1.5. No longer the case, so updating is the best solution!

[Xueyuan Kang]

Thank you all for your suggestions! 

Yes, I will try to update to 1.6 to see whether it can fix. 

[Xueyuan Kang]

Hello all,

I have updated ATS to version 1.6, and I encountered an input-related issue. When I use the PK tree shown below, the setup works fine in version 1.5, but it fails in 1.6.

PK tree:

The error message is:

Evaluator "water_flux" @ "reactive_transport_coupler_next" cannot be created in State.
Verify that (1) SetEvaluator is called, or (2) the name exists in state->evaluators.

 

Interestingly, if I change the PK type of the top-level “flow and transport” block from “subcycling MPC” to “weak MPC,” the error disappears and the model runs. However, it still exhibits the same “maximum number of reaction time step cuts” issue that I previously encountered.

 

I have two questions:

1. What is the difference between “subcycling MPC” and “weak MPC” in ATS?
   Will this choice lead to noticeably different simulation results?
2. To avoid the reaction-cut issue in ATS 1.6, do I still need to explicitly set a “saturation tolerance” for the surface chemistry PK?

 

I have attached my input XML files for the 1.6 setup.

Any guidance or suggestions would be greatly appreciated.

Best,

Xueyuan

 

You received this message because you are subscribed to a topic in the Google Groups "Amanzi-ATS Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/ats-users/OVx4oJwwdfA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to ats-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/ats-users/1ed387bf-40f8-4402-88d4-5c9caa9af3fbn%40googlegroups.com.

[Sam Shaheen]

To use the subcycling MPC with reactive transport, you just need to add the evaluators that are giving you bugs using an alias, as explained by Ethan here:

https://groups.google.com/g/ats-users/c/l6bF9CD9z54/m/FnxgH1I4AQAJ

[Xueyuan Kang]

Thank you, Sam.

After manually adding the evaluators, the model is able to proceed now. However, I ran into another issue related to the simulation start time.

Specifically, if I do not start at t = 0 but instead start at t = 20 years (after a flow-field spin-up), and use the checkpoint.h5 file from the spin-up run as the initial condition for flow, I encounter the following error (very similar to the one discussed here: https://groups.google.com/g/ats-users/c/AINbJY5-Lpw/m/-Iyix12RBwAJ):

terminate called after throwing an instance of 'DBC::Assertion' what(): Assertion: "t_this >= t_current" failed

This setup was working in ATS v1.5, but after moving to v1.6 it now fails. If I initialize flow from a checkpoint and start both flow and transport from t = 0, the simulation runs fine. Using the same configuration in a flow-only model also works smoothly. The issue only arises when the transport PK is included and the simulation start time is non-zero.

Is starting the simulation at t = 0 now required in this case? If so, I can try modifying the climate forcing data so that it effectively starts from t = 0.

Best,
Xueyuan

[Coon, Ethan]

You shouldn’t have to do that.  Instead use the “coupled flow and transport” MPC, which will add the evaluators for you (amongst other things).  That thread is from when 1.6 was in development.

Really you should run your file through the input converter — $ATS_SRC_DIR/tools/input_converters/xml-1.5-1.6.py.  This would make this change (and others) for you.  

Then compare to one of the regression tests for flow and reactive transport, e.g. 07_reactive_transport/column_infiltrat
ion_alquimia_calcite.xml which should help make sure your physics are all correctly updated now.

Once you’ve done that, if you’re still seeing the time issues, let’s talk more!  Starting at non-0 start time is fine.

Ethan

To view this discussion visit https://groups.google.com/d/msgid/ats-users/5eec63ce-6338-4aa1-9d14-ba0568efa7d2n%40googlegroups.com.