Colloid transport

[Sen]

Dear Dr. Parkhurst
Here again for your help. Actinides are radioactive and toxic and difficult to obtain in the laboratory, and a large part of my work needs to be done with the computational capabilities of PHREEQC.

PHREEQC is able to calculate the migration of pseudocolloids https://www.hydrochemistry.eu/exmpls/colloid_U.phr. My target task is to calculate the long-distance migration of Am in groundwater over long periods of time. I have modified it from your code and redefined the thermodynamic data for Am and the hydrological data for the study area. Here is my code.

Code: [Select]

DATABASE D:\USGS\Phreeqc Interactive 3.4.0-12927\database\wateq4f.dat
SOLUTION_MASTER_SPECIES
         Am         Am+3   0.0   Am    243.0000
         Am(+3)       Am+3   0.0   Am         
SOLUTION_SPECIES
          Am+3 = Am+3
            log_k      0.0
         Am+3 +  H2O - H+  = AmOH+2
     log_k    -7.200
  Am+3 + 2H2O - 2H+  = Am(OH)2+
     log_k    -15.100
  Am+3 + 3H2O - 3H+  = Am(OH)3
     log_k     -26.200
         Am+3 + F-  = AmF+2
     log_k     3.400
           delta_h    2.89 kcal
  Am+3 + 2F-  = AmF2+
     log_k     5.800
           delta_h    10.77 kcal
  Am+3 +  Cl- = AmCl+2
     log_k     0.240
  Am+3 + 2Cl- = AmCl2+
     log_k     -0.810
            delta_h    13.11 kcal
  Am+3 + SO4-2 = AmSO4+
     log_k      3.500
            delta_h    9.56 kcal
  Am+3 + 2SO4-2 = Am(SO4)2-
     log_k      5.000
            delta_h    16.72 kcal
  Am+3 + NO3- = AmNO3+2
     log_k     1.28
            delta_h    0.43 kcal
  Am+3 + 2NO3- = Am(NO3)2+
     log_k     0.88
            delta_h    2.58 kcal
  Am+3 + CO3-2 = AmCO3+
     log_k     8.000
  Am+3 + 2CO3-2 = Am(CO3)2- 
     log_k     12.900
  Am+3 + 3CO3-2 = Am(CO3)3-3 
     log_k     15.000
         Am+3 + HCO3- = AmHCO3+2
     log_k      3.100
         Am+3 +  H3SiO4- = AmSiO(OH)3+2
            log_k      8.13
            delta_h    4.44 kcal

SURFACE_SPECIES
     Hfo_sOH  + H+ = Hfo_sOH2+
       log_k  7.18
     Hfo_sOH = Hfo_sO- + H+
       log_k  -8.82
     Hfo_wOH  + H+ = Hfo_wOH2+
       log_k  7.18
     Hfo_wOH = Hfo_wO- + H+
       log_k  -8.82
#Americium
     Hfo_sOH + Am+3 + 2H2O = Hfo_sOAm(OH)2 + 3H+ 
       log_k    -12.5                 
     Hfo_wOH + Am+3 + 2H2O = Hfo_wOAm(OH)2 + 3H+ 
       log_k   -13.65
     Hfo_sOH + Am+3 = Hfo_sOHAm+3
       log_k   10.97
     Hfo_wOH + Am+3 = Hfo_wOHAm+3     
       log_k    3.86
     Hfo_sOH + Am+3 + NO3- = Hfo_sOHAmNO3+2
       log_k      8.55
     Hfo_wOH + Am+3 + NO3- = Hfo_wOHAmNO3+2
       log_k      8.46
     Hfo_sOH + Am+3 + Cl- = Hfo_sOHAmCl+2
       log_k      8.92
     Hfo_wOH + Am+3 + Cl- = Hfo_wOHAmCl+2
       log_k      8.92
#Calcium
            Hfo_sOH + Ca+2 = Hfo_sOHCa+2
              log_k      4.97
            Hfo_wOH + Ca+2 = Hfo_wOCa+ + H+
              log_k     -5.85
   
SOLUTION_SPECIES
            Na+ = Na+;                 -log_k 0; -gamma 1e10 0
            Ca+2 = Ca+2;               -log_k 0; -gamma 1e10 0
            Cl- = Cl-;                 -log_k 0; -gamma 1e10 0
            H2O + 0.01e- = H2O-0.01;   -log_k -9

SURFACE_MASTER_SPECIES /* define sorbed form of Hfo = Sfo */
 Sfo_s Sfo_sOH; Sfo_w Sfo_wOH

SURFACE_SPECIES
           Sfo_sOH = Sfo_sOH
              log_k       0.0
           Sfo_sOH + H+ = Sfo_sOH2+
              log_k       7.29
           Sfo_sOH = Sfo_sO- + H+
              log_k      -8.93
           Sfo_wOH = Sfo_wOH
              log_k      0.0
           Sfo_wOH + H+ = Sfo_wOH2+
              log_k      7.29
           Sfo_wOH = Sfo_wO- + H+
              log_k      -8.93
# Calcium
           Sfo_sOH + Ca+2 = Sfo_sOHCa+2
               log_k     4.97
           Sfo_wOH + Ca+2 = Sfo_wOCa+ + H+
               log_k     -5.85
#Americium
           Sfo_sOH + Am+3 + 2H2O = Sfo_sOAm(OH)2 + 3H+ 
        log_k      -12.5                 
           Sfo_wOH + Am+3 + 2H2O = Sfo_wOAm(OH)2 + 3H+ 
        log_k      -13.65
           Sfo_sOH + Am+3 = Sfo_sOHAm+3
                log_k     10.97
           Sfo_wOH + Am+3 = Sfo_wOHAm+3
                log_k     3.86

SOLUTION 1-21   #Repository groundwater
        temp        25
        pH          8.9
        pe          -3
        redox       pe   
        units       mol/L
        density     1.0   
        Na          3.23e-2
        Ca          3.42e-4
        K           1.16e-4
        Mg          1.19e-4
        Cl          1.79e-2     charge
        S(6)        1.15e-4
        N(5)        3.67e-4     gfw   62.0
        C(4)        2.31e-3     as    HCO3 
        Si          7.05e-5

EQUILIBRIUM_PHASES 1-21
       CO2(g)   -5.5     
SURFACE 1-21            # define small conc's
 Hfo_w 97.5e-15 600 88e-13 Dw 1e-13
 Hfo_s 2.5e-15
 Sfo_w 97.5e-15 600 88e-13 Dw 0
 Sfo_s 2.5e-15
 -donnan 1e-10
 -equil 1
END   

SOLUTION 0  Palse solution with Am
        temp        25
        pH          8.9
        pe          -3
        units       mol/L
        Ca          5.0e-2 
        Am          1e-7
   

SURFACE 0
 Hfo_w 97.5e-5 600 88e-3 Dw 1e-13
 Hfo_s 2.5e-5
 -donnan 1e-10 #d 1 l 0.99 v 1# 1e-9
 -equil 0
END
PRINT; -reset false; -status false

RATES
Sorb_hfo
 -start
 10 if tot("Ca") > 1e-10 then sar = tot("Na") * 10^1.5 / tot("Ca")^0.5 else sar = 1
 20 if sar < 11 then rate = -40 * tot("Hfo_w") else rate = 400 * tot("Sfo_w")
 30 save rate * time
 40 put(sar, 99)
 -end
KINETICS 1-21
 Sorb_Hfo; -formula Hfo_w 0.975 Hfo_s 0.025 Sfo_w -0.975 Sfo_s -0.025

SELECTED_OUTPUT
   -file           Colloid_Am-transport.sel
   -reset          false
   -solution
   -distance       true
   -time           true
   -step           true
   -pH             true
   -totals         Am
   -molalities     AmSiO(OH)3+2    AmOH+2   Am(OH)2+    AmCO3+    Am(CO3)2-   Am(CO3)3-3   
                   Sfo_sOAm(OH)2   Sfo_wOAm(OH)2     Sfo_sOHAm+3   Sfo_wOHAm+3
                   Hfo_sOAm(OH)2   Hfo_wOAm(OH)2     Hfo_sOHAm+3   Hfo_wOHAm+3

USER_GRAPH
 -head SAR Hfo Am_Hfo Am_tot
 -axis_titles "Pore Volumes" "mmol / L" "Sodium Adsorption Ratio (SAR)"
 -chart_title "Leaching Am with Ferrihydrite"
 -axis_scale x_axis 0  auto
 -axis_scale y_axis 0  auto
 -plot_concentration_vs X

 -start
 10 PV = (step_no + 0.5) / 1 / cell_no
 12 if sim_no > 4 then PV = PV + 1.0
 20 Am_sol = 1e3 * tot("Am")
 30 Am_col = 1e3 * (mol("Sfo_sOAm(OH)2") + mol("Sfo_wOAm(OH)2")+ mol("Sfo_sOHAm+3") + mol("Sfo_wOHAm+3")) 
 40 plot_xy PV, get(99), color = Green, symbol = None, y_axis = 2
 50 plot_xy PV, (tot("Hfo_w") + tot("Hfo_s")) * 1e3, color = Orange, symbol = Circle, symbol_size = 10, line_width = 0
 60 plot_xy PV, Am_col, color = Magenta, symbol = Circle, symbol_size = 6, line_width = 0
 70 plot_xy PV, Am_sol + Am_col, color = Magenta, symbol = None
 -end



TRANSPORT
-cells                20
-length               30        # total  length = length *  cells  # total length is 600m
-shifts               20  1         # total diffusion time is  shifts*time_step
-time                 3.1536e+7       # second   
-dispersivity         0.078           

-punch_cells          20
-punch_frequency      1-20
-print                20 
-print_frequency      10   
-flow_direction       forward    # diffusion
-boundary_cond        flux    flux
-correct_disp         true
-diffc                2.0e-13    #diffusion coefficient
-multi_d              true  1e-9   5.8e-3   0    1
-porosities           20*5.8e-3
END

PRINT; -user_graph false

SURFACE 0
SOLUTION 0  Palse solution with Am
        temp        25
        pH          8.9
        pe          -3
        redox       pe   
        units       mol/L
        density     1.0   
        Na          3.23e-2
        Ca          3.42e-4
        K           1.16e-4
        Mg          1.19e-4
        Cl          1.79e-2     
        S(6)        1.15e-4
        N(5)        3.67e-4     gfw   62.0
        C(4)        2.31e-3     as    HCO3 
        Si          7.05e-5
EQUILIBRIUM_PHASES 0
       CO2(g)   -5.5
END

PRINT; -user_graph true
TRANSPORT
 -shifts  50 1
END
where KINETICS does not make any changes.

But the reworked code does not converge well. I guess there might be a problem with the reworked Transport module. I reworked the model from -time 3.1536e+7 to 3.1536e-4 and it converged. please help me to see what is wrong with the model?

Also, I have been studying USER_GRAPH for a long time and I still find it difficult to learn. Can you help me to change the graph to concentration vs Dis-X?

Thank you for your continuous selfless help!

Best wishes,
Sen

[dlparkhurst]

Tony adjusted the rate of transformation between mobile and immobile pools to occur over the length (1 m) and time (0.05 s) of his calculation. By increasing the time scale, the reaction (Sfo <--> Hfo) occurs over a much longer time. The long times allow Hfo or Sfo to be completely consumed, which causes the rate integration to have issues. I think the rate equation can require more than all of the available Hfo to be removed and the integration cannot get past that point in time. I have modified the rate equation so that the rate gets small when the concentrations of Hfo or Sfo get small. 

The calculation is still slow, so I only used 5 cells. You can expand back to 10 or 20 when you are happy with the way the reactions are working.

Colloids are entering the column because a diffusion coefficient has been defined for Hfo in the SURFACE definition.

The RATES definition has a switch so that when SAR is less than 11, the effect is to decrease colloidal Hfo and increase non-mobile Sfo. When SAR is greater than 11, non-mobile Sfo is transformed into mobile Hfo. Because you have no Na in your Am solution and the time step is long, as the Am solution and SURFACE colloid enter the system, all of the colloid accumulates in the sediment. USER_GRAPH 1 shows that all of the Am is accumulated in the first cell after 1 pore volume.

For the next TRANSPORT calculation, no colloids and no Am are entering the system. The SAR is large, which releases all of the Am from the sediment, and a pulse of dissolved Am moves through the system. After 3 years, the pulse has moved more than halfway through the column as shown in USER_GRAPH 2.

The colloid-sediment transition will be spread out farther in the column if you go back to Tony's time steps, or you slow down the reaction defined in RATES. You will have to decide if the SAR approach is appropriate for Americium.

Code: [Select]

#DATABASE D:\USGS\Phreeqc Interactive 3.4.0-12927\database\wateq4f.dat
SOLUTION_SPECIES
 Na+ = Na+; -log_k 0; -gamma 1e10 0
 Ca+2 = Ca+2; -log_k 0; -gamma 1e10 0
 Cl- = Cl-; -log_k 0; -gamma 1e10 0
 H2O + 0.01e- = H2O-0.01; -log_k -9
SOLUTION_MASTER_SPECIES
         Am         Am+3   0.0   Am    243.0000
         Am(+3)       Am+3   0.0   Am         
SOLUTION_SPECIES
          Am+3 = Am+3
            log_k      0.0
         Am+3 +  H2O - H+  = AmOH+2
     log_k    -7.200
  Am+3 + 2H2O - 2H+  = Am(OH)2+
     log_k    -15.100
  Am+3 + 3H2O - 3H+  = Am(OH)3
     log_k     -26.200
         Am+3 + F-  = AmF+2
     log_k     3.400
           delta_h    2.89 kcal
  Am+3 + 2F-  = AmF2+
     log_k     5.800
           delta_h    10.77 kcal
  Am+3 +  Cl- = AmCl+2
     log_k     0.240
  Am+3 + 2Cl- = AmCl2+
     log_k     -0.810
            delta_h    13.11 kcal
  Am+3 + SO4-2 = AmSO4+
     log_k      3.500
            delta_h    9.56 kcal
  Am+3 + 2SO4-2 = Am(SO4)2-
     log_k      5.000
            delta_h    16.72 kcal
  Am+3 + NO3- = AmNO3+2
     log_k     1.28
            delta_h    0.43 kcal
  Am+3 + 2NO3- = Am(NO3)2+
     log_k     0.88
            delta_h    2.58 kcal
  Am+3 + CO3-2 = AmCO3+
     log_k     8.000
  Am+3 + 2CO3-2 = Am(CO3)2-
     log_k     12.900
  Am+3 + 3CO3-2 = Am(CO3)3-3
     log_k     15.000
         Am+3 + HCO3- = AmHCO3+2
     log_k      3.100
         Am+3 +  H3SiO4- = AmSiO(OH)3+2
            log_k      8.13
            delta_h    4.44 kcal
END
SURFACE_SPECIES
     Hfo_sOH  + H+ = Hfo_sOH2+
       log_k  7.18
     Hfo_sOH = Hfo_sO- + H+
       log_k  -8.82
     Hfo_wOH  + H+ = Hfo_wOH2+
       log_k  7.18
     Hfo_wOH = Hfo_wO- + H+
       log_k  -8.82
#Americium
     Hfo_sOH + Am+3 + 2H2O = Hfo_sOAm(OH)2 + 3H+
       log_k    -12.5                 
     Hfo_wOH + Am+3 + 2H2O = Hfo_wOAm(OH)2 + 3H+
       log_k   -13.65
     Hfo_sOH + Am+3 = Hfo_sOHAm+3
       log_k   10.97
     Hfo_wOH + Am+3 = Hfo_wOHAm+3     
       log_k    3.86
     Hfo_sOH + Am+3 + NO3- = Hfo_sOHAmNO3+2
       log_k      8.55
     Hfo_wOH + Am+3 + NO3- = Hfo_wOHAmNO3+2
       log_k      8.46
     Hfo_sOH + Am+3 + Cl- = Hfo_sOHAmCl+2
       log_k      8.92
     Hfo_wOH + Am+3 + Cl- = Hfo_wOHAmCl+2
       log_k      8.92
#Calcium
            Hfo_sOH + Ca+2 = Hfo_sOHCa+2
              log_k      4.97
            Hfo_wOH + Ca+2 = Hfo_wOCa+ + H+
              log_k     -5.85
END
SURFACE_MASTER_SPECIES /* define sorbed form of Hfo = Sfo */
 Sfo_s Sfo_sOH; Sfo_w Sfo_wOH

SURFACE_SPECIES
           Sfo_sOH = Sfo_sOH
              log_k       0.0
           Sfo_sOH + H+ = Sfo_sOH2+
              log_k       7.29
           Sfo_sOH = Sfo_sO- + H+
              log_k      -8.93
           Sfo_wOH = Sfo_wOH
              log_k      0.0
           Sfo_wOH + H+ = Sfo_wOH2+
              log_k      7.29
           Sfo_wOH = Sfo_wO- + H+
              log_k      -8.93
# Calcium
           Sfo_sOH + Ca+2 = Sfo_sOHCa+2
               log_k     4.97
           Sfo_wOH + Ca+2 = Sfo_wOCa+ + H+
               log_k     -5.85
#Americium
           Sfo_sOH + Am+3 + 2H2O = Sfo_sOAm(OH)2 + 3H+
        log_k      -12.5                 
           Sfo_wOH + Am+3 + 2H2O = Sfo_wOAm(OH)2 + 3H+
        log_k      -13.65
           Sfo_sOH + Am+3 = Sfo_sOHAm+3
                log_k     10.97
           Sfo_wOH + Am+3 = Sfo_wOHAm+3
                log_k     3.86
END
SOLUTION 1-21   #Repository groundwater
        temp        25
        pH          8.9
        pe          -3
        redox       pe   
        units       mol/L
        density     1.0   
        Na          3.23e-2
        Ca          3.42e-4
        K           1.16e-4
        Mg          1.19e-4
        Cl          1.79e-2     charge
        S(6)        1.15e-4
        N(5)        3.67e-4     gfw   62.0
        C(4)        2.31e-3     as    HCO3
        Si          7.05e-5
END
SURFACE 1-21                        # define small conc's
 Hfo_w 97.5e-15 600 88e-13 Dw 1e-13
 Hfo_s 2.5e-15
 Sfo_w 97.5e-15 600 88e-13 Dw 0
 Sfo_s 2.5e-15
 -donnan 1e-10
 -equil 1
END
EQUILIBRIUM_PHASES 1-21
       CO2(g)   -5.5     
END   
SOLUTION 0  Palse solution with Am
        temp        25
        pH          8.9
        pe          -3
        units       mol/L
        Ca          5.0e-2
        Am          1e-7
END
SURFACE 0
 Hfo_w 97.5e-5 600 88e-3 Dw 1e-13
 Hfo_s 2.5e-5
 -donnan 1e-10 #d 1 l 0.99 v 1# 1e-9
 -equil 0
END
#PRINT; -reset false; -status false

RATES
Sorb_hfo
 -start
 10 if tot("Ca") > 1e-10 then sar = tot("Na") * 10^1.5 / tot("Ca")^0.5 else sar = 1
 12 k1 = -40e0
 14 k2 = 400e0
 20 if sar < 11 then rate = k1 * tot("Hfo_w") else rate = k2 * tot("Sfo_w")
 22 fSw = tot("Sfo_w") / (1e-18 + tot("Sfo_w"))
 24 fSs = tot("Sfo_s") / (1e-18 + tot("Sfo_s"))
 26 fHw = tot("Hfo_w") / (1e-10 + tot("Hfo_w"))
 28 fHs = tot("Hfo_s") / (1e-10 + tot("Hfo_s"))
 30 moles = rate * time
 32 if (moles < 0) then moles = moles * fHw * fHs
 34 if (moles > 0) then moles = moles * fSw * fSs
 36 SAVE moles
 40 put(sar, 99)
 -end

KINETICS 1-21
 Sorb_Hfo; -formula Hfo_w 0.975 Hfo_s 0.025 Sfo_w -0.975 Sfo_s -0.025
#-cvode
END

TRANSPORT
    -cells                 5
    -shifts                5
    -time_step             31536000 1 # seconds
    -lengths               30
    -dispersivities        5*0.078
    -correct_disp          true
    -diffusion_coefficient 2e-13
    -print_cells           1-5
    -print_frequency       1
    -punch_cells           1-5
    -punch_frequency       5
    -multi_d               true 1e-09 0.0058 0 1

USER_GRAPH 1
    -headings               dist Am_col Am_dis Am_sediment
    -axis_titles            "Distance, m" "mmol / L" "Sodium Adsorption Ratio (SAR)"
    -chart_title            "One pore volume, 5 years"
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
 20 Am_sol = 1e3 * tot("Am")
 30 Am_col = 1e3 * SURF("Am", "Hfo")
 40 GRAPH_X DIST
 60 GRAPH_Y Am_col, Am_sol
 70 GRAPH_Y SURF("Am", "Sfo") * 1e3
-end
    -active                 true
END
USER_GRAPH 1
-detach
END
SURFACE 0
SOLUTION 0  Palse solution with Am
        temp        25
        pH          8.9
        pe          -3
        redox       pe   
        units       mol/L
        density     1.0   
        Na          3.23e-2
        Ca          3.42e-4
        K           1.16e-4
        Mg          1.19e-4
        Cl          1.79e-2     
        S(6)        1.15e-4
        N(5)        3.67e-4     gfw   62.0
        C(4)        2.31e-3     as    HCO3
        Si          7.05e-5
EQUILIBRIUM_PHASES 0
       CO2(g)   -5.5
END

TRANSPORT
 -shifts  3 1
 -punch_frequency 3
USER_GRAPH 2
    -headings               Dist Am_col_5y Am_dis_5y Am_sediment_5y
    -axis_titles            "Distance, m" "mmol / L" "Sodium Adsorption Ratio (SAR)"
    -chart_title            "1.6 Pore volumes, 8 years"
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
 20 Am_sol = 1e3 * tot("Am")
 30 Am_col = 1e3 * SURF("Am", "Hfo")
 40 GRAPH_X DIST
 60 GRAPH_Y Am_col, Am_sol
 70 GRAPH_Y SURF("Am", "Sfo") * 1e3
-end
    -active                 true
END


[Sen]

Dear Dr. Parkhurst
After adjusting, I think I got a curve that matches better with the literature. But what does dist Am_col , Am_dis , Am_sediment stand for respectively, can you help me to annotate it?
Once again, I would like to express my gratitude for your help and dedication all along.

best wishes，
Sen

[dlparkhurst]

Look at the USER_GRAPH definition and figure out what is being plotted.