Problem with MIX when changing proportion

[FlouretAlex]

Hello Phreeqc users,

I have a problem as I don't know how to figure out a problem with to variation of the same script but not the same approach of coding it.

Script one is using M_soil/V_solution ratio for simplicity of having everything by liter to represent the experimentation.

Script 2 is using the mass of water and soil for the calculation of the exchanger and surface and volume of water for the soil solution.

My problem is that when I run script 2 the there is a difference of morality with script 1 in solution creating a big difference of cation present in the solution as you can see below.

Result of the first mix for script 1:

Code: [Select]

Reaction step 1.

Using mix 3.

Mixture 3.

  9.000e-02 Solution 0                                                        
  1.000e+00 Solution 3                                                        
-9.000e-02 Solution 100                                                        

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

Cl                1.090e-05   1.090e-05
Cs                3.350e-04   3.350e-04
K                 1.090e-05   1.090e-05
Mg                8.176e-04   8.176e-04
S                 8.176e-04   8.176e-04


Result of the first mix for script 2:

Code: [Select]

Reaction step 1.

Using mix 3.

Mixture 3.

  9.000e-02 Solution 0                                                        
  1.000e+00 Solution 3                                                        
-9.000e-02 Solution 100                                                        

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

Cl                1.922e-04   9.494e-07
Cs                3.350e-04   1.655e-06
K                 1.922e-04   9.494e-07
Mg                1.441e-02   7.121e-05
S                 1.441e-02   7.121e-05

My problem is clearly coming from the MIX data block (or maybe a unit problem with solution 0-1).

But I don't know how to deal with this. Maybe using another way to simulate the evaporation from the soil?

Here script 1

Code: [Select]

DATABASE C:\Users\International\Documents\PhD\AlexF\Modelo alex\Fitted\000_Database_PHREEQC\PKDLM_BRGM_database_phreeqc_ThermoddemV1.10_06Jun2017.dat

SURFACE 10
-equilibrate 10
-sites_units density 
Illitex 2.07E-05 97 3695.29 # 0.5%
Montx 3.80E-07 800 3695.29 # 1.1%
Kaolinitex 1.46E-06 10 3695.29 # 0.7%
# ^site/nm^2 ^SSA (m^2/g) ^M/V (g/L)
EXCHANGE  10
-equilibrate 10
Illi 6.37E-03
Mont 3.39E-02
Kao 7.19E-04

# Equilibrage spike
SOLUTION 10
 units mol/L
pH 9.32
Na 5.00E-03
K 1.60E-03
Ca 4.50E-03
Mg 1.00E-03
N(5) 1.60E-03
Cl 5.00E-03
S(6) 1.00E-03
C(4) 4.50E-03
EQUILIBRIUM_PHASES   10
Fix_H+ -9.32 HCl 10
save surface 3
save exchange 3
END

#Contaminated soil ==> RIZO
SOLUTION 3 # ==> stagnant cell
units mol/L
pH 9.32
Cl 1.00E-05
K 1.00E-05
Mg 7.50E-04
S(6) 7.50E-04
Cs 3.35E-04
use surface 3 #number of surface and exchange need to be the same as the cell here solution 3 ==> cell 3
use Exchange 3
EQUILIBRIUM_PHASES   3
Fix_H+ -9.32 HCl  10

END

SOLUTION 0-1 # inlet solution 0 # define solution 1 since PHREEQC wants at least 1 cell in the colum for transport ==> mobile cell.
 units mol/L
pH 9.32
Cl 1.00E-05
K 1.00E-05
Mg 7.50E-04
S(6) 7.50E-04

Solution 100 #for evaporation
END

MIX 3 # simulate % water evaporation of solution 3 then replaced by solution 0
0 0.09
100 -0.09
3 1
END

 SELECTED_OUTPUT
-reset false
-file C:\Users\International\Documents\PhD\AlexF\Modelo alex\Fitted\RHIZO\RHIZO_moy\RHIZOsol\Sol-S\2_Fit\test_MV\data_RHIZO_MV.sel

-state TRUE
-solution TRUE
-water TRUE
-step TRUE
-time TRUE
-high_precision TRUE

-totals Cs K Na Ca Mg Illitex Montx Kaolinitex Illi Mont Kao
-molalities IllitexOCs+0.5 IllitexOK+0.5 IllitexONa+0.5 IllitexOCa+1.5 IllitexOMg+1.5 IllitexO-0.5 IllitexOH+0.5 MontxOCs+0.5 MontxOK+0.5 MontxONa+0.5 MontxOCa+1.5 MontxOMg+1.5 MontxO-0.5 MontxOH+0.5 KaolinitexOCs+0.5 KaolinitexOK+0.5 KaolinitexONa+0.5 KaolinitexOCa+1.5 KaolinitexOMg+1.5 KaolinitexO-0.5  KaolinitexOH+0.5 CsIlli KIlli NaIlli CaIlli2 MgIlli2 CsMont KMont NaMont CaMont2 MgMont2 CsKao KKao NaKao CaKao2         

USER_PUNCH
-headings nb_cel T_exp log10(Csaq) log10(Kd) Cssurf Csexch Ksurf Kexch Nasurf Naexch Casurf Caexch Mgsurf Mgexch Csverif Csini Kini Naini Caini Mgini Csaq Kaq Naaq Caaq Mgaq
-start
10 Msoil = 18.31/1000 # soil_masse(Kg)
20 Vwater = 4.94/1000 # volume batch (L)
30 MVw = Msoil/Vwater #Kg/L
40 VwM = Vwater/Msoil #L/Kg
50 Csini = sys("Cs") # (mol/L)
60 Caaq = tot("Ca") #(mol/L)
70 Csaq = tot("Cs") #(mol/L)
80 Kaq = tot("K") #(mol/L)
90 Naaq = tot("Na") #(mol/L)
100 Mgaq = tot("Mg") #(mol/L)

110 surfCs = (surf("Cs","Illitex")+ surf("Cs","Montx")+ surf("Cs","Kaolinitex")) #(mol/L)
120 exchCs = (mol("CsIlli")+ mol("CsMont")+ mol ("CsKao")) #(mol/L)
130 surfK = (surf("K","Illitex")+ surf("K","Montx")+ surf("K","Kaolinitex")) #(mol/L)
140 exchK = (mol("KIlli")+ mol("KMont")+ mol ("KKao")) #(mol/L)
150 surfNa = (surf("Na","Illitex")+ surf("Na","Montx")+ surf("Na","Kaolinitex")) #(mol/L)
160 exchNa = (mol("NaIlli")+ mol("NaKao")+ mol ("NaMont")) #(mol/L)
170 surfCa = (surf("Ca","Illitex")+ surf("Ca","Montx")+ surf("Ca","Kaolinitex")) #(mol/L)
180 exchCa = (mol("CaIlli2")+ mol("CaKao2")+ mol ("CaMont2")) #(mol/L)
190 surfMg = (surf("Mg","Illitex")+ surf("Mg","Montx")+ surf("Mg","Kaolinitex")) #(mol/L)
200 exchMg = (mol("MgIlli2")+ mol("MgMont2")+ mol ("MgKao2")) #(mol/L)

210 Kdsurf_Cs = (surfCs/Csaq)*VwM #(L/Kg)
220 Kdexch_Cs = (exchCs/Csaq)*VwM #(L/Kg)
230 Kdsurf_k = (surfK/Kaq)*VwM #(L/Kg)
240 Kdexch_k = (exchK/Kaq)*VwM #(L/Kg)
250 Kdsurf_Na = (surfNa/Naaq)*VwM #(L/Kg)
260 Kdexch_Na = (exchNa/Naaq)*VwM #(L/Kg)
270 Kdsurf_Ca = (surfCa/Caaq)*VwM #(L/Kg)
280 Kdexch_Ca = (exchCa/Caaq)*VwM #(L/Kg)
290 Kdsurf_Mg = (surfMg/Mgaq)*VwM #(L/Kg)
300 Kdexch_Mg = (exchMg/Mgaq)*VwM #(L/Kg)

310 Cssurf = Csaq*Kdsurf_Cs #(mol/Kg)
320 Csexch = Csaq*Kdexch_Cs  #(mol/Kg)
330 Ksurf = Kaq*Kdsurf_K #(mol/Kg)
340 Kexch = Kaq*Kdexch_K  #(mol/Kg)
350 Nasurf = Naaq*Kdsurf_Na #(mol/Kg)
360 Naexch = Naaq*Kdexch_Na  #(mol/Kg)
370 Casurf = Caaq*Kdsurf_Ca #(mol/Kg)
380 Caexch = Caaq*Kdexch_Ca  #(mol/Kg)
390 Mgsurf = Mgaq*Kdsurf_Mg #(mol/Kg)
400 Mgexch = Mgaq*Kdexch_Mg  #(mol/Kg)

410 FESCs = (surf("Cs","Illitex")) #(mol/L)
420 FESK = (surf("K","Illitex")) #(mol/L)
430 EXILLICs = mol("CsIlli") #(mol/L)
440 EXILLIK = mol("KIlli") #(mol/L)
450 Csfix = surfCs + exchCs #(mol/L)
460 Kd = (Csfix/Csaq)*VwM #(L/Kg)
470 Cssorb = Csaq*Kd #(mol/Kg)
480 Illi_EX_CS = EXILLICs*VwM #(mol/Kg)
490 Illi_EX_K = EXILLIK*VwM #(mol/Kg)
500 Cs_FES = FESCs*VwM #(mol/Kg)
510 K_FES = FESK*VwM #(mol/Kg)

520 Csverif = Csaq + Cssorb*MVw #(mol/L)
530 dt= 3600 # s
540 T_R= 40821.11 # = cel renew time (s)
550 nb_cel= dt*(step_no)/ T_R #nb renewed cel
560 T_exp = dt*(step_no) #experiment time (s)

570 Kini = sys("K") # (mol/L)
580 Naini = sys("Na") # (mol/L)
590 Caini = sys("Ca") # (mol/L)
700 Mgini = sys("Mg") # (mol/L)

800 punch nb_cel T_exp log10(Csaq) log10(Kd) surfCs exchCs surfK exchK surfNa exchNa surfCa exchCa surfMg exchMg Csverif Csini Kini Naini Caini Mgini Csaq Kaq Naaq Caaq Mgaq
 -end
 
TRANSPORT
-cells 1
-lengths 0.01 # 1 cm sol
-shifts 576 # number of time_steps total time is 24 days
-time_step 3600 # each time_step (s) the MIX is performed
-flow_direction  diffusion_only
-boundary_conditions constant closed
-diffusion_coefficient 2.3E-08
-stagnant 1
-punch_cells 3
-punch_frequency 1 # sample every time step
END
Here script 2

Code: [Select]

DATABASE C:\Users\International\Documents\PhD\AlexF\Modelo alex\Fitted\000_Database_PHREEQC\PKDLM_BRGM_database_phreeqc_ThermoddemV1.10_06Jun2017.dat

SURFACE 10
-equilibrate 10
-sites_units density 
Illitex 2.07E-05 97 18.31 # 0.5%
Montx 3.80E-07 800 18.31 # 1.1%
Kaolinitex 1.46E-06 10 18.31 # 0.7%
# ^site/nm^2 ^SSA (m^2/g) ^M (g)
EXCHANGE  10
-equilibrate 10
Illi 3.16E-05
Mont 1.68E-04
Kao 3.56E-06

# Equilibrage spike
SOLUTION 10
-water 4.94E-03 #kg
 units mol/L
pH 9.32
Na 5.00E-03
K 1.60E-03
Ca 4.50E-03
Mg 1.00E-03
N(5) 1.60E-03
Cl 5.00E-03
S(6) 1.00E-03
C(4) 4.50E-03
EQUILIBRIUM_PHASES   10
Fix_H+ -9.32 HCl 10
save surface 3
save exchange 3
END

#Contaminated soil ==> RIZO
SOLUTION 3 # ==> stagnant cell
-water 4.94E-03 #kg
units mol/L
pH 9.32
Cl 1.00E-05
K 1.00E-05
Mg 7.50E-04
S(6) 7.50E-04
Cs 3.35E-04
use surface 3 #number of surface and exchange need to be the same as the cell here solution 3 ==> cell 3
use Exchange 3
EQUILIBRIUM_PHASES   3
Fix_H+ -9.32 HCl  10

END

SOLUTION 0-1 # inlet solution 0 # define solution 1 since PHREEQC wants at least 1 cell in the colum for transport ==> mobile cell.
units mol/L
pH 9.32
Cl 1.00E-05
K 1.00E-05
Mg 7.50E-04
S(6) 7.50E-04

Solution 100 #for evaporation
END

MIX 3 # simulate % water evaporation of solution 3 then replaced by solution 0
0 0.09
100 -0.09
3 1
END

SELECTED_OUTPUT
-reset false
-file C:\Users\International\Documents\PhD\AlexF\Modelo alex\Fitted\RHIZO\RHIZO_moy\RHIZOsol\Sol-S\2_Fit\test_MV\data_RHIZO_M.sel

-state TRUE
-solution TRUE
-water TRUE
-step TRUE
-time TRUE
-high_precision TRUE

-totals Cs K Na Ca Mg Illitex Montx Kaolinitex Illi Mont Kao
-molalities IllitexOCs+0.5 IllitexOK+0.5 IllitexONa+0.5 IllitexOCa+1.5 IllitexOMg+1.5 IllitexO-0.5 IllitexOH+0.5 MontxOCs+0.5 MontxOK+0.5 MontxONa+0.5 MontxOCa+1.5 MontxOMg+1.5 MontxO-0.5 MontxOH+0.5 KaolinitexOCs+0.5 KaolinitexOK+0.5 KaolinitexONa+0.5 KaolinitexOCa+1.5 KaolinitexOMg+1.5 KaolinitexO-0.5  KaolinitexOH+0.5 CsIlli KIlli NaIlli CaIlli2 MgIlli2 CsMont KMont NaMont CaMont2 MgMont2 CsKao KKao NaKao CaKao2       

USER_PUNCH
-headings nb_cel T_exp log10(Csaq) log10(Kd) Cssurf Csexch Ksurf Kexch Nasurf Naexch Casurf Caexch Mgsurf Mgexch Csverif Csini Kini Naini Caini Mgini Csaq Kaq Naaq Caaq Mgaq
-start
10 Msoil = 18.31/1000 # soil_masse(Kg)
20 Vwater = 4.94/1000 # volume batch (L)
30 MVw = Msoil/Vwater #Kg/L
40 VwM = Vwater/Msoil #L/Kg
50 Csini = sys("Cs") # (mol/L)
60 Caaq = tot("Ca") #(mol/L)
70 Csaq = tot("Cs") #(mol/L)
80 Kaq = tot("K") #(mol/L)
90 Naaq = tot("Na") #(mol/L)
100 Mgaq = tot("Mg") #(mol/L)

110 surfCs = (surf("Cs","Illitex")+ surf("Cs","Montx")+ surf("Cs","Kaolinitex")) #(mol/L)
120 exchCs = (mol("CsIlli")+ mol("CsMont")+ mol ("CsKao")) #(mol/L)
130 surfK = (surf("K","Illitex")+ surf("K","Montx")+ surf("K","Kaolinitex")) #(mol/L)
140 exchK = (mol("KIlli")+ mol("KMont")+ mol ("KKao")) #(mol/L)
150 surfNa = (surf("Na","Illitex")+ surf("Na","Montx")+ surf("Na","Kaolinitex")) #(mol/L)
160 exchNa = (mol("NaIlli")+ mol("NaKao")+ mol ("NaMont")) #(mol/L)
170 surfCa = (surf("Ca","Illitex")+ surf("Ca","Montx")+ surf("Ca","Kaolinitex")) #(mol/L)
180 exchCa = (mol("CaIlli2")+ mol("CaKao2")+ mol ("CaMont2")) #(mol/L)
190 surfMg = (surf("Mg","Illitex")+ surf("Mg","Montx")+ surf("Mg","Kaolinitex")) #(mol/L)
200 exchMg = (mol("MgIlli2")+ mol("MgMont2")+ mol ("MgKao2")) #(mol/L)

210 Kdsurf_Cs = (surfCs/Csaq)*VwM #(L/Kg)
220 Kdexch_Cs = (exchCs/Csaq)*VwM #(L/Kg)
230 Kdsurf_k = (surfK/Kaq)*VwM #(L/Kg)
240 Kdexch_k = (exchK/Kaq)*VwM #(L/Kg)
250 Kdsurf_Na = (surfNa/Naaq)*VwM #(L/Kg)
260 Kdexch_Na = (exchNa/Naaq)*VwM #(L/Kg)
270 Kdsurf_Ca = (surfCa/Caaq)*VwM #(L/Kg)
280 Kdexch_Ca = (exchCa/Caaq)*VwM #(L/Kg)
290 Kdsurf_Mg = (surfMg/Mgaq)*VwM #(L/Kg)
300 Kdexch_Mg = (exchMg/Mgaq)*VwM #(L/Kg)

310 Cssurf = Csaq*Kdsurf_Cs #(mol/Kg)
320 Csexch = Csaq*Kdexch_Cs  #(mol/Kg)
330 Ksurf = Kaq*Kdsurf_K #(mol/Kg)
340 Kexch = Kaq*Kdexch_K  #(mol/Kg)
350 Nasurf = Naaq*Kdsurf_Na #(mol/Kg)
360 Naexch = Naaq*Kdexch_Na  #(mol/Kg)
370 Casurf = Caaq*Kdsurf_Ca #(mol/Kg)
380 Caexch = Caaq*Kdexch_Ca  #(mol/Kg)
390 Mgsurf = Mgaq*Kdsurf_Mg #(mol/Kg)
400 Mgexch = Mgaq*Kdexch_Mg  #(mol/Kg)

410 FESCs = (surf("Cs","Illitex")) #(mol/L)
420 FESK = (surf("K","Illitex")) #(mol/L)
430 EXILLICs = mol("CsIlli") #(mol/L)
440 EXILLIK = mol("KIlli") #(mol/L)
450 Csfix = surfCs + exchCs #(mol/L)
460 Kd = (Csfix/Csaq)*VwM #(L/Kg)
470 Cssorb = Csaq*Kd #(mol/Kg)
480 Illi_EX_CS = EXILLICs*VwM #(mol/Kg)
490 Illi_EX_K = EXILLIK*VwM #(mol/Kg)
500 Cs_FES = FESCs*VwM #(mol/Kg)
510 K_FES = FESK*VwM #(mol/Kg)

520 Csverif = Csaq + Cssorb*MVw #(mol/L)
530 dt= 3600 # s
540 T_R= 40821.11 # = cel renew time (s)
550 nb_cel= dt*(step_no)/ T_R #nb renewed cel
560 T_exp = dt*(step_no) #experiment time (s)

570 Kini = sys("K") # (mol/L)
580 Naini = sys("Na") # (mol/L)
590 Caini = sys("Ca") # (mol/L)
700 Mgini = sys("Mg") # (mol/L)

800 punch nb_cel T_exp log10(Csaq) log10(Kd) surfCs exchCs surfK exchK surfNa exchNa surfCa exchCa surfMg exchMg Csverif Csini Kini Naini Caini Mgini Csaq Kaq Naaq Caaq Mgaq
 -end

TRANSPORT
-cells 1
-lengths 0.01 # 1 cm sol
-shifts 576 # number of time_steps total time is 24 days
-time_step 3600 # each time_step (s) the MIX is performed
-flow_direction  diffusion_only
-boundary_conditions constant closed
-diffusion_coefficient 2.3E-08
-stagnant 1
-punch_cells 3
-punch_frequency 1 # sample every time step
END 

I hope this will be clear enough I also attach the database If somebody wants to give it a try =)
Thank you for your help!

[dlparkhurst]

Looks like you scaled the EXCHANGE definition, but not the SURFACE definition in script 2.

Also, I think you need to be more careful in setting up the initial calculations. Separating the calculations with END statements will make your calculation sequence clearer. I think this is how script 1 actually begins.

Code: [Select]

SOLUTION 10
 units mol/L
pH 9.32
Na 5.00E-03
K 1.60E-03
Ca 4.50E-03
Mg 1.00E-03
N(5) 1.60E-03
Cl 5.00E-03
S(6) 1.00E-03
C(4) 4.50E-03
END
SURFACE 10
-equilibrate 10
-sites_units density
Illitex 2.07E-05 97 3695.29 # 0.5%
Montx 3.80E-07 800 3695.29 # 1.1%
Kaolinitex 1.46E-06 10 3695.29 # 0.7%
# ^site/nm^2 ^SSA (m^2/g) ^M/V (g/L)
END
EXCHANGE  10
-equilibrate 10
Illi 6.37E-03
Mont 3.39E-02
Kao 7.19E-04
END
USE solution 10
USE exchange 10
USE surface 10
EQUILIBRIUM_PHASES   10
Fix_H+ -9.32 HCl 10
save surface 3
save exchange 3
END


[dlparkhurst]

Note there are some keywords, documented in the RELEASE.txt file (http://water.usgs.gov/water-resources/software/PHREEQC/Phreeqc_ReleaseNotes.txt), that may be useful: MIX_EXCHANGE and MIX_SURFACE. You could use these to scale the moles of surface sites and exchangers.

These datablocks have not been used extensively, so reply if you have any issues.

[FlouretAlex]

Hello and thank you for your response!

Thank for the tip with the definition of the initial calculation it's indeed clearer like this.
For the problem of surface scaling, I think you didn't notice "-sites_units" was set as density and not absolute thus scaling was made by the change of mass. Thus to make it clearer you can find below the two scripts with "-site_units" set with absolute and see that the quantity of mol was scaled.

Script 1 with the mass set as if the soil solution was 1L ==> M/V = 18.31/4.94E-03 = 3695.29 g/L

Code: [Select]

ATABASE C:\Users\International\Documents\PhD\AlexF\Modelo alex\Fitted\000_Database_PHREEQC\PKDLM_BRGM_database_phreeqc_ThermoddemV1.10_06Jun2017.dat

SOLUTION 10
 units mol/L
pH 9.32
Na 5.00E-03
K 1.60E-03
Ca 4.50E-03
Mg 1.00E-03
N(5) 1.60E-03
Cl 5.00E-03
S(6) 1.00E-03
C(4) 4.50E-03
END

SURFACE 10
-equilibrate 10
-sites_units absolute
Illitex 1.23E-05 97 3695.29 # 0.5%
Montx 1.87E-06 800 3695.29 # 1.1%
Kaolinitex 8.96E-08 10 3695.29 # 0.7%
#^mol ^SSA (m^2/g) ^M/V (g/L)
END
EXCHANGE  10
-equilibrate 10
Illi 6.37E-03
Mont 3.39E-02
Kao 7.19E-04
END
# surface and exchange equilibration with natural soil water
USE solution 10
USE exchange 10
USE surface 10
EQUILIBRIUM_PHASES   10
Fix_H+ -9.32 HCl 10
save surface 3
save exchange 3
END

#Contaminated soil ==> RIZO
SOLUTION 3 # ==> stagnant cell
units mol/L
pH 9.32
Cl 1.00E-05
K 1.00E-05
Mg 7.50E-04
S(6) 7.50E-04
Cs 3.35E-04
use surface 3 #number of surface and exchange need to be the same as the cell here solution 3 ==> cell 3
use Exchange 3
EQUILIBRIUM_PHASES   3
Fix_H+ -9.32 HCl  10

END

SOLUTION 0-1 # inlet solution 0 # define solution 1 since PHREEQC wants at least 1 cell in the colum for transport ==> mobile cell.
 units mol/L
pH 9.32
Cl 1.00E-05
K 1.00E-05
Mg 7.50E-04
S(6) 7.50E-04

Solution 100 #for evaporation
END

MIX 3 # simulate % water evaporation of solution 3 then replaced by solution 0
0 0.09
100 -0.09
3 1
END

TRANSPORT
-cells 1
-lengths 0.01 # 1 cm sol
-shifts 576 # number of time_steps total time is 24 days
-time_step 3600 # each time_step (s) the MIX is performed
-flow_direction  diffusion_only
-boundary_conditions constant closed
-diffusion_coefficient 2.3E-08
-stagnant 1
-punch_cells 3
-punch_frequency 1 # sample every time step
END

Script 2 with mass and volume set with experimental data :
- M_soil =18.31 g.DW
- V_water = 4.94E-03 kg

Code: [Select]

DATABASE C:\Users\International\Documents\PhD\AlexF\Modelo alex\Fitted\000_Database_PHREEQC\PKDLM_BRGM_database_phreeqc_ThermoddemV1.10_06Jun2017.dat

 SOLUTION 10
-water 4.94E-03 #kg
 units mol/L
pH 9.32
Na 5.00E-03
K 1.60E-03
Ca 4.50E-03
Mg 1.00E-03
N(5) 1.60E-03
Cl 5.00E-03
S(6) 1.00E-03
C(4) 4.50E-03
END
 
SURFACE 10
-equilibrate 10
-sites_units absolute 
Illitex 6.10E-08 97 18.31
Montx 9.24E-09 800 18.31
Kaolinitex 4.44E-10 10 18.31
# ^mol ^SSA (m^2/g) ^M (g)
EXCHANGE  10
-equilibrate 10
Illi 3.16E-05
Mont 1.68E-04
Kao 3.56E-06
END
# surface and exchange equilibration with natural soil water
USE solution 10
USE exchange 10
USE surface 10
EQUILIBRIUM_PHASES   10
Fix_H+ -9.32 HCl 10
save surface 3
save exchange 3
END

# Equilibrage spike

EQUILIBRIUM_PHASES   10
Fix_H+ -9.32 HCl 10
save surface 3
save exchange 3
END

#Contaminated soil ==> RIZO
SOLUTION 3 # ==> stagnant cell
-water 4.94E-03 #kg
units mol/L
pH 9.32
Cl 1.00E-05
K 1.00E-05
Mg 7.50E-04
S(6) 7.50E-04
Cs 3.35E-04
use surface 3 #number of surface and exchange need to be the same as the cell here solution 3 ==> cell 3
use Exchange 3
EQUILIBRIUM_PHASES   3
Fix_H+ -9.32 HCl  10

END

SOLUTION 0 # inlet solution 0 # define solution 1 since PHREEQC wants at least 1 cell in the colum for transport ==> mobile cell.
units mol/L
pH 9.32
Cl 1.00E-05
K 1.00E-05
Mg 7.50E-04
S(6) 7.50E-04

Solution 1

Solution 100 #for evaporation
END

MIX 3 # simulate % water evaporation of solution 3 then replaced by solution 0
0 0.09
100 -0.09
3 1
END

TRANSPORT
-cells 1
-lengths 0.01 # 1 cm sol
-shifts 576 # number of time_steps total time is 24 days
-time_step 3600 # each time_step (s) the MIX is performed
-flow_direction  diffusion_only
-boundary_conditions constant closed
-diffusion_coefficient 2.3E-08
-stagnant 1
-punch_cells 3
-punch_frequency 1 # sample every time step
END
The only change I see between the two scripts is the quantity of mol from the solution 0 that is added in the mix between Solution 0; 100 and 3.

Indeed as shown in the previous post the molality of the différent element found in mixture 3 are different between the two scripts.

Mixture 3 from script 1

Code: [Select]

-----------------------------------------
Beginning of batch-reaction calculations.
-----------------------------------------

Reaction step 1.

Using mix 3.

Mixture 3.

  9.000e-02 Solution 0                                                        
  1.000e+00 Solution 3                                                        
-9.000e-02 Solution 100                                                        

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

Cl                1.090e-05   1.090e-05
Cs                3.350e-04   3.350e-04
K                 1.090e-05   1.090e-05
Mg                8.176e-04   8.176e-04
S                 8.176e-04   8.176e-04

Here we can see that mix did add 9% of KCl and MgSO4 from solution 0 and added to solution 3.

On the other hand with script 2 we have:

Code: [Select]

-----------------------------------------
Beginning of batch-reaction calculations.
-----------------------------------------

Reaction step 1.

Using mix 3.

Mixture 3.

  9.000e-02 Solution 0                                                        
  1.000e+00 Solution 3                                                        
-9.000e-02 Solution 100                                                        

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

Cl                1.922e-04   9.494e-07
Cs                3.350e-04   1.655e-06
K                 1.922e-04   9.494e-07
Mg                1.441e-02   7.121e-05
S                 1.441e-02   7.121e-05

Here mix also add 9% of solution 0 to solution 3 but the problem is that solution 3 does not have a volume of 1L thus the quantity of KCl and MgSO4 is not proportional in comparison to the first script (hence the change of morality for K, Cl, Mg and S. Thus, I think that the problem is more a problem with Mix than with the surface and exchange quantity.

But maybe I'm completely wrong ^^.

[dlparkhurst]

Sorry, looking more closely at your simulation, I have no idea what MIX does in a diffusion_only simulation. I am not good with MIX and TRANSPORT calculation in general.

If you are using diffusion only, I would simply evaporate from a cell with KINETICS and let diffusion move water and solutes.

I would also start with a very simple system, without EXCHANGE or SURFACE to see what is going on just with solutes.