pH charge balance during transport with different surface amounts

[Tim]

Dear phreeqc users,

I would appreciate your thoughts on the following observation(s). I've got a diffusion only transport model with 1 M NaCl and a pH of 8.2 in the 10 model cells and 0.1 M in the boundary cell. I'm initially equilibrating the model cells with a surface assemblage varying the amount of surface sites: (1) a high surface site density (2.26 sites/nm?) and (2) a low surface site density (0.0001 sites/nm?) using the donnan approach and ionic-strength dependent calculation of the donnan layer thickness.
What I don't fully understand is what's happening with the pH in the solution in those two cases after the first (and the following) transport steps. It increases on the model boundary in case 1 and decreases in case 2.

I understand that Phreeqc charge balances the solution in transport simulations by using the pH. My first question is: is there a way to easily extract if phreeqc actually did charge balance the pH in a transport step and by which degree or if the pH changed as a result of H+ movement between the diffuse layer and the solution?

Besides this more general question, there is a number of observations I have a hard time to reconcile. In both cases, ionic strength is decreasing in the model continuously, the same amount of moles of Na+ and Cl- is leaving the model towards the boundary - makes sense to me. But Na+ and Cl- are not equally reduced in the diffuse layer in the two cases. Slightly more Cl- than Na+ is leaving the DL in case 1 (many surfaces) and vice versa in case 2. Right now, I guess this could be the reason for the pH to charge balance in the one direction in case 1 (needing to reduce H+ to compensate slight surplus of Na+ in the DL) and the other in case 2 (needs to add H+ to compensate the slight surplus of Cl- in the DL).

However, from the manual (version 3, page 332) it reads "If the identifier -Donnan is added in keyword SURFACE, the composition of the electrical double layer is calculated and stored with each of the surfaces in the assemblage; this explicitly calculated electrical double layer makes each surface a neutral entity." - but if it would be a neutral entity, there would not be the need for a charge balance in the solution, right?

I noticed that the surface charge after the first transport step does not change in case 1, but does for case 2 because there is some protonation happening on the surfaces, therefore the negative charge decreases (but what triggers the protonation in this case?). Ultimately, I can't put together a coherent overall picture from these pieces of observations, so I would be very glad if someone could shed some light on the matter. Please find my piece of code below.

Best,
Tim

Code: [Select]

DATABASE ~/phreeqc.dat

SURFACE_MASTER_SPECIES
        Ill_s           Ill_sOH     

SURFACE_SPECIES               
### strong-binding sites ###
       
   Ill_sOH = Ill_sOH
   log_K   0.0
       Ill_sOH + H+ = Ill_sOH2+
        log_K = 4.9                                           
   Ill_sOH = Ill_sO- + H+
        log_K = -6.8                                         

SOLUTION 1-10 initial composition - higher ionic strength
     pe           6           
     density      1           
     temp         13             
     units        mol/kgw
     pH           8.2     
     Na           1
     Cl           1
      -water       1

SAVE SOLUTION 1-10
END

USE SOLUTION 1
SURFACE 1
 -equil 1               
 -sites_units density
 -donnan debye_lengths 1 limit_ddl 0.6   
# surface  density SSA (m2/g)  mass (g/kgw)         
  Ill_sOH  2.26    100.        3950.30       # case 1: high site density - many surface sites
#  Ill_sOH  0.0001    100.        3950.30    # case 2: low site density - few surface sites

COPY SURFACE 1 2-10
SAVE SURFACE 1-10
END

SOLUTION 11 Boundary condition - lower ionic strength
     pe           6             
     density      1         
     temp         13             
     units        mol/kgw
     pH           8.2 
     Na           0.1
     Cl           0.1
     -water       1
END 

KNOBS
 -logfile            false
 -iterations         10000     
 -convergence_tolerance  1E-12
 -step_size          2       
 -pe_step_size       2   

TRANSPORT
 -cells               10
 -shifts              1000           
 -time_step           4733640000         # = 150 a 
 -flow_direction      diffusion_only
 -boundary_conditions 1 1          # 1=Dirichlet
 -lengths             2
 -diffusion_coefficient 5.32e-11   # diffusion coefficient
 -punch_cells         1-10
 -punch_frequency     1
 -print_cells         9 10
 -print_frequency     1

END

[dlparkhurst]

Sorry to take a while to respond. I spent some time experimenting and trying to decide how to reply.

I think you should probably use -multi_D for multicomponent diffusion if you want to look at transport processes in detail. Allowing ions to diffuse individually (subject to charge balance) should be more realistic. In addition, there are functions MCD_JCONC and MCD_JTOT for the fluxes of ions in each cell, which makes some processes easier to visualize.

My approach would be to ignore the surface, at least initially. The same steady state in solution should evolve whether the surfaces are present or not. If you find the steady state concentrations, then the surfaces will ultimately be in equilibrium with the solution concentrations. I leave it to you to consider the evolution of the initial condition to the steady state.

It seemed like there were some very small fluxes that required very long times to equilibrate H+/OH- at a pH of 8, so I imposed a gradient from pH 3 to pH 4 to speed up the calculations. Also, using the Debye length made it harder to find steady state, so I used the default thickness of 1e-8 m.

Here is a calculation that produces steady state without the surface, and then adds the surface. That should give you the end point of where concentrations and surface compositions are headed. You can look at the effect of surface sites and the time evolution.

Code: [Select]

#DATABASE ~/phreeqc.dat
SURFACE_MASTER_SPECIES
        Ill_s           Ill_sOH     
SURFACE_SPECIES               
### strong-binding sites ###
       
   Ill_sOH = Ill_sOH
        log_K   0.0
   Ill_sOH + H+ = Ill_sOH2+
        log_K = 4.9                                           
   Ill_sOH = Ill_sO- + H+
        log_K = -6.8                                         
END
SOLUTION 0-10 initial composition - higher ionic strength
     pe           6           
     density      1           
     temp         13             
     units        mol/kgw
     pH           3 #8.2     
     Na           1 charge
     Cl           1
      -water       1
END
SOLUTION 11 Boundary condition - lower ionic strength
     pe           6             
     density      1         
     temp         13             
     units        mol/kgw
     pH           4 #8.2
     Na           0.1 charge
     Cl           0.1
     -water       1
END
TRANSPORT
    -cells                 10
    -shifts                200
    -time_step             3.15e9 # 100 years
    -flow_direction        diffusion_only
    -boundary_conditions   constant constant
    -lengths               10*2
    -diffusion_coefficient 5.32e-11
    -thermal_diffusion     4   5.32e-11
    -print_cells           10
    -print_frequency       10
    -punch_cells           1-10
    -punch_frequency       1
    -multi_d               true 5.32e-11 1 0 0
    #-implicit true
USER_GRAPH 1
    -headings               dist H+ OH- Na+ Cl-
    -axis_titles            "Time, years" "Flux H+ OH-" "Flux Na+ and Cl-"
    -chart_title            "Cell 1"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
  5 if cell_no <> 1 then goto 100
 10 GRAPH_X total_time / 3.15e7
 30 GRAPH_Y MCD_JTOT("H+"), MCD_JTOT("OH-")
 40 GRAPH_SY MCD_JTOT("Na+"), MCD_JTOT("Cl-")
100 end
  -end
    -active                 true

USER_GRAPH 2
    -headings               dist H+ OH- Na+ Cl-
    -axis_titles            "Time, years" "Flux H+ OH-" "Flux Na+ and Cl-"
    -chart_title            "Cell 10"
    -axis_scale sy_axis     auto auto auto auto #log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
5 if cell_no <> 10 then goto 100
10 GRAPH_X total_time / 3.15e7
30 GRAPH_Y MCD_JTOT("H+"), MCD_JTOT("OH-")
40 GRAPH_SY MCD_JTOT("Na+"), MCD_JTOT("Cl-")
100 end
  -end

USER_GRAPH 3
    -headings               dist pH LA(Na+) LA(Cl-) LA(H+) LA(OH-)
    -axis_titles            "Distance" "pH"
    -chart_title            "All cells"
    -axis_scale y_axis      2.5 4.5 auto auto
    -axis_scale sy_axis     auto auto auto auto #log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
5 if ABS(step_no mod 200) > 1e-6 then goto 100
10 GRAPH_X dist
30 GRAPH_Y -LA("H+")
40 GRAPH_SY LA("Na+"), LA("Cl-"), LA("H+"), LA("OH-")
100 end
  -end

USER_GRAPH 4
    -headings               dist H+ OH- Na+ Cl-
    -axis_titles            "Distance" "Flux H+ OH-" "Flux Na+ and Cl-"
    -chart_title            "All cells"
    -axis_scale sy_axis     auto auto auto auto #log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
5 if ABS(step_no mod 200) > 1e-6 then goto 100
10 GRAPH_X dist
30 GRAPH_Y MCD_JTOT("H+"), MCD_JTOT("OH-")
40 GRAPH_SY MCD_JTOT("Na+"), MCD_JTOT("Cl-")
100 end
  -end
END
SURFACE 1
 -equil 1               
 -sites_units density
 #-donnandebye_lengths 1 limit_ddl 0.6   
# surface  density SSA (m2/g)  mass (g/kgw)         
  Ill_sOH  2.26    100      3950.30       # case 1: high site density - many surface sites
 # Ill_sOH  0.0001    100       3950.30    # case 2: low site density - few surface sites
END
SURFACE 2
 -equil 2               
 -sites_units density
 #-donnandebye_lengths 1 limit_ddl 0.6   
# surface  density SSA (m2/g)  mass (g/kgw)         
  Ill_sOH  2.26    100      3950.30       # case 1: high site density - many surface sites
 # Ill_sOH  0.0001    100       3950.30    # case 2: low site density - few surface sites
END
SURFACE 3
 -equil 3             
 -sites_units density
 #-donnandebye_lengths 1 limit_ddl 0.6   
# surface  density SSA (m2/g)  mass (g/kgw)         
  Ill_sOH  2.26    100      3950.30       # case 1: high site density - many surface sites
 # Ill_sOH  0.0001    100       3950.30    # case 2: low site density - few surface sites
END
SURFACE 4
 -equil 4             
 -sites_units density
 #-donnandebye_lengths 1 limit_ddl 0.6   
# surface  density SSA (m2/g)  mass (g/kgw)         
  Ill_sOH  2.26    100      3950.30       # case 1: high site density - many surface sites
 # Ill_sOH  0.0001    100       3950.30    # case 2: low site density - few surface sites
END
SURFACE 5
 -equil 5
 -sites_units density
 #-donnandebye_lengths 1 limit_ddl 0.6   
# surface  density SSA (m2/g)  mass (g/kgw)         
  Ill_sOH  2.26    100      3950.30       # case 1: high site density - many surface sites
 # Ill_sOH  0.0001    100       3950.30    # case 2: low site density - few surface sites
END
SURFACE 6
 -equil 6               
 -sites_units density
 #-donnandebye_lengths 1 limit_ddl 0.6   
# surface  density SSA (m2/g)  mass (g/kgw)         
  Ill_sOH  2.26    100      3950.30       # case 1: high site density - many surface sites
 # Ill_sOH  0.0001    100       3950.30    # case 2: low site density - few surface sites
END
SURFACE 7
 -equil 7             
 -sites_units density
 #-donnandebye_lengths 1 limit_ddl 0.6   
# surface  density SSA (m2/g)  mass (g/kgw)         
  Ill_sOH  2.26    100      3950.30       # case 1: high site density - many surface sites
 # Ill_sOH  0.0001    100       3950.30    # case 2: low site density - few surface sites
END
SURFACE 8
 -equil 8             
 -sites_units density
 #-donnandebye_lengths 1 limit_ddl 0.6   
# surface  density SSA (m2/g)  mass (g/kgw)         
  Ill_sOH  2.26    100      3950.30       # case 1: high site density - many surface sites
 # Ill_sOH  0.0001    100       3950.30    # case 2: low site density - few surface sites
END
SURFACE 9
 -equil 9             
 -sites_units density
 #-donnandebye_lengths 1 limit_ddl 0.6   
# surface  density SSA (m2/g)  mass (g/kgw)         
  Ill_sOH  2.26    100      3950.30       # case 1: high site density - many surface sites
 # Ill_sOH  0.0001    100       3950.30    # case 2: low site density - few surface sites
END
SURFACE 10
 -equil 10             
 -sites_units density
 #-donnandebye_lengths 1 limit_ddl 0.6   
# surface  density SSA (m2/g)  mass (g/kgw)         
  Ill_sOH  2.26    100      3950.30       # case 1: high site density - many surface sites
 # Ill_sOH  0.0001    100       3950.30    # case 2: low site density - few surface sites
END
TRANSPORT
    -cells                 10
    -shifts                200
USER_GRAPH 3
    -headings               dist pH LA(Na+)_surf LA(Cl-)_surf LA(H+)_surf LA(OH-)_surf

USER_GRAPH 4
    -headings               dist H+_surf OH-_surf Na+_surf Cl-_surf

USER_GRAPH 5
    -headings               dist pH Ill_sOH Ill_sOH2+ Ill_sO-
    -axis_titles            "Time, years" "pH" "Moles"
    -chart_title            "Cell 10"
    -axis_scale y_axis      2.5 4.5 auto auto
    -axis_scale sy_axis     auto auto auto auto #log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
5 if cell_no <> 10 then goto 100
10 GRAPH_X total_time / 3.15e7
20 GRAPH_Y -LA("H+")
30 GRAPH_SY MOL("Ill_sOH")*TOT("water"), MOL("Ill_sO-")*TOT("water"), MOL("Ill_sOH2+")*TOT("water")

100 end
  -end


[dlparkhurst]

And, as for charge balance, it is only a convenience to associate pH with charge balance. In fact all equations are simultaneous, and all unknowns participate in the solution, and charge balance is required. You can think of the numerical solution as balancing a chemical equation. To balance the simplest equation, it is necessary to conserve H, O, and charge, for example.

Code: [Select]

H2O = H2 + 1/2O2
or
H2O = H+ + OH-