Setting oxygen fugacity in a relatively complex reactive transport model

[Mattm]

Hello, 

I'm a grad student attempting to put together a reactive transport model in X1t. The goal is to construct a 2D reactive transport model to simulate groundwater chemical evolution in a substrate in my study area. My supervisor recommended starting in React and then working my way into X1t and then X2t, when I can get the previous simulations to run and follow the chemical trends I'm after.

My model includes clays, feldspars, carbonates, chlorite, and pyrite. I have XRD results to support wt% values for mineral species, and I'm using a series of kinetic reaction equations composed by my supervisor (which he gathered from reputable research). I got the model to run and follow the ion trends I have for 11 groundwater samples oriented along a flow path, using React. Now, I'm attempting to put together a 1D reactive transport model with the same setup. 

For my initial solution, I'm using a mixture of the 11 samples I have (seems like an okay idea), and I'm flushing into the model a groundwater sample that represents a relatively fresh (low TDS) groundwater sample from elsewhere in my dataset.

I got the model to run without pyrite with some okay results. My issue is that I can't get the model to run if I include kinetic pyrite (0.03 wt%). I think it's because I'm not starting with good initial redox conditions and the model can't converge because oxygen fugacity is shooting off somewhere, but I'm not sure. Inevitably X1t crashes, unless I start with an initial oxygen partial pressure of approximately 2 atm! I guess my foremost question is how you recommend setting initial redox conditions. The aquifer I'm attempting to simulate is overlain by glaciolacustrine sediments, often tens of metres thick, so dissolved oxygen content is low. But, the lower of a value I use for initial oxygen fugacity, the quicker the model crashes. 

I can get the model to run with "simple" pyrite, but I'm not sure if it is logical to keep adding pyrite at every step of the simulation.

One other thing I'm wondering about is how to plot saturation index plots in Gtplot from data in GSS. This seems like a trivial task, and I've put together dozens of scatter plots in Gtplot, but I can't get the same values I have in my GSS spreadsheet to plot in Gtplot. I have Log(Q/K) values in GSS that hover around saturation with respect to calcite, but when I plot them in Gtplot, no matter what style of axes I use, it shows SI values up to 10.  

Attached is my X1t file, GSS spreadsheet of the samples oriented along the east-west flow path I'm trying to simulate, and I'll include a picture of the trends I'm trying to simulate. I've been at this for a few weeks, so any help would be greatly appreciated. 

 

Reactive transport model 1D.x1t Groundbirch PV samples.gss

Edited March 18, 2021 by Mattm

[Jia Wang]

Hello Mattm,

I would suggest you start with just a single sample fluid and not a mixture of the fluid if that's how you initially set up your simulation in React. Once you are successful with that, then you might move to a composition of a mixed fluid. 

It seems like you are using customized rate laws for your kinetic reactions. I can't run your simulation directly because the basic scripts didn't attach with the input file. Were you able to include pyrite in your React run? I would suggest going back to React for troubleshooting this. The other thing to check is whether the program is able to calculate the equilibrium state based on your initial constraints, before any reaction takes place. You can check this by doing a 'go initial' run in React. If this doesn't complete successfully, then I suggest you going back to your basis pane to make sure that the constraints set are correct. Also, does X1t show you an error when you try to run your input file? Or, are you able to start the run but the progress is so slow that it seems like it is not progressing? 

As is right now, your x1t input file is considering all your redox reactions to be in equilibrium. The oxidation state is set by the concentration of O2(g) partial pressure you have designated in the system. The Eh from this is set for all redox couples in the simulation. Unless you specify otherwise, the program prescribes redox equilibrium between species of differing oxidation states. You can disable this for each redox reaction (e.g Fe++ / Fe+++) by decoupling the redox pair in Config -> Redox Couples... If you are unsure about the partial pressure but have reliable dissolved oxygen measure, you might be better off starting your model with the measurement instead of guessing. If you are interested in enabling redox disequilibrium in GWB, please see section 7.3 in the GWB Essentials user guide. 

With regards to plotting, note that the mineral saturation is shown as log (Q/K) in your GSS spreadsheet but is plotted as Q/K in Gtplot. If you change the unit in GSS to Q/K (right click on log and uncheck), then you will see the value will match what is plotted in Gtplot.

Hope this helps,
Jia Wang 

[Mattm]

Hi Jia, 

Thanks for the reply. It turns out I had a few settings set up kind of wrong. But, ultimately, I determined that the reason the software kept crashing (without reading any errors - just closing down) was because my CPU was slightly overclocked. I returned my CPU to factory settings, then the software now runs the models! Not sure how that works.

Cheers,