ChemPlugin Time Marching Loop

[elg0086]

I have a question about time stepping in ChemPlugin. The example scripts provided on the Chemplugin website use the "c.reportTimeStep" command but my script contains time-dependent conditional statements within a WHILE loop. In order to have better control over the time step, "deltat" is set to a fixed integer. As I understand it, the c.reportTimeStep command will use time step criteria established for an instance (i.e. "delxi," "step increase," etc..) to determine an appropriate step. However, based on the explanation of the courant number in the user guide:

 

"The Courant condition requires that a time step not exceed the time required to displace all the fluid from a ChemPlugin instance"

 

Following this rule as a guide, is it appropriate to define a time step as a fixed integer (i.e. 4000 seconds)?

while True:
    deltat = 4000

if current_time >= acid_start and current_time < close_time
    gw_link.FlowRate(0, "m3/s")
    gw_link.Transmissivity(0, "m3/s")
    acid_link.FlowRate(discharge, "m3/s")
    acid_link.Transmissivity(trans, "m3/s")
for c in cp:
    deltat = (deltat)
for c in cp:
     if c.AdvanceTimeStep(deltat):end_run(0)
for c in cp:
    if c.AdvanceTransport():end_run(-1)
for c in cp:
    if c.AdvanceChemical():end_run(-1)

 

Regards,

Erik

[Brian Farrell]

Hi Erik,

If you were writing a simple model that considered only advective transport of a non-reacting solute, you might figure the limiting time step from only the Courant condition. And if the velocity didn’t change with time, you could hard-code that value into your time marching loop. When you construct a reactive transport model that accounts for a variety of other processes, though, such as diffusion, heat transfer, and kinetic reactions, you need to account for the stability of the solution to each of their governing equations. ChemPlugin’s ReportTimeStep() member function in fact does this. The X1t and X2t programs use similar logic, which is described in section 2.20 Time marching in the GWB Reactive Transport Modeling Guide.

In your program, it might be ok to set your own time step without querying ChemPlugin, as long as it’s smaller than the largest possible time step that ChemPlugin would allow. But really, you should always have ChemPlugin instances report the stable time step. You can compare with your own desired time step, if you’d like, then use the lowest of the values.

Based on a previous conversation, it sounds like you might be worried that kinetic reactions are making the simulation take too long. Simply ignoring the reported stability limits would be a bad idea. Instead, you should consider whether you really need so many kinetic reactions, especially really, really fast kinetic reactions. Use kinetics for the more slowly reacting minerals, and equilibrium for the others.

Regards,

Brian Farrell
Aqueous Solutions