X2t seems to be running but get stuck

[iamseanty]

Hi there,

I am using X2t to do reactive transport modelling. My model runs successfully at 25C, 100% finished. 

But after I increase temperature to 80C, the model progress stays at 18.32% and won't proceed anymore... It's strange that the model doesn't crash and looks still busy running... Do you have any idea about this issue? 

Thanks and regards

Yuan

[Brian Farrell]

Hi Yuan,

You can report the factor limiting X2t's time steps by checking "explain steps" in the Config > Stepping dialog, and then "Follow Output" on the Results pane. This might give you a clue about why your run is failing. For more information, please see 9.40 explain_step in the GWB Command Reference. If that doesn't offer any good clues, you can attach your script and we'll try to take a look.

Regards,

Brian Farrell
Aqueous Solutions

[iamseanty]

Hi Brian,

Thanks for your reply. I checked the details as you suggested and I found that this is all because of the carbonates are consumed (very close to zero) at 80C at ~18% of the simulation, which forces the model to constrain time step in response of the change of mass of dolomite. 

Here I attach two simplified X2t input files for comparison. They are mostly the same except for the way of expression of dolomite dissolution kinetics: test1.x2t uses the default built-in rate law and test2.x2t uses a script, of which the rate formula looks the same as test1 input file.

The difference is that: test1 can run till the end of simulation without any issue; however, test2 runs to ~0.09% of the simulation and gets stuck then. That is when dolomite in node1 is close to zero... 

Do you have any idea of the cause of such discrepancy from the two input files? 

Thanks 

Yuan

test_1.x2t

test_2.x2t

[Brian Farrell]

Hi Yuan,

I don't have your custom thermo dataset, but I was able to check out your script using the installed thermo.com.V8.R6+.tdat. I think the difference between your custom rate law and the buildt-in rate law for Dolomite involves the surface area term. The surface area carried in the built-in rate law is the product of a specific surface area in cm2/g (set with the "surface" keyword) and the current mass in g of the mineral (you have set the "volume%" keyword, and the program converts this to g). In your custom rate law, though, you include only the specific surface area (a value of 530). You need to multiply that by the current mass of the mineral in your custom rate law. You can use the helper functions mw() and mass() to return the mole weight and mass in moles, respectively, of any species, so your rate law might include: *530*mass("Dolomite")*mw("Dolomite")*. Alternatively, you can use the internal parameters mw and mass within a rate law to refer to that specific mineral: *530*mass*mw*. If you want to include the parameter for the specific surface area, it might look like *sparea*mass*mw*. Or, if you want to jump straight to the actual surface area, just use the "surface" internal parameter: *surface*. Unfortunately, "surface" has different meanings in the context of the built-in rate law's keywords and the internal parameters for the custom rate laws. For more information, please see Table 5.1 and Table 5.2 in the GWB Reaction Modeling Guide.

Regards,

Brian

[iamseanty]

Hi Brian,

Thanks for your email. Yes, we used thermo.com.V8.R6+.tdat. 

We modified test2.x2t input file following your suggestions. And this time the model ends up with convergence issue at 0.095%, this is where dolomite is mostly consumed (~1.67e-62). So this is like a puzzle again why test1 can run without any convergence problem and test2 can't? 

I attached the modified test2.x2t for your consideration. Test1.x2t is unchanged. Could you please have a check? 

Thanks 

Yuan

test_2_mod.x2t

[iamseanty]

Hi Brian, if you have time, could you please have a look?

The functionality of custom-built kinetics are critical to our project.. and we would like to resolve this problem with your help. 

Thanks and regards

Yuan 

[Brian Farrell]

Hi Yuan,

I'm looking into it. Thanks for your patience. In the meantime, does the built-in rate law work for your purposes?

Regards,

Brian

[iamseanty]

Hi Brian, 

The built-in rate law doesn't work for our case, because the dissolution rates of the carbonates in our system are dependent on three promoting species with three sets of parameters (pre-exp, Ea and power)... we have to use script (or something similar, like equation or script file) to account for these effects. 

At 25C, the script works because the carbonates are not fully consumed in the system, and the model finished without any problem. When T was increased to 80C or 100C, sulfides and carbonates react a lot faster, causing carbonates running very low after several days and model crashed soon due to convergence issue during the simulation. The testing x2t input file just provided an example for this. 

We have to simulate higher temperature for this project, so your help to solve this problem is greatly appreciated. 

Thanks and regards

Yuan

[Brian Farrell]

Hi Yuan,

A numerical model isn’t guaranteed to handle everything that’s thrown at it. Your simulation with the built-in rate law runs to completion, but it still struggles quite a bit because the problem isn’t especially well posed. The dolomite reacts orders of magnitude faster than any of the other kinetic minerals. This forces really small time steps. And because the mineral dissolves completely, in a moving front as the model marches through time, you have incredibly high reaction rates adjacent to zero rates. 

Considering these substantial differences, there’s really no point in using kinetics for the dolomite. An equilibrium model should suffice. Your simulation should run faster, too, because the other kinetic minerals don’t require such small time steps. Or, if you really want to continue using kinetics, you should decrease the rate slightly. Decreasing the pre-exponential factor by even an order of magnitude lets the model run to completion. There’s essentially no difference in the output, considering the current course grid, which can accentuate differences. You could even set in your custom rate law a reasonable maximum reaction rate, and take the minimum of that and your calculated rate. The point is that there’s little difference in the output whether you use an amazingly fast or just a really fast rate for dolomite dissolution, apart from the numerical instabilities you introduce with the amazingly fast rate.

It might seem wrong to disregard or tinker with the kinetics here, but keep in mind that there are always a lot of variables that go into reactive transport models. Certain conditions make the models very sensitive to some parameters and insensitive to others. The precise value you set for the diffusion coefficient, for example, won’t change your output compared to using the default value because the term is insignificant wherever groundwater flows at an appreciable rate. Sometimes you  need to take a step back and determine what is most important.

By the way, the values you’ve set for Courant and Xstable are far too large. Their default values, 1, are the largest they should ever be. You should only set smaller values for these parameters as needed. In this simulation the kinetics of dolomite dissolution end up limiting time steps anyway, but in general you shouldn’t ever increase those values. For more information, please see the GWB Command Reference. 

Hope this helps,
Brian