Difference between assigned and calculated Mass of water in PhreeqcRM

[sharmilan]

Dear Phreeqc exporters,

I am using Parallel Computation with PhreeqcRM. For the simple case, I have assigned only pure water (ph-7, temperature-20 C) without any other phases or gases. And set saturation-1, porosity-1, and representative volume-1. From the assigned details of the solution, the mass of water should be 1kg in each cell results, but in the cells result,  the mass of water is showing the value of 5.287e-01 kg. could anyone explain why there is a difference between assign value and the cell results?

here, I am attching assign solution result and cell result.

Thank you
sharmilan.

[dlparkhurst]

Not sure how to figure out what you have done. There are three types of IPhreeqc instances in PhreeqcRM--the workers, the initial phreeqc, and a utility.

The strategy is to define initial conditions by running a file (or string) in the initial phreeqc instance. If you define SOLUTION 1 in the initial phreeqc instance, I am sure you get the original result of 1 kg water.

If you then run your second RUN_CELLS simulation with -cell 1 with the initial phreeqc instance, I am sure you will still get 1 kg water.

For debugging, you should start with OpenMP, 1 worker, and a minimum number of cells.

After running the initial conditions in the initial phreeqc instance, you need to set the saturation, RV, density, and porosity. Next the strategy is to transfer initial conditions to the worker instances using InitialPhreeqc2Module. You should be able to use OpenFiles, SetPrintChemistryMask, and SetPrintChemistryOn, to turn on printing of the chemistry output file. At that point (after InitialPhreeqc2Module and before any transport), the PhreeqcRM method RunCells should produce the initial condition in the output file. If that is what you have done, and you do not get 1 kg water, I guess some serious debugging would be in order.

[sharmilan]

Dear Sir

Thank you for a quick response. I have found this problem in the initial phreeqc instance. here, I am attaching my Matlab script and output log file. I used 10 nodes (10cell) for my model. I think, now you will get everything in order from attached files (my Matlab script and output log file). could you explain where is a problem in assigned values or initail setting?

and

in my case, I have to include some mol of phases through REACTION keywords in the initial phreeqc instance. but there is no any option to include REACTION keyword in the initial phreeqc instance. but in PhrreqcRM, there is KINETICS option. So, could you give some ideas to how to set KINETICS and RATE data block for include exactly X mol of some phase in our system?

And

For second RUN_CELLS simulation, are we directly use REACTION keywords to include Y mol of phase, otherwise do we need to use MODIFY KINETIC keywords?

Thank you
sharmilan.

[dlparkhurst]

Dunno. I don't have access to Matlab, so I implemented your code in C++ in the attached advection_cpp.cpp and main.cpp. You can replace the files with the same names in Tests of the PhreeqcRM distribution to be able to run it.

The output file sharmilan.chem.txt shows mass of water of 0.9982 and a volume of 1 L using the database file phreeqc.dat. So, I don't know why you get mass of water 0.52. My only thought is that the order of 2D arrays ic may be different for Matlab arrays compared to C++ arrays.

Also attached is a source file for PhreeqcRM that may help with debugging. You need to replace the file in the src directory and recompile the library. It should provide some debug lines that show how the initial phreeqc solution is converted in InitialPhreeqc2Module. For the C++ version, a typical output is as follows:

Code: [Select]

ixyz (cell):         9
v:                   0.998203
f1:                  1
cell_porosity_local: 1
cell_saturation_local: 1
cell_rv_local:         1
current_v:             1.0018
n_old2:                -1


[sharmilan]

Dear Sir

Thank you for your quick response and effort, I have found this problem when using my database (CEMDATA18.dat) for my analysis. However, when I am using phreeqc.dat database, I also can get the same results as you mentioned in the previous mail (mass of water-0.998203 kg).  could you explain what is a problem in CEMDATA18.dat file and how to resolve this problem in CEMDATA18.dat (i have to use CEMDATA18.dat for my reactive transport model)? 
And also, could you give some idea for two more questions about REACTION which are mentioned in my previous mail?

here, I am attaching CEMDATA18.dat database for cement substances.

Thank you
Sharmilan.S

[dlparkhurst]

Change your definition of O in SOLUTION_MASTER_SPECIES to the following:

Code: [Select]

O H2O 0.0 O 15.9994

--In my case, I have to include some mol of phases through REACTION keywords in the initial phreeqc instance. but there is no any option to include REACTION keyword in the initial phreeqc instance. but in PhrreqcRM, there is KINETICS option. So, could you give some ideas to how to set KINETICS and RATE data block for include exactly X mol of some phase in our system?

(1) You could do the REACTION in the initial phreeqc instance and save the results to transfer to the workers:

Code: [Select]

USE solution 1
USE reaction 1
USE equilibrium_phases 1
SAVE solution 1
SAVE equilibrium_phases 1
END

(2) You could include a fake phase in EQUILIBRIUM_PHASES that has a very large log K and the amount of the reaction you want. It would dissolve the first time you ran RunCells in PhreeqcRM.

(3) You could have a fast KINETIC reaction. In the following, NaCl dissolves quickly and stops when it runs out of reactant (1 mole).

Code: [Select]

RATES
NaCl
10 k = 1
20 moles = k*TIME
30 SAVE moles
SOLUTION 1
KINETICS
NaCl
-formula NaCl
-M 1
END

If you want to add the reaction a second time, it depends on your conceptual model. Your KINETIC reaction could continuously add the reaction at a constant rate. If you are simply changing the solution at a boundary, you can do that in the initial phreeqc instance as a different SOLUTION and change your boundary to this solution at the proper time.

Otherwise, you could use EQUILIBRIUM_PHASES_MODIFY (option 2 above) or KINETICS_MODIFY (option 3). Note that if you are running in parallel, the cells will be split among multiple processors or processes, but I think if you generate the _MODIFY data blocks for all cells and run it with the workers (RunCells), it should work; I think the _MODIFY blocks for cells not present in a processor or process will be ignored. Alternatively, there is a way to get the individual worker instances and modify just the cells contained in each instance.

[dlparkhurst]

Or (4), you could define REACTION definitions for all of the cells to which you want the reaction added, followed by RunString for the workers. The next RunCells would include the reactions, but then you would need to do another RunString to delete the reactions if you do not want them to apply to subsequent executions of RunCells.

[sharmilan]

Dear sir

First of all, Thank you for noticing a problem in the database, it is a small mistake in my database, but it took two weeks of my time to solve this mass of water problem (mass of water is an important parameter because it makes the differences in the concentration of species between value through SELECTED_OUTPUT and RM_GetSpeciesConcentrations). Now, everything fine and my model is working well.

And, some more questions,

1) when I increase the value of the RM_SetRepresentativeVolume as 1000, we need to multiply mols by SOLN_VOL (k*TIME* SOLN_VOL) in the RATE data block. at the time I could not transfer exactly X mols through KINETICS data block. so, is there any method to include exactly X mols through KINETICS data block?

2) when I add the reaction a second time, do we need to assign REACTION data block of each cell separately with END? or assign REACTION data block together as MODIFY_keywords without END between each cell?

3) In each time step, I have to change surface area, surface site density and EDL water thickness base on the transported moisture level in each node. which is based way to modify or change those values without numerical convergence problems (more stable way)? 

3) in my reactive transport model, I have to add some EQ and SS phases in intermediate run steps. So, I have to add some more lines to previous SELECTED_OUTPUT. could you give an idea of how to add some more lines in the previously assigned SELECTED_OUTPUT data block or rewrite new SELECTED_OUTPUT data block?

Thank you
Sharmilan.S

[dlparkhurst]

(1) I would not increase the RV to 1000. Phreeqc works best with about a liter of water (plus or minus an order of magnitude). I think you may have convergence problems when using a large RV.

To specify an exact M for KINETICS, I would use an RV of 1.0 and SetUnitsKinetics = 0. In this case the moles of kinetic reaction per cell is M, as defined in KINETICS in the initial phreeqc instance, times RV. If you use a different RV, you can adjust accordingly.

(2) I think the processing of REACTION blocks would be faster if you do not use ENDs.

(3) I don't have much experience adjusting these parameters independently; I think you will just have to try it. There is an option to make the surface dependent on the amount of an equilibrium phase or a kinetic reactant. The surface area and number of sites are proportional to the amount of one of these reactants. Perhaps you could make use of this feature (use a ficticious kinetic reaction?).  I would try a simple SURFACE first, before adding too much complexity.

(3) Why can't you just define all of the needed fields in the initial SELECTED_OUTPUT/USER_PUNCH definitions? Any solid solution components not present in a cell will be set to zero in the selected-output file. You can use RunString with the workers to define another selected output file (new user number), but then you have to merge results from two files.

[sharmilan]

Dear Sir,

Thank you for your quick response,  for my case, I have to do the analyzes for a long time duration (nearly 50 years) and it takes too much of time. Due to the crashing of our reaction model at the intermediate step by assigned solution ions concentration and gas concentration (values from the transport model), is there any method to restart or reinitialize the reaction model by using dump or output file of the previous step of crashing step?  if it is possible, do we need to redefine the initial phreeqc instance or directly can we load the data to worker instance?

Thank you
Sharmilan.S

[dlparkhurst]

One note on one of my previous comments. There is a bug that will be fixed in the next version. If you link a surface to a kinetic reaction or an equilibrium phase and update the moles of reaction or phase with _MODIFY, the surface sites and area are not updated correctly.

Your reaction scheme sounds extra complicated. Do you have everything working correctly with simpler reactions before you add all of the complexity?

Yes, you can use DumpModule to save the state of the simulation. The simplest and safest way to get things started correctly again is to read the dump file with the initial phreeqc instance, and then use InitialPhreeqc2Module. There are other options if you are not using parallelization, or you do not allow load balancing with parallelization.

[dlparkhurst]

One additional comment. If you are updating the number of surface sites, and your surface is something like Hfo_wOH, it is important to update Hfo_w, O, and H rather than just Hfo_w. The following shows that adding just Hfo_w generates H2(aq) by the reaction:

Hfo_w + H2O = Hfo_wOH + 1/2H2(aq)

REACTION (and KINETICS) do not account for any charge on the reactants. Your _MODIFY blocks must explicitly account for a charge-balanced reactant. This kind of unwanted redox reaction often leads to convergence problems.

If you add Hfo_wOH in the REACTION, you get a more reasonable reaction.

Code: [Select]

SOLUTION 1
pH 5
P 1
END
SURFACE 1
Hfo_wOH 0.001 1 600
-eq 1
END
USE solution 1
USE surface 1
REACTION 1
Hfo_w 1
0.001
END

[sharmilan]

Dear Sir

Thank you for your reply and your ideas. in my transport model, the moisture, ions and gas transport are included. and, chemical equilibrium will be performed for transported mass (ions and gases). as you mentioned in the previous reply, I have adjusted the surface site mol through REACTION keyword and modify EDL thickness by using SURFACE _MODIFY keyword to avoid convergence problem. And also, transported gases and void volume through GAS_PHASE_MODIFY keyword. Mass of water is set by using RM_SetSaturation and RM_SetPorosity (RV=1).

but at 9th step, I have faced convergence problem, could you find out what is a problem in assigned value or in my the script? 

here, I am attaching output file of Phreeqc_RM ( i used 31 transport nodes (0-30cells) and set to print the results of cell 0 and 1 on output file).

Note: to set up initial phreeqc instance, I used some of the dump data and added additional data by using modify the keyword. And they are printed in the long list of lines in the output file, therefore, I have deleted those lines in the attached output files.

Thank you
Sharmilan. S

[dlparkhurst]

There are splitters and there are lumpers. You appear to be an extreme splitter. I lean to simplicity. Sorry, no way I can debug your system. It is way too complicated for me.

Just a couple comments.

Why do you have CO2(g) in both GAS_PHASE and EQUILIBRIUM_PHASES? Normally, only one is used. EQUILIBRIUM_PHASES is simpler mathematically.

The only "ERROR" (which indicates a calculation failed) is related to trying to finding the charge on a surface. I would remove your SURFACE and see if everything else works reasonably. Then add it back in slowly, say without donnan f, and without solid solutions, gas phase, and with a simplified EQUILIBRIUM_PHASES.

I note that your reaction adds Surf_s. As I said in the previous post, you should add Surf_sOH. Maybe that will make a difference or maybe not. You are generating sulfide, which I suspect is not reasonable.

I would also run just three or four cells with only one thread. With your cells split among 8 threads it is difficult to unravel the sequence.

[sharmilan]

Thank you for your response.

As you mentioned in the reply, for simplification of my chemical equilibrium model. I have run the system without gas phase and surface complexation, my chemical equilibrium model is working without any convergence problem. Then, I have added the gas phase to the system (without surface complexation) but chemical equilibrium model is again crashed with a convergence problem at 9th-time step. could you find out where is a problem in my chemical equilibrium model?

here, I am attaching output file (i used one thread)

thank you
S.Sharmilan

[dlparkhurst]

At pH 13 the partial pressure of CO2(g) is going to be negligible. My guess is that the numerical method has trouble with moles in the gas phase that are on the order of 1e-16.

So the GAS_PHASE, at least in this configuration, contributes nothing to the calculation, but makes it numerically difficult. I would simply remove it.

[sharmilan]

Thank you for your response.

Mol of gas is in the order of 1e-16 only in start_time 0 (that is before gas transport)  but in other time steps, transported mol of gas is in the order of 1e-3 (1.459e-03 mol). Therefore, I could not neglect or remove  GAS_PHASE in my chemical equilibrium model.

Thank you
S.Sharmilan

[dlparkhurst]

The amount remaining in the gas phase is 1e-16, so you could simply add CO2 to your reaction.

[sharmilan]

Dear Sir

Thank you for your reply.  for more simplification to avoid convergence problem, I have removed gas transport in Transport model and GAS_PHASE in my chemical equilibrium model. only dissolved CO2 ions and other ions in boundary cell 0 (dissolution of CO2 has performed through EQULIBRIUM_PHASE keyword in boundary cell 0) is set as a boundary condition for my transport model. but the same numerical issue is coming in the output file. Could you give any idea on this issue?

here, I am attaching output file (i used one thread)

Thank you
Sharmilan. S

[dlparkhurst]

Earlier, you said it ran without GAS_PHASE and SURFACE. Now that you have removed GAS_PHASE, and it does not work? What is different?

I guess I would start over, transporting without any reactions; add equilibrium_phases; add solid solutions, etc. Once it fails, try to figure out what it is that you added that causes the problem. Make sure the total concentrations of elements look reasonable at each step.

You can also try smaller time steps to see if it makes a difference in transport, plus it should make smaller changes per time step on the chemistry side. You can also try smaller values for -step (say 5) and -pe_step (say 2) in KNOBS.

Looking at the output, mass transfers of gypsum and other equilibrium phases are on the order of 0.1 mol or greater in the first several cells.  Seems like a lot for a 1 day time step.

[dlparkhurst]

I can't really debug your code. You could develop a transport problem in PHREEQC with your reactants. You can either use that to debug your code, or if it fails, it will give me something I can work with.

[sharmilan]

Dear Sir

Thank you for your reply.  when I change boundary ion concentration in the transport model, my chemical equilibrium model does not work after some step (i can be able to notice that when ettringite is dissolved completely (SI: -21.25), the chemical equilibrium gets fail).

in my transport model, there is moisture and ions transport model. but in phreeqc, there is no moisture transport case. Anyway, here, I am sending input data of cell 1 at the failed step (mass of water, transported ions concentration, and EQ and SS phase) and output file of PhreeqcRM (for the cell 1). by using the input value of cell 1 (at the failed step), could you run cell 1 to finding the problem in cell 1?

Thank you
S.Sharmilan

[dlparkhurst]

I can debug a PHREEQC input file, but not what you gave me. My guess is the problem is
in the second SOLID_SOLUTION, where one component is many orders of magnitude less
abundant than the other. Try putting the major component in EQUILBRIUM_PHASES and
ignoring the other component.

It does not appear that you have tried any of my suggestions to simplify to a system
that runs correctly, and then add complexity. You are not helpless here. You can
experiment removing and adding reactions and processes to determine where problems occur.

Best would be to start with a PHREEQC calculation and to simplify your model to include
only the processes in PHREEQC. If you are transporting individual ions, I assume it is a
diffusion calculation. While PHREEQC may not be able to transport moisture, it can
transport individual ions through multicomponent diffusion if Dw is defined for each
species in SOLUTION_SPECIES definitions (phreeqc.dat has definitions). It would be good
to compare a PHREEQC multicomponent diffusion calculation to your transport (without
moisture or reactions). How else are you going to convince people that your model is
working correctly? It looks like your transport does not produce net zero current 
because the charge imbalance increases with each step (see PHREEQC's method in the
mannual if you have questions).

After pure diffusion, you can add EQUILBRIUM_PHASES, SOLID_SOLUTION, and SURFACE in
sequence to see if all of those work in PHREEQC and compare well to your code.

Finally, you would debug your code for the features you are implementing that are not
available in PHREEQC, moisture, etc.

I don't think it is fatal, but your database does not have Fe(2) and Fe(3) defined correctly.
You have an Fe(3) species defined as the "master species", FeO2-. There should be only one
reaction that contains e- to define a "secondary master species" for Fe(2), probably Fe+2.
In your database, e- appears in multiple reactions for Fe(2) species: Fe+2, FeOH+, and
FeCO3(aq) for example. The secondary master species should only include Fe and O and H
if necessary. FeS will not work as a secondary master species for Fe.