Issues on timestep sensitivity and output variables in reactive transport

Hello everyone,

I am currently reproducing the benchmark Transport and Cation Exchange with the original .prj file (no modifications). I encountered two issues and would greatly appreciate any insights:

1. Timestep sensitivity and negative concentration issue :

Non-reactive species (e.g., Cl) are insensitive to the time step, and the results agree well with PHREEQC. Reactive species involved in ion exchange (e.g., K) show strong sensitivity to the time step. However, when using small time steps, the simulation may fail due to negative concentrations, and the computation terminates. As a result, I am unable to reproduce the result shown in the benchmark.

Are there recommended strategies to improve stability? The three numerical stabilization methods either have no effect or result in completely incorrect results.

2. Output of sorption sites

Is it possible to output sorption sites (including ion exchange and surface complexation), similar to equilibrium_reactants, as cell-based or nodal values in the result files? It seems that currently this can only be achieved using “USER_PUNCH”.

Thank you very much for your help!