PHREEQC versus theoretical precipitation via ICE tables

[freiburgermsu]

Dear PHREEQCusers,

I am verifying the precipitation calculations of PHREEQC for my Supporting Information. I created a theoretical ICE table of the attached images based upon a model system of precipitating Gypsum (https://github.com/freiburgermsu/ROSSpy/tree/main/examples/scaling/2021-10-27-ROSSpy-red_sea-transport-pitzer-scaling-all_distance-LinPerm/input.pqi). The predicted precipitation of Gypsum, based upon the initial elemental concentrations is 0.020285 molal (0.352 moles), however, the PHREEQC prediction is 0.823 moles over all cells of the transport column(https://github.com/freiburgermsu/ROSSpy/tree/main/examples/scaling/2021-10-27-ROSSpy-red_sea-transport-pitzer-scaling-all_distance-LinPerm/selected_output.pqo). Did I incorrectly interpret the PHREEQC results, or must something else be considered?

I appreciate any ideas or suggestions.

Thank you :)
  Andrew

[dlparkhurst]

OK, I get 0.822 mol of gypsum in your column at the end of the simulation.

I don't know where you get 0.020285 molal (0.352 moles) for the "predicted precipitation of gypsum" ... "based on the elemental concentrations". Moles is probably the better unit if you you are removing water. Note the mass of water in the TRANSPORT calculation varies from 17.38 to 15.59, and you are running more than 4 pore volumes through the column, with a different H2O extraction in each cell.

The volume of water reaching each cell precipitates a given amount of gypsum depending on the water extraction and the water extraction in the previous cells. The infilling solution reaches cell 1 51 times, cell 2 50 times, ... cell 12 40 times.

[freiburgermsu]

Dear David,

I forgot to attach the worked steps in my previous post; they are now attached.

I determined the expected precipitation after the entire column of simulated cells, based upon the initial concentrations. Is this the wrong approach?

I am eager for your perspective of the work.

Thank you :)
    Andrew

[dlparkhurst]

The log K is the same as in pitzer.dat. The equilibrium relation is

Code: [Select]

log K = LA(Ca+2) + LA(SO4-2) + 2*LA(H2O)

where LA represents log10 activity. Your approximation that total concentrations are equal to activity is not going to be consistent with a PHREEQC calculation. A PHREEQC calculation for the conditions you cite is given below. 

Code: [Select]

SOLUTION
-units mol/kgw
S(6) 0.1
Ca   0.02
EQUILIBRIUM_PHASES
Gypsum 0 0
END

The result with pitzer.dat is precipitation of 0.012 mol of gypsum, compared to your calculation of 0.02 mol. The difference is in the activity coefficients, activity = gamma*molality. Further, this simplified calculation ignores the effects of other specific-ion interactions that are considered when you add the complete solution composition. When your original, complete solution composition is used, PHREEQC calculates that no gypsum precipitates at all.

In addition, regardless of the details, the ICE calculation is not the same as the TRANSPORT calculation. In the transport calculation, a solution with 17.38 kgw enters cell 1; 9.06 moles of water is removed; and gypsum is allowed to precipitate if supersaturated (but does not). At step 2, (ignoring dispersion) this solution is advected to cell 2, 8.95 moles of water is removed and still no gypsum precipitates. During the same time step, another solution with 17.38 kgw enters cell 1, another 9.06 moles of water is removed from cell 1. And so on for 51 steps. Considering that PHREEQC does not calculate any gypsum precipitation until the sixth cell, your ICE calculation is not relevant to the transport calculation.

In a purely advective system, once solution 0 has reached a cell, the amount of gypsum precipitated in a cell is the same for each subsequent step. Solution 0 reaches cell 1--51 times, cell 2--50 times, ..., and cell 12--40 times. So, the sum of delta gypsum at time step 51 for each cell multiplied by the corresponding number of steps will give the total for the column. The simulation below produces the following output, where the last two column are the same, but one multiplies the delta gypsum by the appropriate number of shifts, and the other is the total gypsum accumulated in each cell:

Code: [Select]

    delta_gyp       shifts delta*shifts    equi(gyp)
  0.0000e+00   5.1000e+01   0.0000e+00   0.0000e+00
  0.0000e+00   5.0000e+01   0.0000e+00   0.0000e+00
  0.0000e+00   4.9000e+01   0.0000e+00   0.0000e+00
  0.0000e+00   4.8000e+01   0.0000e+00   0.0000e+00
  0.0000e+00   4.7000e+01   0.0000e+00   0.0000e+00
  2.4088e-04   4.6000e+01   1.1080e-02   1.1080e-02
  3.3186e-03   4.5000e+01   1.4934e-01   1.4934e-01
  3.2829e-03   4.4000e+01   1.4445e-01   1.4445e-01
  3.2469e-03   4.3000e+01   1.3962e-01   1.3962e-01
  3.2109e-03   4.2000e+01   1.3486e-01   1.3486e-01
  3.1747e-03   4.1000e+01   1.3016e-01   1.3016e-01
  3.1384e-03   4.0000e+01   1.2554e-01   1.2554e-01
Column sum:                8.3504e-01   8.3504e-01

Here is the script.

Code: [Select]


SOLUTION 0 red_sea
temp 24.5 #average of Al-Taani et al., 2014 and Longinelli and Craig, 1967..
pH 8.22 charge #None
pe 0.2679    #Al-Taani et al., 2014 // 4.00 is the default (?) // 4.00 is the default (?)
units ppm
Mn 0.000306 #Al-Taani et al., 2014 for the Northern Gulf of Aqaba
Fe 0.006281 #Al-Taani et al., 2014 for the Northern Gulf of Aqaba
B 1.344 #Al-Taani et al., 2014
Cl 24756 #https://www.lenntech.com/composition-seawater.htm in the Red Sea, and Longinelli and Craig, 1967
Na 16417.2 #https://www.lenntech.com/composition-seawater.htm in the Red Sea, and Longinelli and Craig, 1967 describes [Na]=15834
S(6) 9500 #Longinelli and Craig, 1967 and Llyod, 1967
Ca 774 #Abdel-Aal et al., 2015
K 301 #Abdel-Aal et al., 2015
Mg 1646 #Abdel-Aal et al., 2015
Sr 8.3 #Bernat, Church, and Allegre, 1972 from the Mediterranean
Ba 0.011 #Bernat, Church, and Allegre, 1972 from the Mediterranean
Li 0.228 #Stoffyn-Egli and Mackenzie, 1984 for the Mediterranean Sea
-water 17.378153556381264
END

SOLUTION 1-12 Initial solution in the RO module
temp 24.5
units ppm
Mn 0
Fe 0
B 0
Cl 0
Na 0
S(6) 0
Ca 0
K 0
Mg 0
Sr 0
Ba 0
Li 0
-water 17.378153556381264
END
EQUILIBRIUM_PHASES 1-12
gypsum 0 0
END

REACTION 1
H2O -1; 9.057149147747937
REACTION 2
H2O -1; 8.952016761697115
REACTION 3
H2O -1; 8.846884375646292
REACTION 4
H2O -1; 8.74175198959547
REACTION 5
H2O -1; 8.636619603544647
REACTION 6
H2O -1; 8.531487217493824
REACTION 7
H2O -1; 8.426354831443
REACTION 8
H2O -1; 8.321222445392177
REACTION 9
H2O -1; 8.216090059341354
REACTION 10
H2O -1; 8.110957673290532
REACTION 11
H2O -1; 8.005825287239709
REACTION 12
H2O -1; 7.9006929011888865
END
SELECTED_OUTPUT 2
-file RedSea.sel
USER_PUNCH 2
-headings delta_gyp shifts  delta*shifts equi(gyp)
10 IF TOTAL_TIME <= 0 THEN GOTO 200
20 d = EQUI_DELTA("gypsum")
30 n = 52 - CELL_NO
40 sum_gyp = d * n
50 PUT(GET(1) + sum_gyp, 1)
55 PUT(GET(2) + EQUI("gypsum"), 2)
60 PUNCH d, n, sum_gyp, EQUI("gypsum")
70 IF CELL_NO < 12 THEN GOTO 300
80 PUNCH EOL$
90 PUNCH "Column sum:             ", GET(1), GET(2)
110 GOTO 300
200 PUNCH NO_NEWLINE$
300 END

TRANSPORT
-cells 12
-shifts 51
-lengths 0.08466666666666667
-time_step 3.946329703005324 # this satisfies the Courant condition with a feed velocity of 2.575E-1 m/s
-initial_time 0
-boundary_conditions constant flux # Dirichlet boundary condition
-diff 0
-disp 0
-punch_cells 1-12
-punch_frequency 51
-print_frequency 51
END


[freiburgermsu]

Dear David,

I did not realize that activity coefficients would deviate significantly far from 1 for divalent ions in a concentrated geochemical solution. The significance of other ionic interactions is curious as well. These myriad influences upon the geochemical system would certainly explain deviation from the theoretical and over-simplified ICE table.

Why do the values in the "equi(gyp)" column deviate slightly from the "Gypsum" column of the SELECTED_OUTPUT file, since both column sum precipitation in a specific cell over all timesteps?

Thank you for your support :)
  Andrew

[dlparkhurst]

I suspect you are comparing two different simulations with slightly different conditions. In the latest simulation, I reset the temperature of the initial solutions (1-12) to 24.5 C to make the simulation isothermal.

The biggest activity coefficient effect is simply the concentration of Ca+2 and SO4-2. Even without addition of NaCl, the mean activity coefficient is about 0.2 compared to 1.0 in your ICE calculation. Addition of NaCl further reduces the mean activity coefficient and activity of water such that the saturation index of gypsum switches from positive to negative.

Code: [Select]

SOLUTION
-units mol/kgw
S(6) 0.1
Ca   0.02
REACTION
NaCl 1
1 in 20 steps
USER_GRAPH 1
    -headings               rxn Mean_Gamma LA(H2O) SI(Gypsum)
    -axis_titles            "NaCl added, moles" "Gamma, LA(H2O), or SI"
    -axis_scale x_axis      auto 1 auto auto
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X RXN
20 GRAPH_Y SQRT(GAMMA("Ca+2")*GAMMA("SO4-2")), LA("H2O")
30 GRAPH_Y SI("Gypsum")
  -end
END


[freiburgermsu]

Dear David,

I was inspired by your post. I added this USER_PRINT block to the PQI file

Code: [Select]

USER_PRINT
10 PRINT "gamma_Ca: ", GAMMA("Ca+2")
20 PRINT "gamma_SO4: ", GAMMA("SO4-2")
and I discovered that gamma_Ca =0.19 and gamma_SO4 = 0.06 throughout the simulation. I updated the ICE table calculations and acquired 0.18 moles of Gypsum precipitation, which is much closer to the 0.21 moles that is predicted from the pure system of only Ca and SO4. Perhaps the remaining ~17% of error is attributable to the ionic interactions.

I didn't notice the temperature adjustment. Indeed, perhaps the slightly lower temperature in your simulation and the high heat capacity of water encouraged more precipitation than in the original simulation.

Thank you :)
    Andrew