Reactive transport with sorption: convergence problem

[llanduyt]

Hello,

In the framework of my master thesis I am trying to model the effect of pH and pe on the retardation of Sb by Hfo in river systems. In this perspective, I have set up a generic scenario in which I vary pH and pe.

I am experiencing convergence problems for low pH values and have not been able to find a solution yet.. I have been verifying amongst others the Courant condition but since I am not modelling diffusion and I calculated time_step based on the length of the cells and the velocity I wanted to imply, I think this condition is satisfied. Are there other conditions that should be taken into account? For example, I am modelling 1m x 1m x 250 m cells but have not explicitly specified the amount of water. Is this necessary? Doing this does not seem to solve my convergence problem..

When I decrease the Hfo concentration or increase the pH, no convergence problems occur. This seems to imply the convergence problem is due to a (too) high retardation, but I have not been able to find a practical solution..

Thanks in advance for any help or tips!

Greetings,

Lisa

My script (in the adapted database, I deleted all the speciation reactions of Sb so that only Sb(OH)3 and Sb(OH)6- occur):
I want to imply a velocity of 1.2 m/s, which leads to a time_step of 208.33 s for cells with a length of 250 m.

database C:\phreeqc\database\adapted_minteq.dat
   
PHASES
Fix_H+
  H+ = H+
  log_k 0.
Fix_pe
  e- = e-
  log_k 0.
 
SOLUTION 0
pH 3
Na 10
Cl 10
Sb 1E-03 # 10^-3 mmol = 1 µmol
EQUILIBRIUM_PHASES
Fix_H+ -3 NaOH 
   -force_equality true
Fix_pe -2 O2(g)   10
   -force_equality true
END

SOLUTION 1-40
pH 3
Na 10
Cl 10
EQUILIBRIUM_PHASES
Fix_H+ -3 NaOH   
   -force_equality true
Fix_pe -2 O2(g)   10
   -force_equality true
SURFACE 1-40
Hfo_s      5.9533E-05   600 1.057901145
Hfo_w      0.002381319
end

TRANSPORT
 -cells       40
 -shifts       5000
 -time_step      208.33
 -flow_direction     forward
 -boundary_conditions    flux
 -lengths       250
 -dispersivities    0
 -diffusion_coefficient 0
 -stagnant       0
 -print_cells      10
 -print_frequency   10
 -punch_cells      10
 -punch_frequency   10
 
 SELECTED_OUTPUT
 -file advectiontest_sorb.sel
 -reset false
 -step
 -high_precision true

 USER_PUNCH
 -headings timing Sb_diss Sb_sorb Sb_sorb2 Sb_tot pe pH
 -start
10 microSb=TOT("Sb")*1e+06
20 sorbedSb=(MOL("Hfo_sSbO(OH)4")+MOL("Hfo_sOHSbO(OH)4-")+MOL("Hfo_wSbO(OH)4")+MOL("Hfo_wOHSbO(OH)4-")+MOL("Hfo_sSbO(OH)2")+MOL("Hfo_wSbO(OH)2"))*1e+06
30 sorbedSb2 = SURF("Sb", "Hfo")*1e+06
40 totSb=SYS("Sb")*1e+06
50 PUNCH total_time/3600
60 PUNCH microSb
70 PUNCH sorbedSb, sorbedSb2
80 PUNCH totSb
90 PUNCH -la("e-")
100 PUNCH -la("H+")
-end
 
END

[dlparkhurst]

I just ran your input file with minteq.dat without a problem. Perhaps I need your database file, or there is some modification that is needed to demonstrate the failure.

I do not think the following has the effect you expected. There is an initial solution calculation with pH 3 and pe 4; then a reaction calculation that adjusts the pe to -2 and keeps the pH 3. However, because there is no SAVE statement the reaction calculation is not preserved.

SOLUTION 0
pH 3
Na 10
Cl 10
Sb 1E-03 # 10^-3 mmol = 1 µmol
EQUILIBRIUM_PHASES
Fix_H+ -3 NaOH 
   -force_equality true
Fix_pe -2 O2(g)   10
   -force_equality true
END

I would probably use the following, which will be similar, but conceptually, not quite the same.

SOLUTION 0
pe -2
pH 3
Na 10 charge
Cl 10
Sb 1E-03 # 10^-3 mmol = 1 µmol

[llanduyt]

Dear Mr. Parkhurst,

First of all, thank you very much for the quick reply.

The database file was too large to attach it to the post, but I now added it in a zipped format. The most important changes I made to the database are the following:

Sb(OH)6- + 2e- + 3H+ = Sb(OH)3 + 3H2O
   log_k   25.15 # originally 24.31
   delta_h   0   kJ
   -gamma   0   0
      #                  Id:   7407410
      #        log K source:   Bard85                       
      #      Delta H source:   MTQ3.11                       
      #T and ionic strength:   

Hfo_sOH + Sb(OH)6- + H+ = Hfo_sSbO(OH)4 + 2H2O   
   log_k   4.4779
   delta_h   0   kJ
      #                  Id:   8117410
      #        log K source:                                
      #      Delta H source:                                
      #T and ionic strength:   
Hfo_wOH + Sb(OH)6- + H+ = Hfo_wSbO(OH)4 + 2H2O
   log_k   4.4779
   delta_h   0   kJ
      #                  Id:   8127410
      #        log K source:                                
      #      Delta H source:                                
      #T and ionic strength:   

Hfo_sOH + Sb(OH)6- = Hfo_sOHSbO(OH)4- + H2O
   log_k   3.7350
   delta_h   0   kJ
      #                  Id:   8117411
      #        log K source:                                
      #      Delta H source:                                
      #T and ionic strength:   

Hfo_wOH + Sb(OH)6- = Hfo_wOHSbO(OH)4- + H2O
   log_k   3.7350
   delta_h   0   kJ
      #                  Id:   8127411
      #        log K source:                                
      #      Delta H source:                                
      #T and ionic strength:   

Hfo_sOH + Sb(OH)3 = Hfo_sSbO(OH)2 + H2O
   log_k   7
   delta_h   0   kJ
      #                  Id:   8117410
      #        log K source:                                
      #      Delta H source:                                
      #T and ionic strength:   

Hfo_wOH + Sb(OH)3 = Hfo_wSbO(OH)2 + H2O
   log_k   7
   delta_h   0   kJ
      #                  Id:   8127410
      #        log K source:                                
      #      Delta H source:                                
      #T and ionic strength:

I think it might be especially the reaction for Sb(3) sorption which is causing the failure.

I have never used the SAVE statement yet. I have read the explanation in the manual concerning the SAVE statement but it also tells me ‘ADVECTION, TRANSPORT, and RUN_CELLS calculations automatically save results after each calculation and the SAVE keyword has no effect for these calculations’. Do you mean I should use the SAVE statement after the SOLUTION data block and then use the USE statement before the TRANSPORT data block?

I actually just wanted to fixate both pH and pe since I don't have data about redox couples or others. When just specifying the pe and pH in the solution block, pe is still adjusted due to reaction with the surface – but I noticed in the output that during transport the reaction calculation that adjusts pe and pH is also not executed (defining the Fix_pe phase doesn’t seem to have any effect).. So is there a way to fix pe also during transport, or is this too forced? I furthermore thought (and experienced) that by indicating Fix_pe -2, a pe of 2 is imposed. Is this incorrect?

In the meanwhile I am still having the convergence problem. I noticed that for a pH of 3, I can reach  convergence if I define a pe of 1 or less. Could it have something to do with the charge balance (due to the high positive charge of the Hfo-surface)?

Thank you very much for the help!

Greetings,

Lisa

[dlparkhurst]

The problem is with very small concentrations of Sb. I don't think pe or redox is the problem. However, I am not an advocate of fixing the pe. I have no problem with fixing the partial pressure of oxygen; you will tend to stay within the stability field of water. But, if you fix pe and pH changes substantially, then you could range out of the stability field of water (check out a pe/pH diagram for the stability field of water).

Try using KNOBS; -tol 1e-20. Hopefully that works, but calculations with very small element concentrations is problematic with PHREEQC.

Concerning SAVE:

If you have the following (see below), then solution 0 (your infilling solution for TRANSPORT) has a pH of 3 and a pe of 4. The solution composition determined by the equilibration with EQUILIBRIUM_PHASES is not saved. If you add "SAVE solution 0" before the END statement, then solution 0 has a pH of 3 and a pe of 2.

SOLUTION 0
pH 3
Na 10
Cl 10
Sb 1E-03 # 10^-3 mmol = 1 µmol
EQUILIBRIUM_PHASES
Fix_H+ -3 NaOH 
   -force_equality true
Fix_pe -2 O2(g)   10
   -force_equality true
END

[llanduyt]

Thank you very much! By using KNOBS the convergence problem is solved. Saving solution 0 also solved the fluctuations of pe in the beginning of the simulation.

I have one last question about fixing the partial pressure of oxygen. Would this than be with the following statement?

EQUILIBRIUM_PHASES
O2(g)   -61.9951
   -force_equality true

Is it correct that I calculate the partial pressure I have to impose for a given pe (e.g. if I want to fix pe 3) with the following formula?
pO2 = -85.9951 + 4*(pH + pe)
based on the equation (assuming other oxidation reactions are negligible):
2H2O =  O2 + 4H+ + 4e-     
    log_k      -85.9951
        -analytic   38.0229    7.99407E-03   -2.7655e+004  -1.4506e+001  199838.45
Or should I also take into account the oxidation reactions of the other species present in the system?

Thank you very much for all the clarifications.

[dlparkhurst]

I still think you should not fix the pe. Fixing the pO2 will suffice to set the redox environment.

By including

EQUILIBRIUM_PHASES
O2(g)   -61.9951  10

you will add or remove O2 to maintain the specified log partial pressure of O2. Many other redox reactions may occur, but O2 will continue to be added or removed until the partial pressure is satisfied, or all 10 moles of O2 has reacted.