PHREEQC R package

[FlouretAlex]

Hello everyone,

I am not sure if here is the best place to ask this but I am looking for examples using the PHREEQC R package.
Indeed I would like to model Cs sorption on clay and adjust the different parameters that specify sportion capacity of my clay using EXCHANGE_MASTER_SPECIES and SURFACE_MASTER_SPECIES.

My idea was to use FME package with PHREEQC package but I don't see how I can do this. If some of you have some script examples it would help me a lot.

Regards

[dlparkhurst]

I haven't used FME, but I assume it must be something like PEST where parameters in a template file are repeatedly replaced with values from the parameter estimation program. In concept, you would need to define SELECTED_OUTPUT/USER_PUNCH to write the values that you want to fit, and then use RunFile to run the parameter-replaced file (or RunString for a string) and calculate the SELECTED_OUTPUT/USER_PUNCH values, and then use methods to retrieve the SELECTED_OUTPUT/USER_PUNCH values to return to FME.

Sorry, I do not have access to a good R example. Here is an overly complicated example that evaporated waters by factors of 1, 2, 4, and 8 for open and closed systems and merged mineral precipitation amounts into an R data frame. It was used as part of this report https://ngwa.onlinelibrary.wiley.com/doi/10.1111/gwat.12367. Perhaps it has some pieces that would be useful to you.

Code: [Select]

#
# Runs PHREEQC speciation calculations on Brackish Water database
# Data frame bw is read from csv file
# Data frame df contains the selected output
# Data can be merged based on unique_id
# Uses packages phreeqc, data.table
# Package psych is useful for summary stats (describe)
# Timing: ~25 min to run 90,000 solutions
#

######################
# libraries
library(data.table)
library(phreeqc)
library(psych)  # optional


######################
# load database
phrLoadDatabase("C:/BrackishWater/pitzer.dat")

######################
# add selected_output definition
fSelOut <- function(vec) {
    return (
        c(
            vec,
            'SOLUTION 0                                                         ',
            'END                                                                ',
            'EQUILIBRIUM_PHASES 1;                                              ',
            '   Anhydrite  0     0;                                      ',
            '   Barite     0     0;                                      ',
            '   Calcite    0     0;                                      ',
            '   Celestite  0     0;                                      ',
            '   Chalcedony 0     0;                                      ',
            '   Gypsum     0     0;                                      ',
            '   Halite     0     0;                                      ',
            '   CO2(g)     -3.4 10;                                      ',
            'REACTION_TEMPERATURE 1;                                    ',
            '   15;                                                      ',
            'REACTION_PRESSURE 1;                                        ',
            '   1;                                                      ',
            'EQUILIBRIUM_PHASES 11;                                              ',
            '   Anhydrite  0     0;                                      ',
            '   Barite     0     0;                                      ',
            '   Calcite    0     0;                                      ',
            '   Celestite  0     0;                                      ',
            '   Chalcedony 0     0;                                      ',
            '   Gypsum     0     0;                                      ',
            '   Halite     0     0;                                      ',
            '   Glauberite 0     0;                                      ',
            '   Polyhalite 0     0;                                      ',
            '   Sylvite    0     0;                                      ',
            'REACTION_TEMPERATURE 2;                                    ',
            '   25;                                                      ',
            'REACTION_PRESSURE 2;                                        ',
            '   20;                                                      ',           
            'SELECTED_OUTPUT 1                                                  ',
            '    -file                 selected_output_1.sel            ',
            '    -reset                false                            ',
            '    -solution             true                            ',
            '    -pH                   true                            ',
            '    -temperature          true                            ',
            '    -charge_balance       true                            ',
            '    -percent_error        true                            ',
            '    -saturation_indices   Anhydrite Aragonite Barite        ',
            '                          Calcite Celestite Chalcedony    ',
            '                          CO2(g) Gypsum Halite Quartz      ',
            'USER_PUNCH                                                ',
            '-heading charge_balance_eqL            sum_ions                    ',
            '-heading Alk_eqL        B_molL         Ba_molL        Br_molL;     ',
            '-heading C_molL         Ca_molL        Cl_molL        Fe_molL;     ',
            '-heading K_molL         Li_molL        Mg_molL        Mn_molL;     ',
            '-heading Na_molL        SO4_molL       Si_molL        Sr_molL; ',
            '-heading BOH3_act       Ba_act         Br_act         HCO3_act; ',
            '-heading Ca_act         Cl_act         Fe_act         K_act; ',
            '-heading Li_act         Mg_act         Mn_act         Na_act; ',
            '-heading SO4_act        H4SiO4_act     Sr_act; ',
            '-heading CO2_molL       BOH3_molL      BOH4_molL      B3O3OH4_molL;',
            '-heading B4O5OH4_molL;  ',
            '-heading density_calc   sc_25C_calc    sample_id      tds_mgL_calc;',
            '-heading soln_vol_calc  kg_H2O_calc    kg_soln_calc   mu;         ',
            '-heading Cl_Br          sar            sar_ion        sar_act; ',
            '-heading osmotic_c      osmotic_p; ',
            '10     conv = TOT("water") / SOLN_VOL ',
            '20     PUNCH CHARGE_BALANCE / SOLN_VOL ',
            '30     PUNCH CHARGE_BALANCE / (PERCENT_ERROR/100) / SOLN_VOL ',
            '40     PUNCH ALK            * conv ',
            '50     PUNCH TOT("B")       * conv ',
            '60     PUNCH TOT("Ba")      * conv ',
            '70     PUNCH TOT("Br")      * conv ',
            '80     PUNCH TOT("C")       * conv ',
            '90     PUNCH TOT("Ca")      * conv ',
            '100    PUNCH TOT("Cl")      * conv ',
            '110    PUNCH TOT("Fe")      * conv ',
            '120    PUNCH TOT("K")       * conv ',
            '130    PUNCH TOT("Li")      * conv ',
            '140    PUNCH TOT("Mg")      * conv ',
            '150    PUNCH TOT("Mn")      * conv ',
            '160    PUNCH TOT("Na")      * conv ',
            '170    PUNCH TOT("S(6)")    * conv ',
            '180    PUNCH TOT("Si")  * conv ',
            '190    PUNCH TOT("Sr")  * conv ',
            '200    PUNCH ACT("B(OH)3") ',
            '210    PUNCH ACT("Ba+2") ',
            '220    PUNCH ACT("Br-") ',
            '230    PUNCH ACT("HCO3-") ',
            '240    PUNCH ACT("Ca+2") ',
            '250    PUNCH ACT("Cl-") ',
            '260    PUNCH ACT("Fe+2")    ',
            '270    PUNCH ACT("K+")    ',
            '280    PUNCH ACT("Li+")    ',
            '290    PUNCH ACT("Mg+2")    ',
            '300    PUNCH ACT("Mn+2")  ',
            '310    PUNCH ACT("Na+")    ',
            '320    PUNCH ACT("SO4-2")    ',
            '330    PUNCH ACT("H4SiO4") ',
            '340    PUNCH ACT("Sr+2") ',
            '350    PUNCH MOL("CO2")         * conv ',
            '360    PUNCH MOL("B(OH)3")      * conv ',
            '370    PUNCH MOL("B(OH)4-")     * conv ',
            '380    PUNCH MOL("B3O3(OH)4-")  * conv ',
            '390    PUNCH MOL("B4O5(OH)4-2") * conv ',
            '400    PUNCH RHO ',
            '410    PUNCH SC/(1 + 0.021 * (TC - 25)) ',
            '420    PUNCH trim(Description) ',
            '430    kgs = SOLN_VOL * RHO ',
            '440    grams = (kgs - TOT("water"))*1000      ',
            '450    PUNCH 1000*grams/SOLN_VOL ',
            '460    PUNCH SOLN_VOL ',
            '470    PUNCH TOT("water") ',
            '480    PUNCH kgs ',
            '490    PUNCH MU ',
            '500    if TOT("Br") > 0 THEN ClBr = TOT("Cl")/TOT("Br") ',
            '510    PUNCH ClBr ',
            '520    PUNCH TOT("Na")/SQRT(TOT("Ca") + TOT("Mg")) ',
            '530    PUNCH MOL("Na+")/SQRT(MOL("Ca+2") + MOL("Mg+2")) ',
            '540    PUNCH ACT("Na+")/SQRT(ACT("Ca+2") + ACT("Mg+2")) ',
            '550    PUNCH OSMOTIC ',
            '560    REM pi = -RTln(a(H2O))/V(H2O) ',
            '570    p = -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325 ',
            '580    PUNCH p '
            )
        )
}
input <- vector()
input <- fSelOut(input)
phrRunString(input)

input2 <- paste(
'SELECTED_OUTPUT 2;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 2;                                              ',
'  -headings pH_open;                                                   ',
'  -headings pp_Anhydrite_open pp_Barite_open      pp_Calcite_open; ',   
'  -heading  pp_Celestite_open pp_Chalcedony_open  pp_CO2_open;         ',
'  -headings pp_Gypsum_open    pp_Halite_open      pp_Ferrihydrite_open;',
'  -headings pp_mgL_open       os_p_open           sample_id;           ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  60 PUNCH EQUI("CO2(g)");                                  ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH TRIM(Description)                                    '
)
phrRunString(input2)

input3 <- paste(
'SELECTED_OUTPUT 3;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 3;                                              ',
'  -headings pH_open_1x;                                                ',
'  -headings pp_Anhydrite_open_1x       pp_Barite_open_1x;              ',   
'  -headings pp_Calcite_open_1x         pp_Celestite_open_1x;         ',
'  -headings pp_Chalcedony_open_1x;                                     ',
'  -headings pp_CO2_open_1x             pp_Gypsum_open_1x;              ',       
'  -headings pp_Halite_open_1x          pp_Ferrihydrite_open_1x;        ',
'  -headings pp_mgL_open_1x             os_p_open_1x;                   ',
'  -headings sample_id;                                                 ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  60 PUNCH EQUI("CO2(g)");                                  ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH TRIM(Description)                                    '
)
phrRunString(input3)

input4 <- paste(
'SELECTED_OUTPUT 4;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 4;                                              ',
'  -headings pH_open_2x;                                                ',
'  -headings pp_Anhydrite_open_2x       pp_Barite_open_2x;              ',   
'  -headings pp_Calcite_open_2x         pp_Celestite_open_2x;         ',
'  -headings pp_Chalcedony_open_2x;                                     ',
'  -headings pp_CO2_open_2x             pp_Gypsum_open_2x;              ',       
'  -headings pp_Halite_open_2x          pp_Ferrihydrite_open_2x;        ',
'  -headings pp_mgL_open_2x             os_p_open_2x;                   ',
'  -headings sample_id;                                                 ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  60 PUNCH EQUI("CO2(g)");                                  ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH TRIM(Description)                                    '
)
phrRunString(input4)

input5 <- paste(
'SELECTED_OUTPUT 5;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 5;                                              ',
'  -headings pH_open_4x;                                                ',
'  -headings pp_Anhydrite_open_4x       pp_Barite_open_4x;              ',   
'  -headings pp_Calcite_open_4x         pp_Celestite_open_4x;         ',
'  -headings pp_Chalcedony_open_4x;                                     ',
'  -headings pp_CO2_open_4x             pp_Gypsum_open_4x;              ',       
'  -headings pp_Halite_open_4x          pp_Ferrihydrite_open_4x;        ',
'  -headings pp_mgL_open_4x             os_p_open_4x;                   ',
'  -headings sample_id;                                                 ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  60 PUNCH EQUI("CO2(g)");                                  ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH TRIM(Description)                                    '
)
phrRunString(input5)

input6 <- paste(
'SELECTED_OUTPUT 6;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 6;                                              ',
'  -headings pH_open_8x;                                                ',
'  -headings pp_Anhydrite_open_8x       pp_Barite_open_8x;              ',   
'  -headings pp_Calcite_open_8x         pp_Celestite_open_8x;         ',
'  -headings pp_Chalcedony_open_8x;                                     ',
'  -headings pp_CO2_open_8x             pp_Gypsum_open_8x;              ',       
'  -headings pp_Halite_open_8x          pp_Ferrihydrite_open_8x;        ',
'  -headings pp_mgL_open_8x             os_p_open_8x;                   ',
'  -headings sample_id;                                                 ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  60 PUNCH EQUI("CO2(g)");                                  ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH TRIM(Description)                                    '
)
phrRunString(input6)

input7 <- paste(
'SELECTED_OUTPUT 7;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 7;                                              ',
'  -headings pH_open_16x;                                               ',
'  -headings pp_Anhydrite_open_16x       pp_Barite_open_16x;            ',   
'  -headings pp_Calcite_open_16x         pp_Celestite_open_16x;         ',
'  -headings pp_Chalcedony_open_16x;                                    ',
'  -headings pp_CO2_open_16x             pp_Gypsum_open_16x;            ',       
'  -headings pp_Halite_open_16x          pp_Ferrihydrite_open_16x;      ',
'  -headings pp_mgL_open_16x             os_p_open_16x;                 ',
'  -headings sample_id;                                                 ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  60 PUNCH EQUI("CO2(g)");                                  ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH TRIM(Description)                                    '
)
phrRunString(input7)

input12 <- paste(
'SELECTED_OUTPUT 12;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 12;                                              ',
'  -headings pH_closed;                                                 ',
'  -headings pp_Anhydrite_closed pp_Barite_closed      pp_Calcite_closed;     ',   
'  -heading  pp_Celestite_closed pp_Chalcedony_closed;                        ',
'  -headings pp_Gypsum_closed    pp_Halite_closed      pp_Ferrihydrite_closed;',
'  -headings pp_mgL_closed       os_p_closed           sample_id;       ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH TRIM(Description)                                    '
)
phrRunString(input12)

input13 <- paste(
'SELECTED_OUTPUT 13;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 13;                                              ',
'  -headings pH_closed_1x;                                              ',
'  -headings pp_Anhydrite_closed_1x       pp_Barite_closed_1x;      ',   
'  -headings pp_Calcite_closed_1x;                                      ',
'  -headings pp_Celestite_closed_1x       pp_Chalcedony_closed_1x;      ',
'  -headings pp_Gypsum_closed_1x;                                       ',       
'  -headings pp_Halite_closed_1x          pp_Ferrihydrite_closed_1x;    ',
'  -headings pp_mgL_closed_1x             os_p_closed_1x;               ',
'  -headings sample_id;                                                 ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH TRIM(Description)                                    '
)
phrRunString(input13)
                                   

[Jeonghwan Hwang]

Hi, my name is Jeonghwan Hwang and I studied about Cs sorption on illite during my Ph.D course.

I usually  used iPhreeqc with matlab to optimize the modeling parameter.
I think that the Pest+Phreeqc is one of the goog approach, too.

If you want to more information about iPhreeqc and sorption, please answer me.

Thank you

[dlparkhurst]

I guess there is a limit on lines, here is some more of the script.

Code: [Select]

input14 <- paste(
'SELECTED_OUTPUT 14;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 14;                                              ',
'  -headings pH_closed_2x;                                              ',
'  -headings pp_Anhydrite_closed_2x       pp_Barite_closed_2x;      ',   
'  -headings pp_Calcite_closed_2x;                                      ',
'  -headings pp_Celestite_closed_2x       pp_Chalcedony_closed_2x;      ',
'  -headings pp_Gypsum_closed_2x;                                       ',       
'  -headings pp_Halite_closed_2x          pp_Ferrihydrite_closed_2x;    ',
'  -headings pp_mgL_closed_2x             os_p_closed_2x;               ',
'  -headings sample_id;                                                 ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH TRIM(Description)                                    '
)
phrRunString(input14)

input15 <- paste(
'SELECTED_OUTPUT 15;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 15;                                              ',
'  -headings pH_closed_4x;                                              ',
'  -headings pp_Anhydrite_closed_4x       pp_Barite_closed_4x;      ',   
'  -headings pp_Calcite_closed_4x;                                      ',
'  -headings pp_Celestite_closed_4x       pp_Chalcedony_closed_4x;      ',
'  -headings pp_Gypsum_closed_4x;                                       ',       
'  -headings pp_Halite_closed_4x          pp_Ferrihydrite_closed_4x;    ',
'  -headings pp_mgL_closed_4x             os_p_closed_4x;               ',
'  -headings sample_id;                                                 ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH TRIM(Description)                                    '
)
phrRunString(input15)

input16 <- paste(
'SELECTED_OUTPUT 16;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 16;                                              ',
'  -headings pH_closed_8x;                                              ',
'  -headings pp_Anhydrite_closed_8x       pp_Barite_closed_8x;    ',   
'  -headings pp_Calcite_closed_8x;                                      ',
'  -headings pp_Celestite_closed_8x       pp_Chalcedony_closed_8x;      ',
'  -headings pp_Gypsum_closed_8x;                                       ',       
'  -headings pp_Halite_closed_8x          pp_Ferrihydrite_closed_8x;    ',
'  -headings pp_mgL_closed_8x             os_p_closed_8x;               ',
'  -headings sample_id;                                                 ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH TRIM(Description)                                    '
)
phrRunString(input16)

input17 <- paste(
'SELECTED_OUTPUT 17;                                          ',
'   -active false;                                          ',
'END;                                                        ',
'USER_PUNCH 17;                                              ',
'  -headings pH_closed_16x;                                             ',
'  -headings pp_Anhydrite_closed_16x       pp_Barite_closed_16x;        ',   
'  -headings pp_Calcite_closed_16x;                                     ',
'  -headings pp_Celestite_closed_16x       pp_Chalcedony_closed_16x;    ',
'  -headings pp_Gypsum_closed_16x;                                      ',       
'  -headings pp_Halite_closed_16x          pp_Ferrihydrite_closed_16x;  ',
'  -headings pp_mgL_closed_16x             os_p_closed_16x;             ',
'  -headings sample_id;                                                 ',
'   5 PUNCH -LA("H+");                                        ',
'  10 PUNCH EQUI("Anhydrite");                                ',
'  20 PUNCH EQUI("Barite");                                  ',
'  30 PUNCH EQUI("Calcite");                                ',
'  40 PUNCH EQUI("Celestite");                              ',
'  50 PUNCH EQUI("Chalcedony");                              ',
'  70 PUNCH EQUI("Gypsum");                                  ',
'  80 PUNCH EQUI("Halite");                                  ',
'  90 PUNCH TOTMOL("Fe");                                    ',
' 100 pp = EQUI("Anhydrite")*GFW("CaSO4")*1000;                  ',
' 110 pp = pp + EQUI("Barite")*GFW("BaSO4")*1000;                ',
' 120 pp = pp + EQUI("Calcite")*GFW("CaCO3")*1000;              ',
' 130 pp = pp + EQUI("Celestite")*GFW("SrSO4")*1000;            ',
' 140 pp = pp + EQUI("Chalcedony")*GFW("SiO2")*1000;            ',
' 160 pp = pp + EQUI("Gypsum")*GFW("CaSO4(H2O)2")*1000;          ',
' 170 pp = pp + EQUI("Halite")*GFW("NaCl")*1000;                ',
' 175 pp = pp + TOTMOL("Fe")*GFW("Fe(OH)3")*1000;                ',
' 180 PUNCH pp;                                                  ',
' 190 PUNCH -8.314*TK*LOG(ACT("H2O"))/(VM("H2O")/1e6)/101325;    ',
' 200 PUNCH , TRIM(Description)                                    '
)
phrRunString(input17)



# EQUILIBRIUM_PHASES   1  # open system
# EQUILIBRIUM_PHASES   11 # closed system
#
# REACTION_TEMPERATURE 1  # 15 C
# REACTION_TEMPERATURE 11 # 15 C
# REACTION_TEMPERATURE 2  # 25 C
#
# REACTION_PRESSURE    1  # 1 atm
# REACTION_PRESSURE    2  # 20 atm
#
# Selected_output number         Calculation
# 1                              Speciation calculation
# 2                              Open system pp
# 3                              Open system pp  1x
# 4                              Open system pp  2x
# 5                              Open system pp  4x
# 6                              Open system pp  8x
# 7                              Open system pp  16x
# 12          closed system pp
# 13          closed system pp  1x
# 14          closed system pp  2x
# 15          closed system pp  4x
# 16          closed system pp  8x
# 17          closed system pp  16x
######################
# pull in data
bw_orig <- read.csv('C:/BrackishWater/2015_12_31/MajorIons_FILT_nodesc_2015_12_31.csv')
#bw_orig <- read.csv('C:/BrackishWater/test.csv')
#uid <- bw_orig[c("sample_id")]
#bw_orig$dup <- duplicated(uid$sample_id)
#bw_unique <- subset(bw_orig, ! ( bw_orig$dup == 'TRUE' ))

names(bw_orig)[names(bw_orig)=="unique_site_ID"] <- "sample_id"

bw <- bw_orig[,c("sample_id",
             "Temp_C",
             "Dens_gmL",
             "ph_use",
             "Alk_mgL",
             "B_ugL",
             "Ba_ugL",
             "Br_mgL",
             "Ca_mgL",
             "Cl_mgL",
             "Fe_mgL",
             "K_mgL",
             "Li_ugL",
             "Mg_mgL",
             "Mn_ugL",
             "Na_mgL",
             "Si_mgL",
             "SO4_mgL",
             "Sr_ugL",
             "SC_pref",
             "TDS_mgL_use")]
             

######################
# Censor data 
# Temperature
bw <- subset(bw, bw$Temp_C < 100.)
bw <- subset(bw, bw$Temp_C > 0.0)             
# pH/alkalinity
bw <- subset(bw, bw$ph_use < 14.0)
bw <- subset(bw, bw$ph_use > 0.0)
bw$Alk_mgL[bw$Alk_mgL < 0.0] <- 0
bw$Alk_mgL[bw$ph_use <= 4.5 & bw$Alk_mgL > 0] <- 0.0
#bw <- subset(bw, ! (bw$ph_use < 4.5 & bw$Alk_mgL > 0))
#bw <- subset(bw, bw$Alk_mgL < 10000.0)
# default T
bw$Temp_C[is.na(bw$Temp_C)] <- 15.0 
# Set missing to 0
bw <- subset(bw, ! is.na(bw$Alk_mgL) )
bw <- subset(bw, ! is.na(bw$Ca_mgL) )
bw <- subset(bw, ! is.na(bw$Cl_mgL) )
bw <- subset(bw, ! is.na(bw$Mg_mgL) )
bw <- subset(bw, ! is.na(bw$Na_mgL) )
bw <- subset(bw, ! is.na(bw$SO4_mgL) )
bw$B_ugL[is.na(bw$B_ugL)] <- 0
bw$Ba_ugL[is.na(bw$Ba_ugL)] <- 0
bw$Br_mgL[is.na(bw$Br_mgL)] <- 0
bw$Fe_mgL[is.na(bw$Fe_mgL)] <- 0
bw$K_mgL[is.na(bw$K_mgL)] <- 0
bw$Li_ugL[is.na(bw$Li_ugL)] <- 0
bw$Mn_ugL[is.na(bw$Mn_ugL)] <- 0
bw$Si_mgL[is.na(bw$Si_mgL)] <- 0
bw$Sr_ugL[is.na(bw$Sr_ugL)] <- 0
# Density
#bw$Dens_gmL[is.na(bw$Dens_gmL)] <- 1 + 0.03476 * bw$Cl_mgL / (35.453 * 1000.0) + 0.13593 * bw$SO4_mgL / (64.064 * 1000.0)
bw$Dens_gmL[is.na(bw$Dens_gmL) | bw$Dens_gmL == 0] <- 1
# remove NA
#bw <- na.omit(bw)

######################
# Format input for PHREEQC
bw$input <- paste('SOLUTION    ', rownames(bw), bw$sample_id, ';',
            '  units mg/L;',
            '  density   ', bw$Dens_gmL,  '             ;',
            '  temp      ', bw$Temp_C,    '             ;',
            '  pH        ', bw$ph_use,    '             ;',
            '  Alkalinity', bw$Alk_mgL,   '             ;',
            '  B         ', bw$B_ugL,     ' ug/L        ;',
            '  Ba        ', bw$Ba_ugL,    ' ug/L        ;',
            '  Br        ', bw$Br_mgL,    '             ;',
            '  Ca        ', bw$Ca_mgL,    '             ;',
            '  Cl        ', bw$Cl_mgL,    '             ;',
            '  Fe        ', bw$Fe_mgL,    '             ;',
            '  K         ', bw$K_mgL,     '             ;',
            '  Li        ', bw$Li_ugL,    ' ug/L        ;',
            '  Mg        ', bw$Mg_mgL,    '             ;',
            '  Mn        ', bw$Mn_ugL,    ' ug/L        ;',
            '  Na        ', bw$Na_mgL,    '             ;',
            '  Si        ', bw$Si_mgL,    '             ;',
            '  S(6)      ', bw$SO4_mgL,   '             ;',
            '  Sr        ', bw$Sr_ugL,    ' ug/L        ;',
            'END',                                     ';',
            'TITLE open pp',                           ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 1',                ';',
            'USE reaction_temperature 1',              ';',
            'USE reaction_pressure 1',                 ';',
            'SELECTED_OUTPUT 1; -active false',        ';',
            'SELECTED_OUTPUT 2; -active true',         ';',
            'END',                                     ';',
            'TITLE open pp 1x',                        ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 1',                ';',
            'USE reaction_temperature 2',              ';',
            'USE reaction_pressure 2',                 ';',
            'REACTION',                                ';',
            '  H2O -1',                                ';',
            '  0',                                     ';',
            'SELECTED_OUTPUT 2; -active false',        ';',
            'SELECTED_OUTPUT 3; -active true',         ';',
            'END',                                     ';',
            'TITLE open pp 2x',                        ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 1',                ';',
            'USE reaction_temperature 2',              ';',
            'USE reaction_pressure 2',                 ';',
            'REACTION',                                ';',
            '  H2O -1',                                ';',
            '  27.75',                                 ';',
            'SELECTED_OUTPUT 3; -active false',        ';',
            'SELECTED_OUTPUT 4; -active true',         ';',
            'END',                                     ';',
            'TITLE open pp 4x',                        ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 1',                ';',
            'USE reaction_temperature 2',              ';',
            'USE reaction_pressure 2',                 ';',
            'REACTION',                                ';',
            '  H2O -1',                                ';',
            '  41.625',                                ';',
            'SELECTED_OUTPUT 4; -active false',        ';',
            'SELECTED_OUTPUT 5; -active true',         ';',
            'END',                                     ';',
            'TITLE open pp 8x',                        ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 1',                ';',
            'USE reaction_temperature 2',              ';',
            'USE reaction_pressure 2',                 ';',
            'REACTION',                                ';',
            '  H2O -1',                                ';',
            '  48.5625',                               ';',
            'SELECTED_OUTPUT 5; -active false',        ';',
            'SELECTED_OUTPUT 6; -active true',         ';',
            'END',                                     ';',
            'TITLE open pp 16x',                       ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 1',                ';',
            'USE reaction_temperature 2',              ';',
            'USE reaction_pressure 2',                 ';',
            'REACTION',                                ';',
            '  H2O -1',                                ';',
            '  52.03125',                              ';',
            'SELECTED_OUTPUT 6; -active false',        ';',
            'SELECTED_OUTPUT 7; -active true',         ';',
            'END',                                     ';',
            'TITLE closed pp',                         ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 11',               ';',
            'USE reaction_temperature 1',              ';',
            'USE reaction_pressure 1',                 ';',
            'SELECTED_OUTPUT 7; -active false',        ';',
            'SELECTED_OUTPUT 12; -active true',        ';',
            'END',                                     ';',
            'TITLE closed pp 1x',                      ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 11',               ';',
            'USE reaction_temperature 2',              ';',
            'USE reaction_pressure 2',                 ';',
            'REACTION',                                ';',
            '  H2O -1',                                ';',
            '  0',                                     ';',
            'SELECTED_OUTPUT 12; -active false',       ';',
            'SELECTED_OUTPUT 13; -active true',        ';',
            'END',                                     ';',
            'TITLE closed pp 2x',                      ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 11',                ';',
            'USE reaction_temperature 2',              ';',
            'USE reaction_pressure 2',                 ';',
            'REACTION',                                ';',
            '  H2O -1',                                ';',
            '  27.75',                                 ';',
            'SELECTED_OUTPUT 13; -active false',       ';',
            'SELECTED_OUTPUT 14; -active true',        ';',
            'END',                                     ';',
            'TITLE closed pp 4x',                      ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 11',               ';',
            'USE reaction_temperature 2',              ';',
            'USE reaction_pressure 2',                 ';',
            'REACTION',                                ';',
            '  H2O -1',                                ';',
            '  41.625',                                ';',
            'SELECTED_OUTPUT 14; -active false',       ';',
            'SELECTED_OUTPUT 15; -active true',        ';',
            'END',                                     ';',
            'TITLE closed pp 8x',                      ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 11',               ';',
            'USE reaction_temperature 2',              ';',
            'USE reaction_pressure 2',                 ';',
            'REACTION',                                ';',
            '  H2O -1',                                ';',
            '  48.5625',                               ';',
            'SELECTED_OUTPUT 15; -active false',       ';',
            'SELECTED_OUTPUT 16; -active true',        ';',
            'END',                                     ';',
            'TITLE closed pp 16x',                     ';',
            'USE solution ', rownames(bw),             ';',
            'USE equilibrium_phases 11',               ';',
            'USE reaction_temperature 2',              ';',
            'USE reaction_pressure 2',                 ';',
            'REACTION',                                ';',
            '  H2O -1',                                ';',
            '  52.03125',                              ';',
            'SELECTED_OUTPUT 16; -active false',       ';',
            'SELECTED_OUTPUT 17; -active true',        ';',
            'END',                                     ';',     
            'DELETE; -solution ', rownames(bw),        ';',
            'SELECTED_OUTPUT 17; -active false',       ';',
            'SELECTED_OUTPUT 1; -active true',         ';'
            )

######################
# Run speciation, check for errors, store result in a list (so)
so <- list(length = nrow(bw))
so2 <- list(length = nrow(bw))
so3 <- list(length = nrow(bw))
so4 <- list(length = nrow(bw))
so5 <- list(length = nrow(bw))
so6 <- list(length = nrow(bw))
so7 <- list(length = nrow(bw))
so12 <- list(length = nrow(bw))
so13 <- list(length = nrow(bw))
so14 <- list(length = nrow(bw))
so15 <- list(length = nrow(bw))
so16 <- list(length = nrow(bw))
so17 <- list(length = nrow(bw))
for (i in 1:nrow(bw))
#for (i in 1:100)
{
out <- tryCatch(
{
phrRunString(bw[i, "input"])
so[[i]]   <- phrGetSelectedOutput()$n1
so2[[i]]  <- phrGetSelectedOutput()$n2
so3[[i]]  <- phrGetSelectedOutput()$n3
so4[[i]]  <- phrGetSelectedOutput()$n4
so5[[i]]  <- phrGetSelectedOutput()$n5
so6[[i]]  <- phrGetSelectedOutput()$n6
so7[[i]]  <- phrGetSelectedOutput()$n7
so12[[i]] <- phrGetSelectedOutput()$n12
so13[[i]] <- phrGetSelectedOutput()$n13
so14[[i]] <- phrGetSelectedOutput()$n14
so15[[i]] <- phrGetSelectedOutput()$n15
so16[[i]] <- phrGetSelectedOutput()$n16
so17[[i]] <- phrGetSelectedOutput()$n17
},
error=function(cond){
message(paste("SOLUTION calculation FAILED for number ", i))
message(paste(bw[i, "input"]))
},
finally={
if (i%%1000 == 0) {
message(paste("Through SOLUTION number ", i))
}
}
)
}

######################
# Convert selected-output list to a single data frame
df = rbindlist(so)
df2 = rbindlist(so2)
df3 = rbindlist(so3)
df4 = rbindlist(so4)
df5 = rbindlist(so5)
df6 = rbindlist(so6)
df7 = rbindlist(so7)
df12 = rbindlist(so12)
df13 = rbindlist(so13)
df14 = rbindlist(so14)
df15 = rbindlist(so15)
df16 = rbindlist(so16)
df17 = rbindlist(so17)

#df$sample_id <- as.numeric(as.character(df$sample_id))
#df2$sample_id <- as.numeric(as.character(df2$sample_id))
#df3$sample_id <- as.numeric(as.character(df3$sample_id))
#df4$sample_id <- as.numeric(as.character(df4$sample_id))
#df5$sample_id <- as.numeric(as.character(df5$sample_id))
#df6$sample_id <- as.numeric(as.character(df6$sample_id))
#df7$sample_id <- as.numeric(as.character(df7$sample_id))
#df12$sample_id <- as.numeric(as.character(df12$sample_id))
#df13$sample_id <- as.numeric(as.character(df13$sample_id))
#df14$sample_id <- as.numeric(as.character(df14$sample_id))
#df15$sample_id <- as.numeric(as.character(df15$sample_id))
#df16$sample_id <- as.numeric(as.character(df16$sample_id))
#df17$sample_id <- as.numeric(as.character(df17$sample_id))

######################
# Merge columns
#df$sample_id <- sub("\\s+$", "", df$sample_id)
#add <- bw[c("sample_id")]
df = merge(bw_orig, df, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df2, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df3, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df4, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df5, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df6, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df7, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df12, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df13, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df14, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df15, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df16, by="sample_id", all.x=TRUE, all.y=FALSE)
df = merge(df, df17, by="sample_id", all.x=TRUE, all.y=FALSE)

plot(df$sc_25C_calc,df$tds_mgL_calc, xlim=c(0,2500),ylim=c(0,2500))

write.csv(df, "c:/BrackishWater/MajorIons_FILT_nodesc_2015_12_31_processed_20160122.csv")



[FlouretAlex]

Hello,

Thank you for your reply I will take a look and see what I can do. I'll post some script as example when I will be able to run it =).

Have a good day