Estimating soil geochemistry

[Pak]

Dear All,

I would like to express a labarotary model in a similar way as [1]. The model is based on a solution block, a surface block with let's say 3 types of surfaces (illite as non-electrostatic model, Hydrous ferric oxides as diffuse layer model, and the WHAM or Tipping model VI as constant capacitance model) for doing so, I use the CD-Music option, and also an exchange block.

I have the data from the solution and some parameter to estimate the amount of exchange sites, surface sites, etc. After estimation, I put them on a Phreeqc file and equilibrate them with the solution.

If I want to keep the same pH and pE of the laboratory, I need to trick Phreeqc by creating a fix_pH and fix_pe phases, and also allowing to add or remove the equations of O2(g) and CaOH. Something like this:

Equilibrium_Phases
fix_pe -10.7478  O2(g)
fix_pH -6.6 CaOH

The problem is that I would like to run this soil using HP1, and if I use the trick of fixing pE and pH, they will be fix for the whole simulation, which I do not desire.

I have tried to run the outputed raw data, namely solution, exachange, surface, and equilibrium phases after simulation. Although, in this new run, I have removed the equilibrium phases, and as expected the results do change instead of be invariant. I wonder if somebody knows a way to keep the system unchanged, and ready to be used without the fixing pE and pH trick.




Thanks



Ref: [1] Hormann, V., and H. W. Fischer. 2013. 'Estimating the distribution of radionuclides in agricultural soils – Dependence on soil parameters', Journal of Environmental Radioactivity, 124: 278-86.

[dlparkhurst]

I do not know how HP1 works. Can you run a PHREEQC input file to establish initial conditions and simply use "DELETE equilibrium_phases n" to remove the  EQUILIBRIUM_PHASES definition? or simply define a new EQUILIBRIUM_PHASES n to replace the fix definitions once the initial conditions are created?

You seem to be able to use _RAW input? Which works? So, you could run PHREEQC, get all the definitions right (including removing the fix equilibrium_phases) and write the _RAW definitions for HP1?

Note, you should use Ca(OH)2, not CaOH.

[Pak]

Sorry for my late reply, I was not available during a while.

I am not sure if HP1 handles raw input I have not try it so far. I thought that maybe there is a way to go from raw input to classic Phreeqc input. Is there a way? (or a better way to approximate the values of raw to classic, at least)

I have some questions about your ideas to solve the problem. I think I do not get what you mean by defining new Equilibrium_Phases. Do you mean to create a new block of equilibrium phases in order to remove the old one? (If that is the case, so far, I have not think about it since it is unclear what are the minerals and abundances right now)

I think I am not using the command Delete correctly. If you look in the attached files, you can find one phreec file as txt and the database. It is not written in raw format. Here I try to get 2 dump files: 1) one of the equilibrated system with the fix pH and pe, and one after by removing the equilibrium phases, but in the dumping files, in both I have equilibrium_phases_raw. Although I can see that some of the values have changed, it is unclear to me why. I notice that in the second simulation (the simulation that should not have equilibrium phases), there is a sort of precipitation of the fixed phases, at least that is what I understand from the second dump file, since initially the amount of moles are 0 but at the end is 10, which I also do not get. Do you know why?

As always thank you for your help

[dlparkhurst]

There is way too much to deal with here.

Let's start with what solution you want to be in equilibrium with your CD-MUSIC surface. Here are three logical options. I have changed you input a little to define the iron as ferrous, which is consistent with your comment, to make a point. I have also used phreeqc.dat just because it is simpler.

Code: [Select]

PHASES  # from PHREEQC.dat
Fix_pH
H+ = H+
log_k 0.0
Fix_pe
e- = e-
log_k 0.0     
END
SOLUTION 1 Top Lysi 01-A
    temp      10
    pH        6.6
    pe        10.7478
    units     mg/l
    Fe(2)     0.0234
    O(0)      1 O2(g)      -0.68
    -water    1 # kg
END
MIX
1 1
SAVE solution 2
END
USE solution 1
EQUILIBRIUM_PHASES 1
fix_pe -10.7478  O2(g)
fix_pH -6.6 Ca(OH)2
SAVE solution 3
END

The script defines three solutions that differ in redox state. SOLUTION 1 is in redox disequilibrium--both Fe(2) and O2(aq) are present. Note that no saturation index is calculated for goethite because no Fe(3) is present. If you equilibrate your surface with this solution, then there will be Fe(2) sorbed.

Solution 2 is solution 1 after it has reacted to redox equilibrium. Fe(3) is now dominant and Fe(2) is negligible; dissolved oxygen is about 0.35 mmol/kgw. pH is 6.1. Goethite saturation index if 6.5. If you equilibrate your surface with this solution, Fe(3) species will be sorbed.

Solution 3 is solution 1 reacted to equilibrium with your fix_pH and fix_pe phases. Fe(3) is again dominant; dissolve oxygen is absent; pH is 6.6 as specified; and Goethite saturation index is 6.9. If you equilibrate your surface with this solution, Fe(3) species will be sorbed.

So, you have some choices. Do you think Fe(2) is present in your solution?, in which case you need a lower pe and absence of dissolved oxygen. Do you want dissolved oxygen in your solution, in which case your pe must be higher, as calculated from dissolved oxygen. If you stay with your pe and pH, then iron is Fe(3) and greatly supersaturated with respect to iron oxyhydroxides and dissolved oxygen is absent.

At the specified pH (6.6) and pe (10.7), all redox elements will be in their oxidized state (Fe(3), U(6), As(5), etc), but dissolved oxygen will be absent.

In general, I do not like to specify pe. You have fixed pH too, but if you do allow pH to vary, it is possible to be outside the stability field of water. Better to specify a partial pressure of O2(g).

Why don't you figure out whether HP1 will allow _RAW input before we discuss your options related to classic versus _RAW input.

[Pak]

I am not sure if Fe(2) is on the solution, since we got only the total amount of Fe, but I am sure of the values of pH and pe (assuming that the accuracy of the instrument is acceptable). So probably we have Fe(3) rather than Fe(2). Somehow is weird that there is not dissolved oxygen, since these values are from the upper part of the column, I had expected that there would some O2 dissolved gas, but as you have showed in equilibrium is not feasible.

I have run a simulation in Hp1 using Raw_ format and classic format, it seems to work, but it gives different results, so I might be done a mistake, or maybe there software right now is not though to do such thing. Anyway, if working with raw is better, I can always couple iphreeqc or phreeqcRM to a Comsol.

[dlparkhurst]

OK. The next item to consider is the use of END statements.

When you have

Code: [Select]

SOLUTION 1
...
EQUILIBRIUM_PHASE 1
...
EXCHANGE 1
-equil 1
...
SURFACE 1
-equil 1
...
END

The exchanger and surface are equilibrated with solution 1 before it has reacted with equilibrium phases and then the solution, exchanger, surface, and equilibrium phases are brought together and reacted. In particular, the  exchanger and surface would be equilibrated with solution 1 before the solution has been reacted to redox equilibrium. In your case, it will not make much difference (all elements in the oxidized state), but conceptually, it could make a big difference. In my previous post, using Fe(2) would mean that Fe(2) species equilibrate with the exchanger and surface, and then, in the reaction calculation (solution+exchange+surface+equilibrium phases) the Fe(2) in solution, exchange, and surface would be oxidized by addition of O2 from equilibrium phases.

You should also look at the Types and Sequence of Calculations section of the manual. If you separate keywords by END, it will make less difference, but if you have a long string of keywords without ENDs, then it becomes important.

The sequence of calculations is
(1) Initial solution,
(2) Initial exchange (-equil)
(3) Initial surface (-equi)
(4) Initial gas phase
(5) Batch reaction (SOLUTION + reactants)
(6) Inverse modeling
(7) ADVECTION
(8) TRANSPORT
(9) RUN_CELLS
(10) COPY
(11) DUMP
(12) DELETE

Therefore, when you have DELETE and DUMP in the same block between END statements, DUMP will occur before DELETE. If you have

Code: [Select]

DELETE ...
END
DUMP
...
END

then DELETE will occur before DUMP.

[Pak]

Yes, thanks for the insight on the preparation of the initial conditions, I did not pay that attention and did know that. And also thanks for my noob error, I forgot to take into account the sequence of calculation, now I can use the delete correctly. Thank you for your patience.

I have noticed that running the solution after removing the equilibrium phases have an effect on the values, for instance my pe goes from 10 (my desire value) to 15. Although, in general, most of the values stay close to the values of the simulation with the fix pH and pe trick. So probably, it is the closest that I can get to the real system which certainly is not in total equilibrium (and probably there are far more uncertanties).

I guess that in this case is better to use _raw format than classic format. I remember to ask time ago about going from _raw to classic format, and if I remember well the problem was on the correct amount of total H, OH, although I am not sure. If there other things that I should keep on mind, if I try to go from _raw to classic?

[dlparkhurst]

I don't want to go into a long exposition on redox. I will say that unlike acid-base reactions, redox reactions may be slow. The obvious example is that O2(aq) and N2(aq) coexist happily, even though NO3- is thermodynamically favored. In considering a natural environment, Berner had a categorization based on the presence of oxygen, nitrate, sulfide, and methane; this approach gets you in the redox ballpark without trying to determine "the" pe. A soil adds further complications because there may be oxic and anoxic micro-zones.

As for classic versus _RAW formats. The _RAW formats are probably easiest to deal with because they can be written and read directly. The following allows you to write a classic definition if for some reason that is preferable. Note that, as an example, I did not charge balance the initial solution. I suppose you could write a program to convert _RAW output to classic generating the same output as this script.

Code: [Select]

PHASES
Fix_pH
    H+ = H+
    log_k     0
Fix_pe
    e- = e-
    log_k     0
END
SOLUTION 1 Top Lysi 01-A
    temp      10
    pH        6.6
    pe        10.7478
    redox     pe
    units     mg/l
    density   1
    C         83.9
    Ca        68.4
    Cl        70.92
    Fe        0.0234
    K         23
    Mg        10.72
    N         121.6
    Na        16.19
    O(0)      1 O2(g)      -0.68
    S         4.3
    -water    1 # kg
END
USE solution 1
EQUILIBRIUM_PHASES 1
    fix_pH    -6.6 Ca(OH)2   10
    fix_pe    -10.7478 O2(g)     10
SAVE solution 1
SELECTED_OUTPUT 101
    -file                 classic.pqi
    -reset                false
USER_PUNCH 101
    -headings
    -start
 10  s$ =      "SOLUTION " + STR_F$(CELL_NO, 5, 0) + EOL$
 20  s$ = s$ + "-units mol/kgw" + EOL$
 30  s$ = s$ + "-water " + STR_E$(TOT("water"), 25, 14) + EOL$
 40  s$ = s$ + "-temp  " + STR_F$(TC, 11, 4) + EOL$
 50  s$ = s$ + "pH     " + STR_F$(-LA("H+"), 21, 14) + EOL$
 60  s$ = s$ + "pe     " + STR_F$(-LA("e-"), 21, 14) + EOL$
110 t = SYS("elements", count, name$, type$, moles)
120 FOR i = 1 to count
130   IF (name$(i) = "H" or name$(i) = "O") THEN GOTO 200
140   IF INSTR(name$(i),")") > 0 THEN GOTO 200
150   s$ = s$ + PAD(name$(i),7) + STR_E$(TOTMOL(name$(i)), 25, 14) + EOL$
200 NEXT i
300 PUNCH s$ + "END" + EOL$
    -end
END
SELECTED_OUTPUT 101
-active false
END
INCLUDE$ classic.pqi


[Pak]

Thanks