The factor influencing pe in PhreeqcRM

[inhee]

Hello. I'm perfoming reactive transport models using COMSOL and PhreeqcRM.

I got a problem about convergence in my model.

I simulated injecting oxidized water into reduced groundwater and wanted to identify uranium dissolution.

However, those error messages appeared in PhreeqcRM over the time step. Finally, failed to converge and got an error.

Code: [Select]

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying reduced tolerance 1e-14 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying reduced tolerance 1e-15 ...

ERROR:                    C has not converged. Total: 5.607351e-04 Calculated: 5.607253e-04 Residual: 9.804688e-09

ERROR:                    S has not converged. Total: 8.643901e-04 Calculated: 8.644072e-04 Residual: -1.709325e-08

ERROR:                    U has not converged. Total: 6.704384e-09 Calculated: 3.089464e-09 Residual: 3.614920e-09

ERROR:                   Mu Ionic strength has not converged. Residual: 1.037241e-08

ERROR:                   pH Charge balance has not converged. Residual: -9.487852e-09

ERROR: Numerical method failed on all combinations of convergence parameters, cell/soln/mix 966

ERROR: PhreeqcRM failed.
ERROR: PhreeqcRM::RunCells


Failure in PhreeqcRM

I tried considering the transport of H, O, H2O, Charge and so on...

But I'm not sure that Charge affects pe of the water in PhreeqcRM. I checked charge was varied when I changed different redox condition of water, though.

In the flow model, Charge was assumed as one of species. It is transported and imported to PhreeqcRM by cells. The value of Charge is the range of -2.8e-5 (reducing) ~ 0.0002 (oxidizing). Is it an appropriate method?

In my model, pe was so fluctuated (3~13, not direct influence of oxidized water) (sometimes pH, too). It may cause convergence error.

What makes the convergence error? Is there any solution to help convergence in Phreeqc? (I tried to increase convergence_tolerance, tolerance, step_size, and pe_step_size, but that can't solve this problem I guess)

Here is the Phreeqc code

Code: [Select]

SOLUTION 1
units mg/l
temp 18.9
pH 6.95
pe 3.6
Alkalinity 97.8 as HCO3
        Ca      36.2
K 0.98
Mg 5.96
Na 8.26
Cl 2.92
Fe 0.067825
S 43.4
U 0.000356
        Mn      0.34512

EQUILIBRIUM_PHASES 1
    Fe(OH)3(a) 0.4

RATES 1
   Fe
    -start
    10 k = 2e13 #2e13
    20 r = (2.91e-9 + k * (MOL("OH-"))^2 * TOT("Do")/2 / 0.0015236 )
    30 if (r > 0.3) then r = 0.3
    40 rate = r * MOL("Fe+2")
    50 moles = rate  * (60* 0.05)
    200 save    moles
    -end
  

[inhee]

I found that Ferrous iron (Fe(2)) oxidation reaction is too fast. So, Fe concentration became negative. It made errors.

I don't know what is the order of caclulation about EQUILIBRIUM_PHASES and KINETICS in Phreeqc. The kinetic rate is too fast (more than 1 Fe(2)).

When I do not consider ferrous oxidation reaction, the error disappeared.

[dlparkhurst]

You have not given the KINETICS definition, so I don't know what is reacting, but the RATES definition has problems.

TOT("Do") is not defined, you probably want to use MOL("O2") to get the molality of O2(aq).

You need to multiply the rate time the internal variable TIME to calculate the moles of the reaction. The integration method uses subtime steps of variable length, which is given by TIME, so you cannot use a fixed time interval (0.05*60) in the RATES definition.

You should test your reaction network with a PHREEQC batch reaction or TRANSPORT calculation to make sure the chemistry is doing what you think it should before you do the PhreeqcRM reactive-transport modeling. I am limited in helping you with your model, but I can help debug a PHREEQC calculation.

Code: [Select]

RATES 1
   Fe
    -start
    10 k = 2e13 #2e13
    20 r = (2.91e-9 + k * (MOL("OH-"))^2 * TOT("Do")/2 / 0.0015236 )
    30 if (r > 0.3) then r = 0.3
    40 rate = r * MOL("Fe+2")
    50 moles = rate  * (60* 0.05)
    200 save    moles
    -end


[inhee]

Thank you for replying, Parkhurst.

I set dissolved oxygen concentration as Do (I defined, but I didn't write upper code) because oxygen made pe too high than I expected.

In addition, although I don't know exact mechanism, rate is calculated 20 times at time step 60 minutes in PhreeqcRM. So, I set RATES like that.

Anyway, I'm curious how pe is determined in PhreeqcRM. Is the way to transport "Charge (a component in PhreeqcRM)" like other species concentration all right?

In my model, pe is not stabilized. It causes the convergence problem and abnormal surge and plunge Fe concentration by EQUILIBRIUM_PHASES Fe(OH)3(a).

[dlparkhurst]

You must use the TIME variable in your rates definition. Do not argue.

PhreeqcRM is designed to have either H, O, and charge transported, or H2O, excess H, excess O, and charge transported. Do not argue. If you charge balance all of your solutions, you should have near-zero charge (<1e-8) that you are transporting in almost all cases.

PHREEQC solves all of the algebraic equations simultaneously, including H, O, charge, and all other mole balance equations and mass-action equations. pe is one of the variables that is solved.

As you said, you calculation is successful if you remove KINETICS; therefore, I assume you RATES or KINETICS definition is the problem. I will help you debug your rates definitions if you make a stand-alone PHREEQC input file and post it. I probably will not be able to help you with problems in PhreeqcRM in your COMSOL application.

[dlparkhurst]

I think you are placing too much importance on what you think the pe should be. If there is even a small amount of DO (say, 1e-8 mol/kgw), the pe will be high.

[inhee]

I appreciate your answer. I wrote the all Phreeqc code and selected output of PhreeqcRM because I can't attach the file in this post.

Code: [Select]

TITLE GS

SOLUTION 1 # groundwater
units mg/l
temp 18.9
pH 6.95
pe 3.6
Do     1.09
Alkalinity 97.8 as HCO3
    Ca      36.2
K 0.98
Mg 5.96
Na 8.26
Cl 2.92
Fe 0.067825
S 43.4
U 0.000356


SOLUTION 2 # surface water
units mg/l
temp 24.5
pH 7.12
pe 10.8
Do     10.66
Alkalinity 29.2 as HCO3
  Ca      29
K 0.77
Mg 9.67
Na 3.91
Cl 3.26
Fe 0.012247
S 83.3
U 0.001602

SOLUTION_MASTER_SPECIES
Do Do-2 0.0 Do 16   

SOLUTION_SPECIES
Do-2 = Do-2
      log_k  0.0

MIX 1
1 0.5
2 0.5
SAVE solution 3
END


EQUILIBRIUM_PHASE 1
Fe(OH)3(a) 0.8

RATES 1
Uraninite
-start
10 k1 = 0.00273
20 k2 = 0.00152
30 k_1 = 5e-5
40 rate =  (k1 * k2 * 22e-6 / 238 * TOT("Do") * MOL("HCO3-")) / (k_1 + k2 * MOL("HCO3-") + k1 * TOT("Do"))
50 moles = rate * TIME * 2.3e-3 #Surface area of UO2
200 save moles
-end

KINETICS 1
Uraninite
-formula U 1
-steps 255600
USE solution 3

ENDMy PhreeqcRM model simulates the dilution of groundwater and surface water like mixing in Phreeqc.

In my model (even though in Phreeqc), charge (maybe Electrical balance in Phreeqc?) is more than 1e-4 in absolute value. Isn't the solution appropriate in terms of charge balance? I don't know the reason.

In MIXING of Phreeqc, pH and pe seems like normal, and Fe concentration makes sense. However, in PhreeqcRM, pH, pe, and Fe concentration is abnormal (negative pe, too high Fe) like the Selected output.

What is the cause? Can diffusion and dispersion in a transport model disrupt charge balance?

#Selected output of PhreeqcRM

Code: [Select]

Cell number 4674
     Selected output:
          1 soln: 4673
          2 pH: 9.4081
          3 pe: -5.3467
          4 Alk(eq/kgw): 0.00081813
          5 H(mol/kgw): 1.1819e-11
          6 U(mol/kgw): 9.4951e-08
          7 Fe(mol/kgw): 0.00013725
          8 Do(mol/kgw): 0.00052216
          9 m_OH-(mol/kgw): 1.2373e-05
          10 m_HCO3-(mol/kgw): 0.00048401
          11 m_Fe+2(mol/kgw): 6.2171e-05
          12 m_Fe(OH)3(mol/kgw): 1.2581e-08

[dlparkhurst]

I don't understand what you are trying to do with Do. I think you want to use O(0) for dissolved oxygen in the SOLUTION definition. The gram-formula weight assigned to O(0) is 16, so in a SOLUTION definition in units of mg/L the mg/L will be divided by 16 to get moles, and the moles of O2(aq) will be 1/2 of the moles of O(0).

I interpret your script that you are oxidizing uranium metal, because you have -formula U 1. I have increased your rate by a factor of 1e9 to show the consumption of U metal and dissolved oxygen, and U(6) is produced. Perhaps this is not what you have in mind.

Code: [Select]

TITLE GS
RATES 1
Uraninite
-start
10 k1 = 0.00273
20 k2 = 0.00152
30 k_1 = 5e-5
#40 rate =  (k1 * k2 * 22e-6 / 238 * TOT("Do") * MOL("HCO3-")) / (k_1 + k2 * MOL("HCO3-") + k1 * TOT("Do"))
40 rate =  (k1 * k2 * 22e-6 / 238 * TOT("O(0)") * MOL("HCO3-")) / (k_1 + k2 * MOL("HCO3-") + k1 * TOT("O(0)"))
50 moles = rate * TIME * 2.3e-3 #Surface area of UO2
      60    moles = moles * 1e9
200 save moles
-end
END
SOLUTION 1 # groundwater
units mg/l
temp 18.9
pH 6.95
pe 3.6
#Do     1.09
O(0)      1.09
Alkalinity 97.8 as HCO3
    Ca      36.2
K 0.98
Mg 5.96
Na 8.26
Cl 2.92
Fe 0.067825
S 43.4
U 0.000356


SOLUTION 2 # surface water
units mg/l
temp 24.5
pH 7.12
pe 10.8
#Do     10.66
O(0)      10.66
Alkalinity 29.2 as HCO3
  Ca      29
K 0.77
Mg 9.67
Na 3.91
Cl 3.26
Fe 0.012247
S 83.3
U 0.001602

MIX 1
1 0.5
2 0.5
SAVE solution 3
END
USE solution 3

EQUILIBRIUM_PHASE 1
Fe(OH)3(a) 0.8

KINETICS 1
Uraninite
-formula U 1
-steps 255600 in 10 steps

USER_GRAPH 1
    -headings               time U O(0)
    -axis_titles            "Time, seconds" "Molality" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y TOT("U"), TOT("O(0)")
  -end
    -active                 true
END