Galena Kinetics

[Scarlet]

 I Am working with kinetic dissolution of Galena,They would first equilibrate with rain water,Can some take a look and see if I got wrong?

SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg

RATES
    Galena
-start
10 a0 = PARM(1)
20 v = PARM(2)
30 rate = (a0/v)*(M/M0)^0.67*10^-13.7*(1-SR("Galena"))
40 SAVE rate*TIME
-end

KINETICS 1
Galena
    -formula  PbS  1
    -m        1
    -m0       0.19
    -parms    8.164 0.584
    -tol      1e-08
-steps       1577000000 in 100 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5

INCREMENTAL_REACTIONS True

USER_GRAPH 1
    -headings               time Pb
    -axis_titles            "years" "mmol PbS /l" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3.1536e7
20 GRAPH_Y 1000*TOT("Pb")

[dlparkhurst]

You may have a problem initially because SR("Galena") will not be defined for pure water. You can either adjust your rate so that it is finite in pure water or add a little lead and sulfate to your initial solution.

[Scarlet]

Thank you , but This galena is part of an arkose that is balanced firstly with rainwater.
SOLUTION 1
    temp      25
    pH        4.5
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    Adularia  0 0.17
    Barite    0 0.17
    Cerussite 0 0.16
    Galena    0 0.19
    Kaolinite 0 0.08
    Kmica     0 0.07
    Quartz    0 0.17
RATES
    Galena
-start
10 a0 = PARM(1)
20 v = PARM(2)
30 rate = (a0/v)*(M/M0)^0.67*10^-13.7*(1-SR("Galena"))
40 SAVE rate*TIME
-end

KINETICS 1
Galena
    -formula  PbS  1
    -m        1
    -m0       0.19
    -parms    8.164 0.584
    -tol      1e-08
-steps       1577000000 in 10 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5

INCREMENTAL_REACTIONS True

USER_GRAPH 1
    -headings               time Pb
    -axis_titles            "years" "mmol PbS /l" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3.1536e7
20 GRAPH_Y 1000*TOT("Pb")
  -end
    -active                 true

[dlparkhurst]

I only meant a tiny amount so that the SR("Galena") would be defined, but I think it probably does not matter.

But, explain what you think will happen with galena as both a KINETIC and an EQUILIBRIUM_PHASE.

Also consider that you are plotting total Pb, but have stipulated that "Cerussite" ("Cerrusite" in current databases). So total dissolved Pb will not be the total amount of galena dissolved, but the combination of galena and cerrusite.

[Scarlet]

Thank you Mr David,

But i have an environment where the ph rain is neutral,and i need to study chemical speciation of heavy metal ( Pb,Cu,Zn...)  in pit lakes.
The arkose wich I chose contains (galena, cerusite ...)  after dissolution in pure water it gives heavy metals wich are present in lakes.
in this exemple i need to represent galena kinetic as a fonction of time.

[dlparkhurst]

You did not explain what would happen if you include galena in both EQUILIBRIUM_PHASES and KINETICS.

[Scarlet]

When I run the program, I have an increase in the pH solution as a function of time,and variation of molality. I want to know if I have not made errors when I wrote the program,or I must add the values of the (cerussite, kaolinite ...) in RATES and KINETIC ?


SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    Adularia  0 0.17
    Barite    0 0.17
    Cerussite 0 0.16
    Galena    0 0.19
    Kaolinite 0 0.08
    Kmica     0 0.07
    Quartz    0 0.17
RATES
    Galena
-start
10 a0 = PARM(1)
20 v = PARM(2)
30 rate = (a0/v)*(M/M0)^0.67*10^-13.7*(1-SR("Galena"))
40 SAVE rate*TIME
-end

KINETICS 1
Galena
    -formula  PbS  1
    -m        1
    -m0       0.19
    -parms    8.164 0.584
    -tol      1e-08
-steps       1577000000 in 10 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5

INCREMENTAL_REACTIONS True

USER_GRAPH 1
    -headings               time Pb
    -axis_titles            "years" "mmol PbS /l" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3.1536e7
20 GRAPH_Y 1000*TOT("Pb")
  -end
    -active                 true

[dlparkhurst]

Think about what will happen when you have galena in both EQUILIBRIUM_PHASES and KINETICS?

[Scarlet]

when I run program ( EQUILIBRIUM_PHASES and KINETICS):

Galena in equilibrium
Saturation index equal to zero
pH Increase in pit  lakes from 7 to 10

[dlparkhurst]

Yes, adding KINETICS does not have an effect because the solution will be in equilibrium with galena because of equilibrium phases.

If you want galena to react kinetically, then you should remove galena from EQUILIBRIUM_PHASES. Similarly for the rest of the phases, you must choose whether you want them to react kinetically (KINETICS) or react instantaneously to equilbrium (EQUILIBRIUM_PHASES).

[Scarlet]

Okay! Thank you Mr. David
Please! Is the "porosity" parameter of the galena considered 0?

[dlparkhurst]

There are porosity parameters in special TRANSPORT calculations, but they would not be zero. There is the area to volume ratio in your galena rate expression, but it would not be zero. So, I which porosity parameter you are talking about?

[Scarlet]

Thank you Mr david,
 I calculate the porosity parameter to delimit the 2 parameters A0 and V, in order to draw a curve of the development of galena kinetics as a function of time.
can I consider the index of saturation of the galena equal to zero when the ph=7?

[rfembilejr]

Hi,

My question is related to this so I will just reply to the old thread.
 
I am modeling galena dissolution in seawater using a batch reactor. We have obtained a rate law and we want to incorporate it in the RATES and KINETICS block. The plot of the code below is pretty much what happened in our actual experiment (pH was stable at 4.0). The reaction time is only for 5 hours.
However, I am not quite sure whether I should start with a seawater devoid of Pb, and how to input dissolved oxygen as 3.5 mg/L.  I also want to see how dissolved oxygen evolves. Ideally, it should be constant ie. 3.5 mg/L at the end of the simulation. Also, I would like to see at what time period will anglesite precipitate. I run a simulation for longer time and it shows anglesite will start to precipitate (SI >0) after four days, however, the parameters may not be all correct.
*The initial moles and surface area of 1.0 were modified to come up with the resulting plot. This value is not correct.

Thanks in advance for any advice.

### Galena dissolution in seawater
### Initial galena solid is 1 g
### Initial conditions are pH = 4.0, DO = 3.5 mg/L, Temp = 20
### Use Gn 36 as example

SOLUTION 1 Final solution (Nordstrom et al. 1979 seawater)
     units   ppm
        pH      4
        #pe      4
        density 1.03
        temp    20.0
        redox   O(0)/O(-2)
        Ca              412.3
        Mg              1291.8
        Na              10768.0
        K               399.1
        Fe              0.002
        Si              4.28
        Cl              19353.0
        Pb              1 ppb  ##this is the background, 14 ppb is the final concentration
        Alkalinity      141.682 as HCO3
        S(6)            2712.0
        O(0)            3.5     O2(g) 3.5  # not sure
   -water 1.03

EQUILIBRIUM_PHASES 1
     Anglesite 0 0


RATES

Galena
-start
1 rem PARM(1) = log10(A/V, 1/dm)
2 rem PARM(2) = exp for (m/mo)
3 rem PARM(3) = exp for H+
4 rem PARM(4) = exp for O2
5 rem PARM(5) = exp for Cl-
10 if (m <= 0) then goto 200
20 if (SI("Galena") >= 0) then goto 200
30 lograte = -11.18 + PARM(1) + PARM(2)*log10(m/m0) + PARM(3)*lm("H+") + PARM(4)*lm("O2") + PARM(5)*lm("Cl-")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 SAVE moles
-end


KINETICS 1

Galena
   -formula PbS
      -m 1                   #0.98 #from calculations similar to that of albite
      -m0 1
      -PARMS 1 0.66 0.41 0.41 0.26   
                        ##0.20 = first parm is from surface area, assumes full liberation
                     ##second parm is assumed to be 2/3 for uniformly dissolving spheres and cubes (Parkhurst & Appelo,1999)
                     ##0.42 is the reaction order of H+ from own calculation
                      ## 1.20 is the reaction order of O2 from own calculation
                     ## 0.34 is the reaction order of Cl- from own calculation
-steps       18000 in 14 steps # 5 hours
-step_divide 1
-runge_kutta 6
-bad_step_max 10000000

INCREMENTAL_REACTIONS True

USER_GRAPH 1
    -headings               Pb pH
    -axis_titles            "Time (Seconds)" "Pb (ppb) ""
    -initial_solutions      false
    -axis_scale x_axis      auto auto auto auto
    -connect_simulations    true
    -plot_concentration_vs  time
  -start
10 PLOT_XY TOTAL_TIME/3600, TOT("Pb")*207.2*1000*1000
20 PLOT_XY TOTAL_TIME/3600, -LA("H+")
1000 REM end
  -end
    -active                 true


USER_PUNCH 1
    -headings Pb_ppb Sulfur_ppb
    -start
10 PUNCH TOT("Pb")*207.2*1000*1000
-end

[rfembilejr]

If I may add to the previous post, how do I write 3.5 mg/L of dissolved oxygen in the input parameters block?

Thanks.

[dlparkhurst]

O(0)            3.5     O2(g) 3.5  # not sure

This input asks for a partial pressure of 2000 atm of O2(g), which gives a very large concentration of O2(aq)

O(0) 3.5

will give you 3.5 ppm O2(aq)

[rfembilejr]

Hi,

Following my previous post. I have a few questions regarding my model below. Initial conditions are specified in the heading. The rates equation is based on our findings (see below), we are testing it if it can simulate our observed Pb trend for seven hours. It seems that it is not replicating our actual results. The Pb concentration is way too low than our observed values. I was thinking it might be some mistake in the parameter values or just the code has some errors. We are also concerned with the saturation of anglesite that is why it is in one of the user graphs. The actual data is attached.

Rgalena = -10^-11.18 *(mH+)^0.41 * (mO2)^0.41 * (mCl-)^0.26

Another rate which we are also testing is: Rgalena = -10^-11.18 *(mH+)^0.41 * (mO2)^1.20 * (mCl-)^0.34
This one was obtained from a modified method.

Note: This experiment was conducted in a batch reactor with seawater circulating over the sample of 1 g galena.

Q: I set the pH of seawater at 4.0. This was done by adding HCl in the solution. I was wondering if this is ok or should I add a <Fix_pH -4.0 HCl 10> block to set the pH to 4.0. The seawater was prepared to have a pH of 4.0 prior to start of run.

Q: I was curious if I put the initial moles of 0.00418 correct. Calculation is from a 1 g of galena sample.

Q: The SA is based on the BET where the raw data is 0.2 m2/g. Do I need to adjust this?


Thank you in advance.

Rodrigo
--------------------------------------------------------------------------
### Galena dissolution in seawater
### Initial galena solid is 1 g
### Initial conditions are pH = 4.0, DO = 3.5 mg/L, Temp = 20

RATES

Galena
-start
1 rem PARM(1) = log10(A/V, 1/dm)
2 rem PARM(2) = exp for (m/mo)
3 rem PARM(3) = exp for H+
4 rem PARM(4) = exp for O2
5 rem PARM(5) = exp for Cl-
10 if (m <= 0) then goto 200
20 if (SI("Galena") >= 0) then goto 200
30 lograte = -11.8 + PARM(1) + PARM(2)*log10(m/m0) + PARM(3)*lm("H+") + PARM(4)*lm("O2") + PARM(5)*lm("Cl-")
40 moles = (10^lograte)*TIME
50 if (moles > m) then moles = m
200 SAVE moles
-end

KINETICS 1

Galena
   -formula PbS
      -m 0.00418                  #1 g PbS(mol/239.266g) = 4.18E-3 mol PbS
      -m0 0.00418
      -PARMS 0.235 0.66 0.41 0.41 0.26   
                        ##0.235 = from BET
                     ##second parm is assumed to be 2/3 for uniformly dissolving spheres and cubes (Parkhurst & Appelo,1999)
                     ##0.41 is the reaction order of H+ from own calculation
                      ## 0.41 is the reaction order of O2 from own calculation
                     ## 0.26 is the reaction order of Cl- from own calculation
-steps       25200 in 14 steps # 7 hours                     
            #432000 is 5 days
                        #25200 is 7 hours
-step_divide 1
-runge_kutta 6
-bad_step_max 100

INCREMENTAL_REACTIONS True

SOLUTION 1 Synthetic seawater
     units   ppm
        pH      4
        density 1.023
        temp    20.0
        redox   O(0)/O(-2)
        Ca              412.3
        Mg              1291.8
        Na              10768.0
        K               399.1
        Cl              19353
        Pb              1 ppb
        S(6)            2712.0
        O(0)            3.5   
   -water 1.03

#EQUILIBRIUM_PHASES 1
     #Anglesite 0 0
     #Cerussite 0 0

USER_GRAPH 1
    -headings               time Pb pH O
    -axis_titles            "Time (Seconds)" "Pb (ppb) " "pH"
    -initial_solutions      false
    -axis_scale x_axis      auto auto auto auto
    -axis_scale sy_axis     1 7  auto auto
    -connect_simulations    true
    -plot_concentration_vs  time
  -start
10 GRAPH_X total_time
20 GRAPH_Y TOT("Pb")*207.2*1000*1000  #, color = Orange, symbol = Square, symbol_size = 6, y-axis = 1
30 GRAPH_SY -LA("H+") #, color = Blue, symbol = Square, symbol_size = 6, y-axis = 1
40 GRAPH_SY TOT("O")/16 #, color = Orange, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end

USER_GRAPH 2
    -headings               time Anglesite (SI)
    -axis_titles            "Time (Seconds)" "Anglesite (SI)"
    -initial_solutions      false
    -axis_scale x_axis      auto auto auto auto
   # -axis_scale y_axis      0 50 5 auto
   # -axis_scale sy_axis     1 7  auto auto
    -connect_simulations    true
    -plot_concentration_vs  time
  -start
10 GRAPH_X total_time
20 GRAPH_Y SI("Anglesite")*207.2*1000*1000  #, color = Orange, symbol = Square, symbol_size = 6, y-axis = 1
1000 REM end
  -end

[dlparkhurst]

Are the units of your equation mol/m^2^s? For pyrite (Willemson and Rimstidt) the log rate K(pyrite) is -8.19. Yours is -11.18, which is three orders of magnitude less. If you multiply your rate by 1000, then you are close to the Pb increase in your figure (which is mol/kgw), so consider why your log K is so small.

Your surface area would be just 0.2 m^2 (for 1 g sample) which you would multiply the rate. To convert the rate to moles, you need to multiply by surface area and time. Note parm(1) is log10, so log10(.2) = -0.7 the parm(1) value. With this adjustment, your rate is about 4 orders of magnitude too slow.

So my guess is your K is wrong.

And the SI is log10(iap/K), which requires no unit conversion, so you should simply plot SI("Anglesite") in USER_GRAPH 2.

[rfembilejr]

Yes, the unit is in mol/m^2 s..This rate was for reaction in seawater that is why it is relatively lower. Similar studies for chalcopyrite, pyrrhotite and pyrite have k values of -9.38, -7.27 an -11.02 respectively. I will have to play with the equation and parameters.

Thanks!

[rfembilejr]

Just some clarification, for pyrite, accdg to Williamson and Rimstidt (1994), shouldn't it -10.19 for the K?

[dlparkhurst]

Check it out. Geochimica et Cosmochimica Acta
Volume 58, Issue 24, December 1994, Pages 5443-5454

[rfembilejr]

My bad. I was mainly pointing to the -10.19 expressed in the example for RATES in the phreeqc manual.
I have encountered this number for pyrite reaction rates in several phreeqc examples. Is it valid, at all?

[dlparkhurst]

As long as all of the units work out to moles.

[rfembilejr]

Copy. Thanks.

[Scarlet]

Hello,

if you allow, I want to know, is that I can use a negative pH on phreeqc because when I perform a galena analysis with hcl, it does not react.
thanks also to see if I made a mistake:

RATES
    Galena
-start
 10 REM parm(1)=log10(A/V,1/dm)
 20 REM parm(2)=exp for (m/m0)
 30 REM parm(3)=exp for H+
 40 REM parm(3)=exp for O2
 50 REM parm(3)=exp for Cl-
 60 IF(M <= 0) THEN GOTO 110
 70 IF(SI("Galena") >= 0) THEN GOTO 110
 80 lograte = -11.8+PARM(1)+PARM(2)*LOG10(M/M0)+PARM(3)*LM("H+")+PARM(4)*LM("O2")+PARM(4)*LM("Cl-")
 90 moles = (10^lograte)*TIME
100 IF(moles > M) THEN moles = M
110 SAVE moles
-end


KINETICS 1
Galena
    -formula  PbS  1
    -m        0.00041
    -m0       0.00041
    -parms    0.235 0.66 0.41 0.41 0.26
    -tol      1e-08
-steps       1 in 14 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500

INCREMENTAL_REACTIONS True


SOLUTION 2
    temp      20
    pH        0
    pe        4
    redox     pe
    units     ppm
    density   1
    Cl        80000
    Pb        1.2 ppb
    S(6)      2712
    -water    1 # kg
EQUILIBRIUM_PHASES 1
    Cerussite 0 0


USER_GRAPH 1
    -headings               time Pb pH O
    -axis_titles            "Time" "Pb" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y TOT("Pb")*207.2*1000*1000
  -end
    -active                 true

[dlparkhurst]

Reconsider your rate formulation, or your solution composition.

Lograte has a term PARM(4)*LM("O2"), but LM("O2") is on the order of -70, hence, the log rate is a large negative number and the rate is a very small number.

Your solution has no oxygen and O2(g) is not included in the equilibrium phases, so there is no oxidant available.

[rfembilejr]

Hi Scarlet,

In addition, there is a disconnect on the assigning of parameters. Check your parameters.
Parm(3) is repeated thrice. I believe the most recent one will overwrite the previous.
I also notice that you are using a logK of -11.8. Note that this was for galena dissolution in seawater. It may not be appropriate for a very acidic solution.

Cheers,
Rodrigo

[Scarlet]

Okay thank you very much.
if you can help me find the equation of galena that I have to put in rate.
knowing that I seek the cinqueque of galena by adding Hcl.in a lake. in ph acide

[Scarlet]

AND WHEN I wanted to train by using your example Mr RODRIGO. I did not have the same graph as you. knowing that I use phreeqc 3.3.12 (Phreeqc.dat)