Desalination PHREEQpy applications -- ROSSpy

[freiburgermsu]

Dear PHREEQCusers,

I am developing ROSSpy (https://pypi.org/project/ROSSpy/) as an application of PHREEQC/py for simulating the geochemical reactive transport of reverse osmosis desalination. I have discovered a few inconsistent results from the PHREEQC calculations, however, which I am optimistic that some of you may be able to help clarify and resolve. I provide links to files in the ROSSpy GitHub repository of simulations and files to minimize the verbosity of this post. The PHREEQ input files for each of the referenced simulations are provided in the respective simulation folders on GitHub. The inconsistent results are detailed in the following sections. 



1) Mass balance
I created an ICE table (https://www.chem.purdue.edu/gchelp/howtosolveit/Equilibrium/Solubility_Products.htm) of Gypsum precipitation from a simple example of ROSSpy (https://github.com/freiburgermsu/ROSSpy/blob/main/examples/other/ICE_table/scaling_ice_table.ipynb). The initial and final elemental concentrations, after converting into moles, are provided in the following JSON format. The values are peculiar for deviating -- where the 1:1 stoichiometry of calcium and sulphate should create identical changes -- and for a positive change in sulfate moles, which suggests that sulfate was somehow added during the transport system.

{'Ca': {'change': -0.0006491063068378433,
        'final_moles': 0.335957721232,
        'initial_moles': 0.33660682753883786},
 'S': {'change': 0.07327891623898242,
       'final_moles': 1.7969614040000002,
       'initial_moles': 1.7236824877610177}}

Is this inconsistent result the consequence of an error of the corresponding PHREEQ input file (https://github.com/freiburgermsu/ROSSpy/blob/main/examples/scaling/2021-10-27-ROSSpy-red_sea-transport-pitzer-scaling-all_distance-LinPerm/input.pqi), or are the PHREEQ calculations being used beyond its intentions?


2) Mineral moles
The mineral and d_mineral -- e.g. Gypsum and d_Gypsum -- columns are inconsistent with each other. The d_mineral column, according to the PHREEQC version 3 documentation (https://pubs.usgs.gov/tm/06/a43/pdf/tm6-A43.pdf), represents the change in moles of the corresponding mineral during a timestep, however, the actual change in moles between timesteps in the mineral column are not equal to the d_mineral column. The Gypsum and d_Gypsum columns, for example, do not align in this simulation output (https://github.com/freiburgermsu/ROSSpy/blob/main/examples/scaling/2021-10-27-ROSSpy-red_sea-transport-pitzer-scaling-all_distance-LinPerm/selected_output.pqo).


3) Unexpected mineral precipitation trend
Mineral precipitation via most databases, including the Pitzer database, perplexingly decreases while elemental concentrations increase in a transport process. The Gypsum precipitation of the following figure evinces this phenomenon, where the precipitation quantity decreases toward the end of the transport processes despite that the elemental concentrations increase (https://github.com/freiburgermsu/ROSSpy/blob/main/examples/scaling/2021-10-27-ROSSpy-red_sea-transport-pitzer-scaling-all_distance-LinPerm/all_minerals.svg). Here is another example of this phenomenon, where halite decreases over the entire reactive transport process instead of increasing with the elemental concentrations (https://github.com/freiburgermsu/ROSSpy/blob/main/examples/sensitivity_analyses/database_selection/2021-10-27-ROSSpy-red_sea-transport-llnl-scaling-all_distance-LinPerm/all_minerals.svg). What causes this discrepancy between the ionic strength of the solution and mineral precipitation? and how can it be resolved?


4) Stagnant layer
We observe unexpected behavior when implementing a stagnant layer [https://github.com/freiburgermsu/ROSSpy/tree/main/examples/sensitivity_analyses/single_vs_dual_domain]. The mineral precipitation from a single domain (https://github.com/freiburgermsu/ROSSpy/blob/main/examples/sensitivity_analyses/single_vs_dual_domain/2021-10-27-ROSSpy-red_sea-transport-pitzer-scaling-all_distance-LinPerm/all_minerals.svg) is orders of magnitude greater than that of the stagnant layer of the dual domain (https://github.com/freiburgermsu/ROSSpy/blob/main/examples/sensitivity_analyses/single_vs_dual_domain/2021-10-27-ROSSpy-red_sea-transport-pitzer-scaling-all_distance-LinPerm-2/all_minerals.svg), despite that the stagnant layer should possess much larger concentration.

We have a related uncertainty in how the value of the Exchange Factor (EF; 1/sec) parameter influences the mobile and immobile phases of a dual domain. Comparisons of the mobile and immobile phase after parameterizing EF=1e-10 vs EF=1e10 for otherwise the same simulation (https://github.com/freiburgermsu/ROSSpy/tree/main/examples/sensitivity_analyses/exchange_factor) suggests that the value of EF is inversely related to the solvent exchange between the mobile and immobile phases. Is this the proper interpretation of the EF parameter?



I am optimistic that members of this community can clarify whether these inconsistencies originate from errors of the input files or instead manifest from limitations of the PHREEQC calculations. I appreciate any and all suggestions towards their resolution. I can arrange a video call to more efficiently communicate and share content with anyone, at your discretion.

Thank you for your time and support :)
  Andrew

[dlparkhurst]

I'll start by answering your first question.

(1) Yes, I do think you have some inconsistencies in your calculations. You need to reconsider your units. For the final number of moles, mol/kgw * kgs is not equal to moles in solution, you must multiply by the mass of water in the solution instead. PHREEQC SOLUTIONs have a specified mass of water, usually one kilogram, but in your case 17.4 kg. You used a solution mass of 17.4 kg.

PHREEQC converts ppm to mol/kgw by subtracting the sum of the masses of solutes from 1 kg of solution to determine the mass of water. The molalities of the elements can then be calculated from the masses of the elements and the mass of water. The results are then multiplied by the mass of water that you defined with -water.

If you look at the calculation for SOLUTION 0, the total number of moles of S in solution 0 is 1.816 and the total molality of S is 0.1045. For Ca, the total number of moles is 0.3545 and the molality is 2.040e-2. These values for the total number of moles of S and Ca in the initial solution differ from your table, but they are the values that are consistent with the solutions at the end of the simulation.

The following simplified file uses PUT to save and GET to retrieve the initial total number of moles (Basic function TOTMOL). It uses USER_PRINT to write to the output file the difference with the total number of moles in each cell between the initial and the final concentrations, and further calculates the ratio of the change in S to the change in Ca, which is very close to 1.0 in all cases. Note that the TRANSPORT calculation will achieve steady state soon after one pore volume has passed through the column. I used 24 steps to allow for two pore volumes. Also, I think "constant flux" or "flux flux" make more sense for the boundary conditions. INCREMENTAL_REACTIONS has no effect on the TRANSPORT calculations.

You can convert PHREEQC results to mol/kgs if necessary (mass of solution = RHO * SOLN_VOL), but if you are going to use PHREEQC, it is easier to work with its natural unit mol/kgw and mass of water [Basic function TOT("water")].

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
USER_PRINT
10 IF GET(1) <= 0 THEN PUT(TOTMOL("S"), 1)
20 IF GET(2) <= 0 THEN PUT(TOTMOL("Ca"), 2)
END
SOLUTION 1-12 Initial solution in the RO module
temp 25
pH 7 charge
-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

#linear_permeate
    #Effluent module 1:
        #Estimated CF: 1.118E0
        #Estimated solution mass: 15.545180409311685
END

TRANSPORT
-cells 12
-shifts 24
-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 flux flux
-punch_cells 1-12
-print_cells           1-12
-punch_frequency       24
-print_frequency        1
#-correct_disp

USER_PRINT
10 d_S = GET(1) - TOTMOL("S")
20 d_Ca = GET(2) - TOTMOL("Ca")
30 PRINT "Delta S, moles:   ", d_S
40 PRINT "Delta Ca, moles:  ", d_Ca
50 PRINT "Delta S/Delta Ca:", STR_F$(d_S/d_Ca, 20, 15)

USER_GRAPH 1
    -headings               dist H2O,kg Cl,mol/kgw
    -axis_titles            "Distance, meters" "Kilograms of water" "Molality"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X DIST
20 GRAPH_Y TOT("water")
30 GRAPH_SY TOT("Cl")
  -end
    -active                 true


[dlparkhurst]

(2) You are writing only step 51 to the selected output file. The d_Gypsum is the change from step 50 to step 51, whereas the value for Gypsum is the total amount of Gypsum precipitated in the cell over the entire simulation.

Here is a file that prints only cell 6 to the selected output file, but it prints each of the 24 transport steps. It takes a few steps before Gypsum precipitates, but you will find that if you sum the d_Gypsum value for a selected step with the Gypsum value for the previous step, you will arrive at the Gypsum value for the selected step.

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
SOLUTION 1-12 Initial solution in the RO module
temp 25
pH 7 charge
-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

#linear_permeate
    #Effluent module 1:
        #Estimated CF: 1.118E0
        #Estimated solution mass: 15.545180409311685
END
SELECTED_OUTPUT 2
-file red_sea.txt
-saturation_indices Gypsum
-equilibrium_phases Gypsum
-time true
-solution
-step
TRANSPORT
-cells 12
-shifts 24
-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 flux flux
-punch_cells 6
-print_cells            6
-print_cells            1-12
-punch_frequency       1
-print_frequency        1
END


[dlparkhurst]

(3) Unexpected mineral precipitation trend

O ye of little faith. There is nothing wrong with the calculation.

The simple answer is that you are removing less water at each successive cell. In a thought experiment, if you removed zero water from a cell, you would precipitate zero gypsum.

The more complicated answer is that you remove some water, and gypsum precipitates to a saturation index of zero, but the fraction of the bare ions Ca+2 and SO4-2, the activities of all species, the activity coefficients, and the activity of water are all changing. The net effects of all of these changes are not obvious.

[dlparkhurst]

(4) I'm not that good at at the stagnant zone calculations. If you provide a simple example as a .pqi file, I will look at it, but no promises.

Make sure you look at the section Transport in Dual Porosity Media in the 1999 manual and example 13 in the 2013 manual.

[freiburgermsu]

Dear Dr. Parkhurst,

I reply to each of your suggestions, respectively, in the following sections.

1) I just executed your PQI file in iPHREEQC and I see the mole quantities that you cite and the 1:1 ratio of ion exchange during the precipitation of Gypsum. I have a few subsequent questions.

a)
How are the 1.82 moles SO4 and 0.35453 moles Ca calculated? Here is my logic, with the liters of solution that is provided by the PHREEQC output file from your PQI file:
SO4: 9500 (ppm SO4 = mg SO4/L solution) / [(96.04 g SO4/1 mole SO4) * (1000 mg / 1 g)] * (17.67 L solution) = 1.748 moles SO4
Ca: 774 (ppm Ca = mg Ca/L solution) / [(40.08 g Ca/1 mole Ca) * (1000 mg / 1 g)] * (17.67 L solution) = 0.3412 moles Ca

What is incorrect in this calculation?

b)
Why do you believe that "constant flux" or "flux flux" are more suitable for representing RO than "constant constant"? Leakage at the boundaries seems like it could undesirably compromise mass balance, although, the end of an RO module may indeed have some backflow. What is the significance of these boundary conditions upon the geochemical predictions?

2) I executed this PQI file as well, and I see the agreement between the Gypsum and d_Gypsum moles. Thank you for articulating the discrepancy and providing the educational PQI file.

3) Does the data from the SELECTED_OUTPUT of my file depict the change in Gypsum moles at each distance instead of the net accumulation of Gypsum over all preceeding distances? We intend for the latter, where the mineral moles is ever-increasing, albeit at perhaps a diminishing rate. The Gypsum column from your SELECTED_OUTPUT in question (2) depicts an ever-increasing moles, which suggests that our TRANSPORT or SELECTED_OUTPUT blocks are incorrect. What should be changed in these blocks to depict the total accumulation of minerals over the distance?

4) The Exchange Factor parameter is clearly and concisely defined in the 1999 manual, with a mathematical explanation.

The Example 13 input files appear to use the MIX block instead of defining a set of mobile cells [1-12], a space cell [13], and then a set of immobile cells [14-25], which we implement in this example (https://github.com/freiburgermsu/ROSSpy/blob/main/examples/sensitivity_analyses/single_vs_dual_domain/2021-10-27-ROSSpy-red_sea-transport-pitzer-scaling-all_distance-LinPerm-2/input.pqi). Is the MIX method preferential, or more accurate, than our method of the sets of cells?

I greatly appreciate your detailed answers, suggestions, and context in the PHREEQC and BASIC calculations. I can elaborate anything at your request.

Thank you :)
  Andrew

[dlparkhurst]

No Dr.

(a) What is incorrect in this calculation?

As my high school chemistry teacher Mr. Mullin taught us, "Follow the units". The units must cancel to arrive at the proper final units. If you are working with brines you need to be careful.

9500 (ppm SO4 = mg SO4/L solution) Wrong!

9500 ppm = 9500 mg SO4 per kg solution.

9500 mg/kgs / (1000 mg/g * 96.04 g/mol SO4)  =  0.0989 mol/kgs

If you multiply by the number of kilograms of solution you will have the total number of moles in that many kilograms of solution. mol/kgs is absolutely independent of temperature and pressure.

If you have some number of liters of solution, then you need to multiply by the density and the number of liters to get the moles in that many liters of solution. mol/L will vary with temperature and pressure.

PHREEQC uses solutions based on kilograms of water, which effectively independent of temperaure and pressure (you could have a little bit of hydrolosis to create or consume some H2O with a closed-system T/P change). To convert from mol/kgs to mol/kgw, the sum of the solutes are subtracted from 1 kg to arrive at the mass of water.

Code: [Select]

1 - (Cl + Na + SO4 + ...) * 1e-6 kg/mg
1 - (24756 + 16417.2 + 9500 ...) * 1e-6 = kg H2O

As I said in the previous post, knowing the moles of a solute in a kilogram of solution and the mass of water in a kilogram of solution, you can convert concentrations to moles/kgw. When you use -water in the SOLUTION definition, you are defining the mass of water in the solution (not mass of solution, not liters of solution). So, if you multiply molality by the number of kilograms of water in solution, you will get the total number of moles of a solute in a solution with that many kilograms of water.

So, for brines mol/kgs, mol/L, and mol/kgw will be distinctly different. Further, the number of moles of a solute in solution will differ depending if you have 1 kg of solution, 1 L of solution, or a solution with 1 kg of water. You seem to want to use mass of solution, but you will need to make conversions to be consistent with PHREEQC which is based on the mass of water in a solution. Either convert mol/kgs to mol/kgw, convert mol/kgw to mol/kgs, or just use PHREEQC units, which will be moles, moles/kgw, and kgw.

(b) You are modeling advection and dispersion when you use TRANSPORT. Water is flowing out of the column, and there is dispersion within and at the ends of the column. Constant and constant would mean that you have an infinite reservoir at the inlet and an infinite reservoir at the outlet, so that any advection and dispersion at the ends of the column makes no difference in the reservoir concentrations. Flux flux would indicate water flowing into the column at one concentration and out at the concentration determined by the reactions.  Figures 2 and 3 show of the 1999 manual show the difference between flux and constant at the inlet. In an advection dominated system, there probably is not much difference between constant or flux at the outlet, but conceptually, flux makes more sense to me for a flow-through reverse-osmosis column.

[dlparkhurst]

(3) The selected output gives the total number of moles in each cell at each time interval that is printed. The moles are not summed from different cells.

You can accumulate the moles of gypsum in the column with USER_PUNCH or USER_GRAPH, although you may to post-process the selected output file to get exactly what you want.

Here is an example with USER_GRAPH.

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
SOLUTION 1-12 Initial solution in the RO module
temp 25
pH 7 charge
-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

#linear_permeate
    #Effluent module 1:
        #Estimated CF: 1.118E0
        #Estimated solution mass: 15.545180409311685
END
USER_GRAPH 1
    -headings               time Total_gyp
    -axis_titles            "Time, seconds" "Gypsum in column, moles" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 IF CELL_NO = 1 THEN PUT(0, 1)
 20 PUT(GET(1) + EQUI("Gypsum"), 1)
 30 IF CELL_NO < 12 THEN GOTO 100
 40 GRAPH_X TOTAL_TIME
 50 GRAPH_Y GET(1)
100 END
  -end
    -active                 true

TRANSPORT
-cells 12
-shifts 24
-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 flux flux
-punch_cells 1-12
-punch_frequency       1
-print_frequency        24
END

(4) The exchange factor applies to a single stagnant cell for each mobile cell. You will get more numerical accuracy if you use multiple stagnant cells for each mobile cell, but at the cost of considerably more complexity in setting MIX definitions and data handling. I would use the simplified version until you can demonstrate the need for more complexity.

[freiburgermsu]

Dear David,

1)
a) I reproduced the proper moles after converting (L of solution) -> (kg of solution) with the solution density from the output of your PQI. I also updated my ICE Table, which revealed a 1:1 ratio of concentration changes. Thank you for your detailed explanation!

b) The comparison of Figures 2 & 3 in the 1999 manual suggests that the inlet boundary condition has a minor influence upon the elemental concentrations. Could a steady stream of feed water through the RO module be conceptually equivalent to an infinite reservoir, since the feed supply is inexhaustible and is assumed to have a constant concentration? I concur with your assessment of flux at the outlet of an RO module.

3) The SELECTED_OUTPUT can then, essentially, print the accumulation of precipitate within a single cell over all time, however, post-processing is required to calculate the accumulation of precipitate over all cells at a single time?

Would the precipitation decrease in the final cells of this example figure(https://github.com/freiburgermsu/ROSSpy/blob/main/examples/scaling/2021-10-27-ROSSpy-red_sea-transport-pitzer-scaling-all_distance-LinPerm/all_minerals.svg) be explained by the rate of precipitation exceeding the concentration rate, thus, the precipitate per cell begins to decrease as the [Ca] and [SO4] concentrations decrease? This interpretation, however, is unsupported by the SELECTED_OUTPUT for this simulation (https://github.com/freiburgermsu/ROSSpy/blob/main/examples/scaling/2021-10-27-ROSSpy-red_sea-transport-pitzer-scaling-all_distance-LinPerm/selected_output.pqo), where only the [Ca] decreases, while the [SO4] increases over these final cells of the transport column. What explains this diverging of the ionic concentrations?

4) We -- and likely our reviewers -- certainly prefer simplicity. We will continue to troubleshoot the dual porosity calculations.

My mentors and I are primarily Civil Engineers, thus geochemistry is slightly beyond our knowledge domain. Would you like to be a co-author in the publication of this software in 1-2 months? We would benefit from your geochemical expertise and perspective, particularly through understanding the fundamentals of PHREEQC and the geochemical reactive transport of RO desalination that may not be completely captured by my 4 inquiries. I can share the manuscript package with you in a few weeks, at your discretion.

Thank you for your assistance :)
  Andrew

[dlparkhurst]

(3) Yes, you have to either write a USER_PUNCH script to accumulate among cells or post-process cell-by-cell selected output results.

As for the diverging concentrations, maybe this is a better explanation (or not). Equilibrium with gypsum requires the ion activity product to be constant [IAP = a(Ca+2)a(SO4-2)a(H2O)^2]. Not correct in detail, but if you assume activity of water of 1, and activity equals concentration, then Ca*SO4 = K.

The concentration of SO4 is larger than Ca, so assume that precipitation of gypsum has little effect on the SO4 concentration. Then if the solution is concentrated by 1 percent, about 2 percent of the resulting Ca concentration must precipitate IAP ~ (1.01SO4)*[0.98*1.01Ca] ~ K, where 0.98 (1/1.01^2) is the fraction of Ca remaining in solution after precipitation. Regardless of the details, the idea is that the precipitation of gypsum is approximately proportional to the Ca concentration. As the concentration of Ca decreases through the cells, the amount of gypsum precipitation also decreases.

Concerning authorship, I have responded by message.

[freiburgermsu]

Dear David,

1)
a) I updated the ICE table, and I noticed a slight discrepancy. The change between the [Ca] and [SO4] are identical, however, the moles of Gypsum precipitation deviates slightly from the changes in the ionic moles:

{'change': {'Ca': -0.018565196026999964, 'Gypsum': 0.019614, 'S': -0.018567569229999803},
 'final': {'Ca': 0.335957721232, 'Gypsum': 0.019614, 'S': 1.7969614040000002},
 'initial': {'Ca': 0.354522917259, 'Gypsum': 0, 'S': 1.81552897323}}

Is this ~7% deviation between d_[Ca] and d_[SO4] versus d_Gypsum within the error of the calculations, or is summing the d_Gypsum column in the SELECTED_OUTPUT file for this simulation (https://github.com/freiburgermsu/ROSSpy/blob/main/examples/other/ICE_table/selected_output.pqo) not the appropriate approach for determining the moles of Gypsum that precipitate in a timestep of the simulation? I confirmed that the sum of d_Gypsum columns were identical for all timestep shifts after SOLUTION 0 completely passes through the column.

3) Your explanation is very helpful. The approximation, to summarize, is that the least concentrated ion in a precipitation equilibrium is reduced in concentration with equal magnitude to all equilibrium ions collectively? Why is this approximation necessary relative to simultaneously reducing each equilibrium ion by its stoichiometric amount -- i.e. reduce [SO4] and [Ca] each by 1% instead of not reducing [SO4] and reducing [Ca] by 2%? I suppose that these mechanics are only important for intermediary predictions of concentrations, since all instances lead to the same final ionic concentrations.


Thank you for your assistance :)
  Andrew

[dlparkhurst]

(1) PHREEQC does not report the dispersive/diffusive fluxes, so it is difficult to calculate the mole balance if you have dispersion or diffusion. You would have to numerically estimate the derivatives in the transport equations if you want a complete mole balance.

Here are a couple of checks on mole balance.

(a) If you set -dispersion 0 (default) and -diffusion_coefficient to 0 (default 0.3e-9) in TRANSPORT, then you have a purely advective system. In that case the mass balance for cell n at time j is as follows:

Code: [Select]

TOTMOL(n-1,j-1) = TOTMOL(n,j) + EQUI_DELTA(n,j)

(b) With a dispersive system, if you run 6 shifts, corresponding to half of one pore volume, virtually all of the inflow will still be in the column. The sum of the flux in equals the sum of  the dissolved S and gypsum in all cells. The following script makes that calculation.

Code: [Select]

SOLUTION 1-12 Initial solution in the RO module
temp 25
pH 7 charge
-water 17.378153556381264
END
KNOBS
-conv 1e-10
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
USER_PRINT
10 PUT(TOTMOL("S"), 1)
20 PUT(TOTMOL("Ca"), 2)
END
USER_PRINT
10 REM
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

#linear_permeate
    #Effluent module 1:
        #Estimated CF: 1.118E0
        #Estimated solution mass: 15.545180409311685
END
USER_PRINT
10 IF CELL_NO = 1 THEN PUT(0, 3)  # sum EQUI
20 IF CELL_NO = 1 THEN PUT(0, 4)  # sum TOTMOL
30 PUT(GET(3) + EQUI("Gypsum"), 3)
40 PUT(GET(4) + TOTMOL("S"), 4)
50 IF (CELL_NO < 12) THEN GOTO 200
60 t_S_in = GET(1)*STEP_NO
70  PRINT "Flux of S into column, moles: ", STR_F$(t_S_in, 20, 10)
80  PRINT "Gypsum in column, moles:      ", STR_F$(GET(3), 20, 10)
90  PRINT "Dissolved S, moles:           ", STR_F$(GET(4), 20, 10)
100 PRINT "Flux - Gypsum - dissolved:    ", STR_F$(t_S_in - GET(3) - GET(4), 20, 10)
200 END

TRANSPORT
-cells 12
-shifts 6
-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 flux flux
-punch_cells 1-12
-print_cells 1-12
-punch_frequency       1
-print_frequency        6
-disp                   0.008
END

3. I was only trying to answer your question as simply as possible. There are probably better ways, but in the end, I rely on the PHREEQC calculation. Here is a batch calculation showing the evolution of your Red Sea water as water is extracted. A concentration factor of 10 results in halite supersaturation.

Code: [Select]

SOLUTION 1 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
END
REACTION 1
H2O -1
48.5625 in 200 step
END
EQUILIBRIUM_PHASES 1
Gypsum 0 0
END
USER_GRAPH 1
    -headings               factor Cl S Ca Gyp
    -axis_titles            "Concentration factor" "Log molality" "Gypsum, moles"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 GRAPH_X (1 / TOT("water"))
 20 GRAPH_Y log10(TOT("Cl"))
 30 GRAPH_Y log10(TOT("S"))
 40 GRAPH_Y log10(TOT("Ca"))
 50 GRAPH_SY (EQUI("Gypsum"))
100 END
  -end
    -active                 true
END
INCREMENTAL_REACTIONS
USE solution 1
USE equilibrium_phases 1
USE reaction 1
END


[freiburgermsu]

Dear David,

1) I did not consider the affects of diffusion fluxes. I appreciate your explanations and the PQI files that exemplify the mole balance in each system. These are very helpful for my understanding.

3) I clearly see the trends in ionic concentrations and Gypsum moles that you described. I will trust the results from PHREEQ.

I just posted an inquiry about the use of PHREEQpy in unix operating systems to this forum (https://phreeqcusers.org/index.php/topic,1864.0.html), although, I am not certain whether you are familiar with the mechanics of this  PHREEQ iteration. I will, nevertheless, contact you again as more questions emerge in the coming weeks.

Thank you for your assistance :)
 Andrew