Relationship between uraninite and gypsum

[Jeonghwan Hwang]

I modeled a situation where a high concentration of uranium was leaked into bentonite in which various minerals were present. Because it is a reducing environment, uranium was set to precipitate as uraninite, and the model was simulated.
As a result of the model, it was confirmed that uraninite was mostly precipitated. But here all the gypsum disappeared, and pyrite was deposited. The pH and eh were 10 and -11, respectively, resulting in slightly strange calculations.
It is clear that the change in 'S' due to the dissolution of the gypsum caused the precipitation of pyrite, which appears to increase the pH.
However, I could not have an idea why the gypsum rapidly dissolved as precipitating the uraninite.  I searched the LLNL database to determine why the uraninite precipitation affects the abnormal reaction of gypsum and pyrite, but I could not clearly understand.
Can I get the advise?
My code is written as below;

=====================================================
Database C:\PHREEQC\database\llnl.dat

PHASES # Sit TDB
Calcite
CaCO3      = 1.000Ca+2     + 1.000CO3-2
     log_k    -8.480     #82PLUM/BUS
     delta_h  -10.620    #kJ/mol        #82PLUM/BUS
     # Enthalpy of formation:           -1207.61        #kJ/mol       
     -analytic -1.03405E+1 0E+0 5.5472E+2 0E+0 0E+0
     
     Dolomite
CaMg(CO3)2      = 1.000Ca+2     + 1.000Mg+2     + 2.000CO3-2
     log_k   -17.120     
     delta_h  -35.960    #kJ/mol       
     # Enthalpy of formation:           -2324.5        #kJ/mol        #95ROB/HEM
     -analytic -2.34199E+1 0E+0 1.87832E+3 0E+0 0E+0
end

EXCHANGE_MASTER_SPECIES
   Z   Z-
   
EXCHANGE_SPECIES #SKB TR-06-16 Input data
   Z- = Z-
      log_k   0.0
   Z- + Na+ = NaZ
      log_K 0.0
   Z- + K+ = KZ
      log_k 0.6
   2Z- + Ca+2 = CaZ2
      log_k 0.41
   2Z- + Mg+2 = MgZ2
      log_k 0.34   
   
SURFACE_MASTER_SPECIES
   Mont_wa   Mont_waOH
   Mont_wb   Mont_wbOH

SURFACE_SPECIES
   Mont_waOH = Mont_waOH
   log_k 0
   Mont_waOH + H+ = Mont_waOH2+
   log_k 4.5
   Mont_waOH = Mont_waO- + H+
   log_k -7.9
   
   Mont_wbOH = Mont_wbOH
   log_k 0
   Mont_wbOH + H+ = Mont_wbOH2+
   log_k 6.0
   Mont_wbOH = Mont_wbO- + H+
   log_k -10.5
   
SOLUTION 1#solnum      
   units mol/L   
   pH   7.2
   pe   -2.42
   Temp   27
   C   2.20E-03 as HCO3
   Ca   2.33E-02
   Cl   1.53E-01
   Fe   3.31E-05
   K   8.75E-04
   Mg   9.30E-03
   Na   8.88E-02
   S   6.80E-03 as SO4
   Si   1.85E-04
   Br      1E-15
   U       0.27
EXCHANGE 1#solnum   
      NaZ   0.8465
      CaZ2   0.1025
      KZ   0.0235
      MgZ2   0.0470   
SURFACE 1 #solnum   
   -sites_units absolute
   Mont_waOH   0.0627   
   Mont_wbOH   0.0627   
   no_edl   
EQUILIBRIUM_PHASES 1#solnum   
   Quartz  0   1.3063
   Pyrite   0.0   0.0092
   Gypsum   0   0.0807
   Calcite  0   0
   Dolomite  0  0
   Siderite  0  0
   Uraninite 0  0
END

selected_output
-file Hwang.txt
USER_PUNCH
-headings Ca Na Mg K S Fe cal dol gyp qrt sid pyr Ura
10 punch TOTMOLE("Ca") TOTMOLE("Na") TOTMOLE("Mg") TOTMOLE("K") TOTMOLE("S") TOTMOLE("Fe") EQUI_DELTA("Calcite") EQUI_DELTA("Dolomite")  EQUI_DELTA("Gypsum") EQUI_DELTA("Quartz") EQUI_DELTA("Siderite") EQUI_DELTA("Pyrite") EQUI_DELTA("Uraninite")
-end
end

Run_cells
-cells 1
end

Solution_Modify 1
-totals
Br   6.33E-09
C   0.001982682
Ca   0.017207166
Cl   0.154461215
Fe(2)   0.000333815
Fe(3)   4.19E-11
K   0.0015062
Mg   0.006799786
Na   0.187869226
S(-2)   0.001460669
S(6)   0.046997382
Si   5.83E-05
U   0.272357417
EXCHANGE_Modify 1
SURFACE_Modify 1
EQUILIBRIUM_PHASES_Modify 1
Run_cells
-cells 1
end    

[dlparkhurst]

I think the problem is your use of SOLUTION_MODIFY. You change the concentrations of most solutes, but you do not adjust total H, total O, and charge balance. I don't know where you are getting the values for SOLUTION_MODIFY, but if you are transporting constituents, you must also transport H, O, and charge balance to define a new solution. If you are using a new solution, you should probably define a new SOLUTION 1 to be sure that all quantities (including total H, total O, and charge) are defined correctly. Bottom line, you cannot redefine a solution the way you used SOLUTION_MODIFY; you are missing the definitions of -total_h, -total_o, and -cb.

Here is a modification of your script that shows what your solution looks like after SOLUTION_MODIFY. It simply gives the equilibrium composition of the solution with no additional reactions. There is almost 1 mol/kgw of H(0) [0.5 mol/kgw H2(aq)], which provides a large reduction capacity. The large reduction capacity reduces all of the sulfate in the system, including all the sulfate in gypsum, and all of the carbon, including the carbon in dolomite.

Code: [Select]

 
SOLUTION 1#solnum     
   units mol/L   
   pH   7.2
   pe   -2.42
   Temp   27
   C   2.20E-03 as HCO3
   Ca   2.33E-02
   Cl   1.53E-01
   Fe   3.31E-05
   K   8.75E-04
   Mg   9.30E-03
   Na   8.88E-02
   S   6.80E-03 as SO4
   Si   1.85E-04
   Br      1E-15
   U       0.27

EQUILIBRIUM_PHASES 1#solnum   
   Quartz  0   1.3063
   Pyrite   0.0   0.0092
   Gypsum   0   0.0807
   Calcite  0   0
   Dolomite  0  0
   Siderite  0  0
   Uraninite 0  0
END


Run_cells
-cells 1
end

Solution_Modify 1
-totals
Br   6.33E-09
C   0.001982682
Ca   0.017207166
Cl   0.154461215
Fe(2)   0.000333815
Fe(3)   4.19E-11
K   0.0015062
Mg   0.006799786
Na   0.187869226
S(-2)   0.001460669
S(6)   0.046997382
Si   5.83E-05
U   0.272357417
USE solution 1
REACTION
END


[dlparkhurst]

One other comment. It doesn't make sense to me to have 0.1 mol/kgw of U(4) in solution unless, possibly, the pH were extremely low. I think the solution would be U(6) to be mobile, and then reduction occurs through some process. In a roll-front deposit, U(6) is present in an oxidizing environment and then gets transported to a reducing environment, usually caused by organic decomposition.

[Jeonghwan Hwang]

I found the H(0) concentration in my model results, and I was able to recognize that this was the problem. Unfortunately, when I modified total_h and total_o in Solution_modify, the problem was not solved. I think that this problem appears to be because the concentration of U is too large (0.27 mol/kgw). Following your additional advice, I plan to modify the model to input the 0.27 mol of uraninite in equilibrium_phases and track the behavior of dissolved uranium concentration from uraninite.

Thanks for the advice.

[dlparkhurst]

If you use SOLUTION_MODIFY, you must account for H, O, and charge; otherwise, your results will not be reliable.

If you are using a transport model, you should use PhreeqcRM, which is designed to be coupled with multicomponent transport models.

[Jeonghwan Hwang]

I had one more question.

I can take out the amounts of H and O using User_punch as below;

Selected_output
-reset false
user_punch
-heading H O
punch 10 TOTMOLE("H") TOTMOLE("O")
-end
End

and I can input the H and O using Solution_modify -total_h and total_o.

However, I cannot have an idea to take out the information about the charge (-cb). I know that the -cb means charge balance.
In addition, I cannot understand how charge balance transports like H and O.

Can I get an advice for that?
Thank you

Sincerely,
Jeonghwan Hwang

[dlparkhurst]

First, you should use PhreeqcRM. I think it would avoid a number of problems. It automatically provides H, O, and charge, converts units, handles density, saves state, and has many other features that you do not have to program, including parallelization.

The Basic function charge_balance provides the charge imbalance in units of equivalents.

Code: [Select]

10 PUNCH TOTMOLE("H"), TOTMOLE("O"), charge_balance

Consider the charge imbalance as some unspecified ion (that can be either positive or negative). It is possible the charge balance is always zero, but if you use real analyses or some surface reactions, the charge imbalance may not be zero and may vary within the domain (your solution 1 has a charge imbalance of  -1.489e-02 eq). If you do not account for it, you will get unusual ph (your high pH result) or redox reactions (your reduction of gypsum). So, in general, it is best to transport the charge to account for any variation in charge balance.

[Jeonghwan Hwang]

Thank you for advice.
I'm agree that I should study 'PhreeqcRM' and very glad to know the importance about charge balance.
I will study the PhreeqcRM using research paper.
Thank you for help!

[Jeonghwan Hwang]

Hello, I have an updated question.

I recently modified the code to transport 'charge', 'total H', 'total O', etc. In addition, O(0) and H(0), H+ activity (for pH), and e- (for pe) were also set to be transported together. Here, pH and pe are updated after transforming the values from H+ and e- to pH and pe at Solution_modify. The types of total ions I factored into transport are shown below. The minerals are not transported and only value changes are updated by using Equilibrium_Phases_Modify.

     SELECTED_OUTPUT
     -reset false
     USER_PUNCH
     headings Br C Ca Cl Fe K Mg Na S Si O(0) H(0) O H charge H+ e- Calcite Dolomite Gypsum Pyrite Siderite
     -start
     10 punch TOT("Br")*1000, TOT("C")*1000, TOT("Ca")*1000, TOT("Cl")*1000, TOT("Fe")*1000, TOT("K")*1000, TOT("Mg")*1000, TOT("Na")*1000, TOT("S")*1000, TOT("Si")*1000, TOT("O(0)")*1000, TOT("H(0)")*1000, TOT("O")*1000, TOT("H")*1000,
     20 punch TOT("charge")*1000
     30 punch ACT("H+")*1000, ACT("e-")*1000,
     40 punch EQUI("Calcite"), EQUI("Dolomite"), EQUI("Gypsum"), EQUI("Pyrite"), EQUI("Siderite"),
     -end

Considering all of the above, it was confirmed that the abnormal reaction between gypsum and pyrite was resolved.

My question is 'are there any components that should not be reflected in transport among the components I considered to be transported above, or components that should be additionally considered?
For example, I wonder if my model considering that the transportation of H+ and e- is conceptually acceptable.

Thank you.

Jeonghwan Hwang.

[dlparkhurst]

PhreeqcRM has proved sufficient for linking PHREEQC reactions with transport. Usually, PhreeqcRM provides concentrations for the total H in the system (including in H2O), the total O, and the charge. [Alternatively, H2O and the non-H2O H, and non-H2O O (plus charge) can be provided for transport by PhreeqcRM.]

I do not think you should transport H+, e-, O(0), or H(0). Assuming no additional chemical reactions (phases, exchange, etc), if you transport these quantities, and compare them to the results of using the transport of total H, total O, and charge followed by speciation in PHREEQC; the values for these quantities provided by transport will not be the same as the speciated values. The speciated values are the correct values, the transported values are estimates at best, but could be very much in error.

[Jeonghwan Hwang]

Your advice has been very helpful.

Additionally, I compared the transport/no transport of pH and pe, and confirmed the results that transporting pH and pe did not significantly affect the simulation results. Therefore, I decided to transport only total_H, total_O and charge.

There's something I'd like to ask for advice on.

I conducted a reactive transport model considering both bentonite and rock using COMSOL (transport) and PHREEQC (geochemical simulation).

In the simulation results, it was confirmed that the pH and pe values of the rock greatly changed at the boundary between bentonite and rock (pH = 7.2 to 3.18 and pe = -2.42 to 18.4).

When I check the reason, it seems that the initial charge balance (cb=0.091717) of bentonite is very large compared to the rock (cb=-5.6757E-12). For this reason, the charge balance of rock increases because of the charge balance transport from bentonite to rock.

The PHREEQC code calculated on the rock side is as follows, and the -cb 0.051 value applied in Solution_modify is the transported value by COMSOL.
When -cb is excluded from Solution_Modify, it can be seen that the change is relatively reduced.

=================================================================
Database D:\Apro\PHREEQC_DB\Sit.dat

SELECTED_OUTPUT
-reset false
-file HWANG_1.txt
USER_PUNCH
headings Br C Ca Cl Fe K Mg Na S Si O H charge_balance H+ e- U
-start
10 Wat_Mass = 11.5999999999999996447286321199 
20 punch TOT("Br")*1000, TOT("C")*1000, TOT("Ca")*1000, TOT("Cl")*1000, TOT("Fe")*1000, TOT("K")*1000, TOT("Mg")*1000, TOT("Na")*1000, TOT("S")*1000, TOT("Si")*1000,
30 punch TOT("O")*1000, TOT("H")*1000,
40 punch charge_balance 
50 punch -LA("H+"), -LA("e-"),
80 punch TOT("U")*1000,
-end
END

SOLUTION 1
units mol/L
pH 7.2
pe -2.42
Temp  15
C 2.20E-03 as HCO3
Ca 2.33E-02
Cl 10 charge
Fe 3.31E-05
K 8.75E-04
Mg 9.30E-03
Na 8.88E-02
S 6.80E-03 as HS
Si 1.85E-04
Br 1E-7
-water 11.5999999999999996447286321199
END

Run_cells
Cells 1
end

SOLUTION_MODIFY 1
-total_o 645.421258848828870213765185326
-total_h 1287.87814051758755340415518731
-cb 0.051387432549277474447269042912
-pH 7.20000000000000106581410364015
-pe -2.42000000000000081712414612412
-totals
Br 7.71679751984970656747972808653e-07
C 0.0397307909093962963176238645246
Ca 0.296448134983714417511890815149
Cl 2.51591540390129653204098758579
Fe 0.000597786632652432827449151542254
K 0.0191920363684371826640262526098
Mg 0.12691882955335270777297296263
Na 2.37468070298041178389780725411
S 0.344608350014510844694370916841
Si 0.00334097997446734188559114286932
U 0
END

Run_cells
Cells 1
end
=================================================================

Is there any way to overcome these problems?
Any suggestions are appreciated.

Sincerely,
Jeonghwan Hwang

[dlparkhurst]

There is not enough information for me to evaluate your issue. My guess is that you are doing something wrong with transport or units.

I don't know why you are transporting a large amount of charge. Some of the solutions must not be charge balanced, so, my first suggestion is to use charge-balanced solutions throughout. With charge-balanced solutions, you should not accumulate any charge in the solutions (except possibly for some SURFACE calculations), and all transport should result in essentially zero charge for all cells. I'm not sure this will solve your problems, but it should either point to an error, or eliminate charge as the issue.

I would set up a 1D transport calculation that can be run by your simulator and by PHREEQC. In that case, you have a way to determine where any differences occur. You can either run your chemistry and gradually simplify to isolate errors, or start simply and gradually increase the complexity of the system.