Applications and Case Studies > Civil engineering

Desalination PHREEQpy applications -- ROSSpy

(1/3) > >>

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: ---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

--- End code ---

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: ---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

--- End code ---

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.

Navigation

[0] Message Index

[#] Next page

Go to full version