What is the use of 'aquatic complexes'
I am wondering how aquatic complexes may change the simulation results.
below are my findings:
1. aquatic complexes, if declared in 'chemical.inp', are assumed to be at equilibrium with basic species. this means that without defining any aquatic complexes, user may be able to estimate the concentration of the aquatic complexes based on the concentration of basic species and the log K of the reactions.
2. aquatic complexes does not need to be claimed in the initial and boundary conditions, and their concentration during simulation can not be exported (so why should i declare them?).
3. I tried to run example 2 (P2_eos9_aquia) by removing the aqueous complex component, and got almost exactly the same results (with changes only on the fourth effective digits).
4. non of the minerals are reacting with aqueous complexes according to the chemical reaction table (therakin.dat), which means aqueous complexes will not change the dissolution and precipitations of the minerals.
I am running a TR simulation and found the simulation fail to converge (by checking the runlog.out) whenever i add any aqueous complexes (e.g., OH). Given that i can not check concentrations of complexes from the output, I do not need to input the initial values, why should i declare any aqueous complexes, and they will result in non-convergence?
I am using TOUGHREACT 2.0 any help on this will be greatly appreciated!
Thanks, Chenming
4 replies
-
Hi Chenming,
The "aqueous complexes" entry in chemical.inp allows the user to specify a specific set of secondary species that are present in the thermodynamic database. If the field is left blank, TOUGHREACT automatically takes all aqueous complexes from the thermodynamic database. If you look in the thermodynamic database, you will see that all secondary species are written in terms of reactions involving primary (basis) species. So they are always used, whether you declare them in chemical.inp or not. The reason to declare them in chemical.inp in some cases is to define the chemical system completely, without the code using many species that could be less important or not relevant for the conditions of interest. Hence, in your test problem you picked the important species, and obtained almost the same results.
The secondary species can be written out by specifying them in solute.inp. We don't write them out by default because there can be so many and the files would get very large. However, by specifying the ones you want to see in solute.inp, they will be written out after the primary species in the Tecplot file and in the time history file (e.g. time.out).
Regarding minerals (and gases) the reactions are also written in terms of primary species activities. However, that does not mean that the secondary species are not involved in the chemical system. The chemical system is solved completely with all primary and secondary species as well as minerals and gases. Before any mineral reactions are solved, the full aqueous chemical system is speciated (look at chump.out), and all the secondary species and their activities are shown in addition to the primary species activities, and the total concentrations for each component.
Hope this helps.
regards,
Eric
-
Chenming,
On your last question of why the simulation does not converge, there could be many reasons, so you should start by looking at the chdump.out, the runlog.out, iter.dat, and flow.out. It can start with problems in flow, leading to errors in transport, and then causing the chemical system to be unreasonable. Often it is just too large a time step. Reactive-transport needs well-controlled time-stepping, rather than the huge increases the flow can be solved with. So limit the timestep by using the Courant limit "rcour" in solute.inp (0.5 is a good value). Also make the maximum time step in flow.inp something reasonable for the problem, and always start with a very small time step, like 1 second. Always make sure you have good flow initial conditions at steady-state, and then run the flow problem so you know it is reasonable. Then run transport with a tracer without mineral reactions, to see that the transport is behaving as expected (the transport is more sensitive to the gridding and the timestep than the flow). Then add the mineral reactions. This work flow will make it much easier to figure out what could be the problem, and it will give much more confidence in the results.
Lastly, TOUGHREACT V2.0 was released in 2012, and we have since released 2 new parallel mostly rewritten versions with numerous new features, most recently V332-OMP in 2017, which are described on the LBNL website.
cheers,
Eric
-
Thank you Eric for your prompt response, appreciated! I will again look into it and get it back ( may take a few days, as I need to re examine the other output files)
-
Now i understood better, as shown in the manual:
""Aque-2. Secondary aqueous species
This input block can be omitted. In this case, all possible secondary aqueous species found in the input thermodynamic database will be automatically selected.""This suggests as a starter, it would be better to leave the secondary aqueous species empty in chemical.inp, which will consider all possible reactions related to second aqueous species and therefore not miss out any reactions listed in therakin.inp that could be potentially important, despite more computational expensive (which is ok for me now).
Following your suggestion, i do have found the way to output secondary aqueous species (i.e., NWAQ in Record_7, and Record_11, which input the indexes of the primary and secondary species that user wants to output). I also noticed that chdump.out lists all the indices of the species, which is quite convenient.
However, I am a bit confused about the units of the data output in OUTplot (shown at the lower section of the figure, filename is tec_conc.dat, sourced from example case P10_eos1_Metal-cycly).
From solute.inp (upper section of the figure), primary species 8, 15, 16,18 are specified both in record_9 and record_11. and therefore they are listed twice in the OUTplot file (bottom section of the figure), but in different numbers, which to me is supposed to be the same. (e.g., t_fe+2 is 0.2034e-18; while fe+2 is 0.1733e-18 for the whole column listed)
so what is the unit difference between t_fe+2 ouput by record_7 and fe+2 output by record_11?
Thanks, Chenming
