Trying to integrate incongruent dissolution into a reactive transport simulation

[wpowell9]

Hello:

I am trying to integrate the incongruent dissolution reaction of hannebachite (calcium sulfite hemihydrate) to gypsum when an acid solution is passed through the solid into a reactive transport equation. The simulation doe snot immediately error out but appears stuck on the first timestep during the equilibration (solution 1) reaction. Gypsum uses simple dissolution rate, KINETICS code for hannebachite is as follows:
    Hannebachite
-start
10 REM store initial amount of hannebachite
20 IF EXISTS(1) = 0 THEN PUT(M, 1)
30 REM calculate moles of reaction
40 SR_han = SR("Hannebachite")
50 IF (M<=0 and SR_han<1) THEN GOTO 130
60 REM PARM(1)=k(dissolution-neutral)
70 REM PARM(2)=A/V
80 REM PARM(3)=correction factor
90 moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_han) * TIME
100 REM Convert to Gypsum
110 IF ABS(SI("Gypsum")) > 1e-3 THEN GOTO 130
130 SAVE moles
-end

How to fix?
 

[dlparkhurst]

If you attach a complete file that demonstrates the problem (and database if necessary), I will take a look.

[wpowell9]

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.6.1-15000\database\phreeqc.dat
TITLE Title: Column Test of sFGD with AMD
      # Use Phreeqc.dat
      # Double-check time step. Convert in Excel to L/S ratio for outputs, redo simulation as 10-cell with outputs at 5 (2 L/S) and 10 (1 L/S)  Add Epsomite and add dynamic kinetics for same. Check for amorphous silica in cement models - may be useful for aluminum values.

PRINT
    -user_print            true
    -selected_output       true
    -status                false
    -warnings              1000
SOLUTION 0 AMD
    temp      25
    pH        3.18
    pe        7.61
    redox     pe
    units     mg/l
    density   1
    Ca        351.254
    Cl        11.7745
    K         3.5
    Mg        59.6634
    Na        4.4
    S(6)      1096.74
    Si        14.1943
    Fe        33.0487
    -water    1 # kg
EQUILIBRIUM_PHASES 3
    CO2(g)    -3.5 1000
    O2(g)     -0.7 1000
    Gypsum     0.0 Ca(SO3):0.5(H2O) 10.0
SAVE solution 0
END
PHASES
Brucite
    Mg(OH)2 + 2H+ = 2H2O + Mg+2
    log_k     16.298
    delta_h   -111.34 kJ
Hannebachite
    Ca(SO3):0.5(H2O) = Ca+2 + 0.5H2O + SO3-2
    log_k     231.533  #log_k and delta_h from Thermochemical Data of Pure Substances, 3rd ed. (Barin, 1995)
    delta_h   1311.701 kJ
Ettringite
    Ca6Al2(SO4)3(OH)12:26H2O + 12H+ = 2Al+3 + 6Ca+2 + 38H2O + 3SO4-2
    log_k     62.5362
    delta_h   -382.451 kJ
    -analytical_expression -1057 -0.11585 59580 385.85 1012.1 0
#Epsomite
#    MgSO4:7H2O  =  + 1.0000 Mg++ + 1.0000 SO4-- + 7.0000 H2O
#    log_k           -1.9623
#    -delta_H   0            # Not possible to calculate enthalpy of reaction   Epsomite
#   Enthalpy of formation:   0 kcal/mol
Magnetite
    Fe3O4 + 8H+ = Fe+2 + 2Fe+3 + 4H2O
    log_k     10.4724
    delta_h   -216.597 kJ
    -analytical_expression -305.1 -0.079919 18709 111.78 292.03 0
Sillimanite
    Al2SiO5 + 6H+ = 2Al+3 + H2O + H4SiO4
    log_k     16.308
    delta_h   -238.442 kJ
    -analytical_expression -71.61 -0.032196 12493 22.449 194.96 0
Portlandite
    Ca(OH)2 + 2H+ = Ca+2 + 2H2O
    log_k     22.5552
    delta_h   128.686 kJ
Magnesite
    MgCO3 +1.0000 H+  =  + 1.0000 HCO3- + 1.0000 Mg++
    log_k       2.2936
    -delta_H   -44.4968 kJ
    -analytical_expression -1.6665e+002 -4.9469e-002 6.4344e+003 6.5506e+001 1.0045e+002
SOLUTION_SPECIES
SO4-2 = SO3-2 + 0.5O2
    log_k     -46.6244
    delta_h   267.985 kJ
    -analytical_expression 13.771 0.0006512 -13330 4.7164 -208 0
SOLUTION 1-6 Presaturated with AMD in Column
    temp      25
    pH        3.18
    pe        7.61
    redox     pe
    units     mg/l
    density   1
    Ca        351.254
    Cl        11.7745
    Fe        33.0487
    K         3.5
    Mg        59.6634
    Na        5.4
    S(6)      1096.74
    Si        14.19
    -water    0.01 # kg
TRANSPORT
    -cells                 6
    -shifts                500
    -time_step             14400 # seconds
    -lengths               6*0.05
    -dispersivities        6*0.03
    -correct_disp          true
    -diffusion_coefficient 3e-09
    -thermal_diffusion     2   3e-09
    -print_cells           3 6
    -punch_cells           3 6
KINETICS 1-6 Kinetic reactions to all cells
Calcite
    -formula  CaCO3  1
    -m        0.0167 # increased from 0.02 to test pH stabilization
    -m0       0.0167
    -parms    1.55e-08 1360 0.1 2.17e-1 0.501 6.71e-9 3.31e-4 1.02e-5 # USGS
    -tol      1e-08
Gypsum
    -formula  CaSO4:2H2O  1
    -m        0.0194
    -m0       0.0194
    -parms    1.62e-05 1360 0.1 1.55
    -tol      1e-08
Magnesite
    -formula  MgCO3  1
    -m        0.0593
    -m0       0.0593
    -parms    4.571e-12 1360 0.1 6.3e+18 4.17e-7 6.91e+13 6.03e-6 4.78e+12 # USGS
    -tol      1e-08
Dolomite
    -formula  CaMg(CO3)2  1
    -m        0.0001
    -m0       0.0001
    -parms    2.951e-10 1360 0.1 6.78e-8 6.46e-4 3.10e-14 7.76e-6 2.58e-12 # USGS
    -tol      1e-08
Brucite
    -formula  Mg(OH)2  1
    -m        0.0343
    -m0       0.0343
    -parms    5.754e-11 1360 0.1 #
    -tol      1e-08
Halite
    -formula  NaCl  1
    -m        0.0427
    -m0       0.0427
    -parms    0.006167 1360 0.1
    -tol      1e-08
Portlandite
    -formula  Ca(OH)2  1
    -m        0.0010
    -m0       0.0010
    -parms    1e-05 1360 0.1
    -tol      1e-08
Ettringite
    -formula  Ettringite  1
    -m        0.0001
    -m0       0.0001
    -parms    5.0e-12 1360 0.1 1.78e-33 # Baur (2004) low end is 5e-12 orig 7.08e-13 (Perkins, 1999)
    -tol      1e-08
Hematite
    -formula  Fe2O3  1
    -m        0.0626
    -m0       0.0626
    -parms    2.512e-17 1360 0.1
    -tol      1e-08
Hannebachite
    -formula  Ca(SO3):0.5(H2O)  1
    -m        0.3871
    -m0       0.3871
    -parms    2.1262e-07 1360 0.1 1120
    -tol      1e-08
Magnetite
    -formula  Magnetite  1
    -m        0.0298
    -m0       0.0298
    -parms    4.571e-11 1360 1
    -tol      1e-08
Quartz
    -formula  SiO2  1
    -m        0.0277
    -m0       0.0277
    -parms    1.023e-16 1360 1
    -tol      1e-08
Sillimanite
    -formula  Sillimanite  1
    -m        0.0463
    -m0       0.0463
    -parms    3.162e-13 1360 1
    -tol      1e-08
 Gibbsite
    -formula  Gibbsite  1
    -m        0.000001
    -m0       0.000001
    -parms    3.162e-14 1360 0.1 2.00e-18 2.24e-10 2.82e-22 2.24e-19 2.82e-13 # orig 1.38e-6 Benezeth, Geochem Acta, 2008
    -tol      1e-08
-steps       5760000 in 500 steps # seconds
-step_divide 1000
-runge_kutta 3
-bad_step_max 50000
-cvode true
-cvode_steps 100
-cvode_order 2
INCREMENTAL_REACTIONS True
RATES
    Calcite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_cc = SR("Calcite")
30 if (M<=0 and SR_cc<1) then goto 350
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 REM PARM(4)=k(precipitation-neutral)
80 REM PARM(5)=k(dissolution-acid)
90 REM PARM(6)=k(precipitation-acid)
100 REM PARM(7)=k(dissolution-alkali)
110 REM PARM(8)=k(precipitation-alkali)
120 IF SR_cc <= 1 THEN
130   IF act("H+")  >= 1e-6 THEN
140     moles = PARM(5) * PARM(2)* PARM(3)*(M/M0)^0.67 * act("H+") * (1-SR_cc) * TIME # pH component )
150     goto 350
160   ELSE IF act("H+")  < 1e-6 AND [H+] >= 1e-8 THEN
170     moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_cc) * TIME # pH component
180     goto 350
190   ELSE IF act("H+")  < 1e-8 THEN
200     moles = PARM(7) * PARM(2)* PARM(3)*(M/M0)^0.67 * -1*act("H+") * (1-SR_cc) * TIME # pH component
210     goto 350
220   END IF
230 ELSE IF SR_cc > 1 THEN
240   IF act("H+")  >= 1e-6 THEN
250     moles = PARM(6) * PARM(2)* PARM(3)*(M/M0)^0.67 * act("H+") * (1-SR_cc) * TIME # pH component 
260     goto 350
270   ELSE IF act("H+")  < 1e-6 AND [H+] >= 1e-8 THEN
280     moles = PARM(4) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_cc) * TIME # pH component
290     goto 350
300   ELSE IF act("H+") < 1e-8 THEN
310     moles = PARM(8) * PARM(2)* PARM(3)*(M/M0)^0.67 * -1*act("H+") (1-SR_cc) * TIME # pH component
320     goto 350
330   END IF
340 END IF
350 save moles
-end
    Brucite
-start
10 REM use the simple dissolution rate from trainsition state theory
20 SR_br = SR("Brucite")
30 if (M<=0 and SR_br<1) then goto 70
40 REM PARM(1)=k
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_br) * TIME
80 SAVE moles
-end
    Magnesite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_ms = SR("Magnesite")
30 if (M<=0 and SR_ms<1) then goto 350
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 REM PARM(4)=k(precipitation-neutral)
80 REM PARM(5)=k(dissolution-acid)
90 REM PARM(6)=k(precipitation-acid)
100 REM PARM(7)=k(dissolution-alkali)
110 REM PARM(8)=k(precipitation-alkali)
120 IF SR_ms <= 1 THEN
130   IF act("H+")  >= 1e-6 THEN
140     moles = PARM(5) * PARM(2)* PARM(3)*(M/M0)^0.67 * act("H+") * (1-SR_ms)^4 * TIME # pH component Benezeth ((10^-14)/act("H+"))
150     goto 350
160   ELSE IF act("H+")  < 1e-6 AND [H+] >= 1e-8 THEN
170     moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_ms) * TIME # pH component Benezeth
180     goto 350
190   ELSE IF act("H+")  < 1e-8 THEN
200     moles = PARM(7) * PARM(2)* PARM(3)*(M/M0)^0.67 * -1*act("H+")4 * (1-SR_ms) * TIME # pH component Benezeth
210     goto 350
220   END IF
230 ELSE IF SR_ms > 1 THEN
240   IF act("H+")  >= 1e-6 THEN
250     moles = PARM(6) * PARM(2)* PARM(3)*(M/M0)^0.67 * act("H+") * (1-SR_ms)^4 * TIME # pH component Benezeth
260     goto 350
270   ELSE IF act("H+")  < 1e-6 AND [H+] >= 1e-8 THEN
280     moles = PARM(4) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_ms) * TIME # pH component Benezeth
290     goto 350
300   ELSE IF act("H+") < 1e-8 THEN
310     moles = PARM(8) * PARM(2)* PARM(3)*(M/M0)^0.67 * -1*act("H+") (1-SR_ms) * TIME # pH component Benezeth
320     goto 350
330   END IF
340 END IF
350 save moles
-end
    Dolomite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_do = SR("Dolomite")
30 if (M<=0 and SR_do<1) then goto 350
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 REM PARM(4)=k(precipitation-neutral)
80 REM PARM(5)=k(dissolution-acid)
90 REM PARM(6)=k(precipitation-acid)
100 REM PARM(7)=k(dissolution-alkali)
110 REM PARM(8)=k(precipitation-alkali)
120 IF SR_do <= 1 THEN
130   IF act("H+")  >= 1e-6 THEN
140     moles = PARM(5) * PARM(2)* PARM(3)*(M/M0)^0.67 * act("H+")^0.5 * (1-SR_do) * TIME # pH component 
150     goto 350
160   ELSE IF act("H+")  < 1e-6 AND [H+] >= 1e-8 THEN
170     moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_do) * TIME # pH component
180     goto 350
190   ELSE IF act("H+")  < 1e-8 THEN
200     moles = PARM(7) * PARM(2)* PARM(3)*(M/M0)^0.67 * -1*act("H+")^0.5 * (1-SR_do) * TIME # pH component 
210     goto 350
220   END IF
230 ELSE IF SR_do > 1 THEN
240   IF act("H+")  >= 1e-6 THEN
250     moles = PARM(6) * PARM(2)* PARM(3)*(M/M0)^0.67 * act("H+")^0.5 * (1-SR_do) * TIME # pH component
260     goto 350
270   ELSE IF act("H+")  < 1e-6 AND [H+] >= 1e-8 THEN
280     moles = PARM(4) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_do) * TIME # pH component
290     goto 350
300   ELSE IF act("H+") < 1e-8 THEN
310     moles = PARM(8) * PARM(2)* PARM(3)*(M/M0)^0.67 * -1*act("H+")^0.5 (1-SR_do) * TIME # pH component
320     goto 350
330   END IF
340 END IF
350 save moles
-end
    Gypsum
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_gy = SR("Gypsum")
30 if (M<=0 and SR_gy<1) then goto 120
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
90 moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_gy) * TIME
120 SAVE moles
-end
    Hannebachite
-start
10 REM store initial amount of hannebachite
20 IF EXISTS(1) = 0 THEN PUT(M, 1)
30 REM calculate moles of reaction
40 SR_han = SR("Hannebachite")
50 IF (M<=0 and SR_han<1) THEN GOTO 130
60 REM PARM(1)=k(dissolution-neutral)
70 REM PARM(2)=A/V
80 REM PARM(3)=correction factor
90 moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_han) * TIME
100 REM Convert to Gypsum
110 IF ABS(SI("Gypsum")) > 1e-3 THEN GOTO 130
130 SAVE moles
-end
    Hematite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_he = SR("Hematite")
30 if (M<=0 and SR_he<1) then goto 70
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_he) * TIME
80 SAVE moles
-end
    Magnetite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_mg = SR("Magnetite")
30 if (M<=0 and SR_mg<1) then goto 70
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_mg) * TIME
80 SAVE moles
-end
    Quartz
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_qz = SR("Quartz")
30 if (M<=0 and SR_qz<1) then goto 70
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_qz) * TIME
80 SAVE moles
-end
    Sillimanite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_sm = SR("Sillimanite")
30 if (M<=0 and SR_sm<1) then goto 70
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_sm) * TIME
80 SAVE moles
-end
    Ettringite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_et = SR("Ettringite")
30 if (M<=0 and SR_et<1) then goto 150
40 REM PARM(1)=k(dissolution)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 REM PARM(4)=k(precipitation)
80 IF SR_et <= 1 THEN
90    moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_et) * TIME
100   goto 150
110 ELSE IF SR_et >= 1 THEN
120   moles = PARM(4) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_et) * TIME
130   goto 150
140 END IF
150 save moles
-end
    Portlandite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_pt = SR("Portlandite")
30 if (M<=0 and SR_pt<1) then goto 70
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_pt) * TIME
80 SAVE moles
-end
    Halite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_ha = SR("Halite")
30 if (M<=0 and SR_ha<1) then goto 70
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_ha) * TIME
80 SAVE moles
-end
    Gibbsite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_gb = SR("Gibbsite")
30 if (M<=0 and SR_gb<1) then goto 350
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 REM PARM(4)=k(precipitation-neutral)
80 REM PARM(5)=k(dissolution-acid)
90 REM PARM(6)=k(precipitation-acid)
100 REM PARM(7)=k(dissolution-alkali)
110 REM PARM(8)=k(precipitation-alkali)
120 IF SR_gb <= 1 THEN
130   IF act("H+")  >= 1e-6 THEN
140     moles = PARM(5) * PARM(2)* PARM(3)*(M/M0)^0.67 * act("H+")^0.992 * (1-SR_gb) * TIME # pH component Benezeth ((10^-14)/act("H+"))
150     goto 350
160   ELSE IF act("H+")  < 1e-6 AND [H+] >= 1e-8 THEN
170     moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_gb) * TIME # pH component Benezeth
180     goto 350
190   ELSE IF act("H+")  < 1e-8 THEN
200     moles = PARM(7) * PARM(2)* PARM(3)*(M/M0)^0.67 * -1*act("H+")^-0.784 * (1-SR_gb) * TIME # pH component Benezeth
210     goto 350
220   END IF
230 ELSE IF SR_gb > 1 THEN
240   IF act("H+")  >= 1e-6 THEN
250     moles = PARM(6) * PARM(2)* PARM(3)*(M/M0)^0.67 * act("H+")^0.992 * (1-SR_gb) * TIME # pH component Benezeth
260     goto 350
270   ELSE IF act("H+")  < 1e-6 AND [H+] >= 1e-8 THEN
280     moles = PARM(4) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_gb) * TIME # pH component Benezeth
290     goto 350
300   ELSE IF act("H+") < 1e-8 THEN
310     moles = PARM(8) * PARM(2)* PARM(3)*(M/M0)^0.67 * -1*act("H+")^-0.784 (1-SR_gb) * TIME # pH component Benezeth
320     goto 350
330   END IF
340 END IF
350 save moles
-end
SELECTED_OUTPUT 1
    -file                 C:\OSU CEG & GEO\Research\6-cell models\Gavin 1314_6cells_base_statbruc_noalk_lowport.sel
    -totals               Ca  Al  Cl  Mg  Na  S(6) Fe  Si
    -equilibrium_phases   Calcite  Halite  Portlandite  Gypsum
                          Dolomite  Magnesite  Brucite  Hannebachite
                          Hematite  Quartz  Manganite  Gibbsite Ettringite #Epsomite
    -saturation_indices   Brucite  Calcite  Gibbsite  Portlandite
                          Gypsum  Dolomite  Magnesite Hannebachite Ettringite #Epsomite
    -kinetic_reactants    Gypsum  Dolomite  Calcite  Brucite
                          Halite  Magnesite  Portlandite  Quartz
                          Gibbsite Hannebachite Ettringite #Epsomite
USER_GRAPH 1 pH and Ca
    -headings               LS Ca SO4 Cl Na Mg Fe Al pH
    -axis_titles            "Liquid to Solid Ratio" "Concentration (mg/L)" "pH"
    -chart_title            "Gavin Fixated FGD Material"
    -axis_scale x_axis      auto auto auto auto log
    -axis_scale y_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X (STEP_NO +0.5) * 14400 / 86400
20 GRAPH_Y TOT("Ca") * 40e3
30 GRAPH_Y TOT("S(6)") * 96e3
40 GRAPH_Y TOT("Cl") * 35.5e3
50 GRAPH_Y TOT("Na") * 23e3
60 GRAPH_Y TOT("Mg") * 24.3e3
70 GRAPH_Y TOT("Fe") * 55.9e3
80 GRAPH_Y TOT("Al") * 27e3
90 GRAPH_SY -LA("H+")
  -end
    -active                 true
USER_GRAPH 2 pH and Minerals
    -headings               LS Halite Magnesite Calcite Brucite Portlandite Gypsum Dolomite Ettringite Gibbsite pH #Epsomite
    -axis_titles            "Liquid to Solid Ratio" "Concentration (mole/L)" "pH"
    -chart_title            "Carmeuse LWM"
    -axis_scale x_axis      auto auto auto auto log
    -axis_scale y_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 GRAPH_X (STEP_NO +0.5) * 14400 / 86400
 20 GRAPH_Y KIN("Halite")
 30 GRAPH_Y KIN("Magnesite")
 40 GRAPH_Y KIN("Calcite")
 50 GRAPH_Y KIN("Brucite")
 60 GRAPH_Y KIN("Portlandite")
 70 GRAPH_Y KIN("Gypsum")
 80 GRAPH_Y KIN("Dolomite")
 90 GRAPH_Y KIN("Ettringite")
100 GRAPH_Y KIN("Gibbsite")
#110 GRAPH_Y KIN("Epsomite")
120 GRAPH_SY -LA("H+")
  -end
    -active                 true
END
END

[dlparkhurst]

Your rates definitions are very logical, but unfortunately they are not correct Basic. There is no statement END IF, and an IF/THEN must be completed on a single line.

Attached is a rewrite of the Gibbsite rate, that should work, barring any of my mistakes.

I suggest that you simplify your problem initially to react SOLUTION 0 with KINETICS 1 and plot the results versus TOTAL_TIME. It should give you a feel for the reactions without the complications of TRANSPORT.

Note that 6 shifts will advect solution 0 through the entire column, so when you do run TRANSPORT, solution 0 is only in contact with the kinetic reactants for 6 time steps, which will limit the extent of reaction in a single parcel of water.

[wpowell9]

Have updated the underlying transport script (sans hannebachite incongruent dissolution) with proper IF statements(all on one line). Now have issues when running with dynamic precipitation/dissolution of gibbsite and/or dolomite - fails to converge. Dolomite works for a simplified case where precipitation and dissolution are dynamic but do not involve rate changes due to pH (semistat_dolm) whereas trying to run full pH rate change case (testfull) fails. Similar problem when dynamic brucite attempted. Codes attached. Wondering direction to go to fix, then can return to sorting out incongruent dissolution issue.

[dlparkhurst]

I just looked at the dynbruc example.

For starters, the gibbsite and dolomite RATES definitions were the old incorrect Basic definitions.

But, I think your main problem is that the parameters change so abruptly. For brucite, the rate change between undersaturation and supersaturation is about 40 orders of magnitude. Doesn't seem reasonable, but even if it is, the numerical method is going to have difficulties. CVODE uses the derivatives of the rates, and brucite will be non-differentiable at SR= 1. Given that brucite approaches equilibrium, there are huge rate variations for stepping from a little bit under to a little bit over saturation (or the opposite). I suspect that the result is an oscillation getting close to equilibrium from one side, and then way over shooting when you step to the other side.

I know it is elegant to use kinetics for all minerals, but the calculation will be more robust and faster if you can substitute equilibrium conditions for brucite, dolomite, and (or) gypsum.

[wpowell9]

Dr. Parkhurst:

Upthread, you noted:

"Note that 6 shifts will advect solution 0 through the entire column, so when you do run TRANSPORT, solution 0 is only in contact with the kinetic reactants for 6 time steps, which will limit the extent of reaction in a single parcel of water. "

The code is simulating a large volume of AMD being flowed over time through a column of solids. After 6 timesteps, the AMD which entered the column in cell 1 exits cell 6. New AMD is inflowing as old AMD exits. If the code as written does not reflect that setup, please state so, otherwise please confirm that code is fit to above-stated purpose.

(i.e. I think what you are trying to say is that any given aliquot of AMD will only react for 6 timesteps. In this case, this is by design, as new AMD is entering as old AMD is leaving)

[dlparkhurst]

I think you understand the way transport works correctly.

[wpowell9]

New code attached with brucite, gibbsite, and gypsum/hannebachite running using equilibrium phases instead of kinetics. New question: Would very much like to set a starting gypsum population, and would like to have per-timestep populations of each along with the kinetics output. How to accomplish?

After getting this to work, will likely set dolomite to equilibrium phases as well as empirical results do not show persistent dolomite formation, but simulation does. Suspect this is due to the quasi-kinetic model currently in place. Need to be able to track formation per-timestep as with gypsum, brucite,e tc.

[dlparkhurst]

A few comments.

You are running 500 steps of KINETICS 1 reacting with SOLUTION 1 before the the TRANSPORT calculation, so KINETICS 1 is different that KINETICS 2-6 at the beginning of TRANSPORT.

Cell 2, only, has EQUILIBRIUM_PHASES during transport.

SOLUTION 0 is equilibrated with brucite, so it has a pH of 9 for your infilling solution. This seems at odds with your stated purpose of "dissolution reaction of hannebachite (calcium sulfite hemihydrate) to gypsum when an acid solution is passed through".

You should probably define SO3-2 as a new redox state, otherwise some of the output looks strange (negative O(0)).

Code: [Select]

SOLUTION_MASTER_SPECIES
S(4) SO3-2 0 12 12

You start with 0.225 kgw in cells 1-6, but solution 0 has 1 kgw, so after 6 shifts, each cell has 1 kgw. I prefer to scale the reactants in the column to 1 kgw, but you could simply make solution 0 the same volume as the cells.

Dolomite can generally be assumed not to precipitate.

SELECTED_OUTPUT and USER_PUNCH can be used to write simulation results to file.

But my basic problem is that this simulation is way too complicated to debug as it is. You said you wanted to dissolve Hannebachite with acid. Here is a simpler simulation that is the way I interpret what you are doing. I modified the RATES definition for Hannebachite to allow dissolution. With this rate, one time step is sufficient to dissolve most of the Hannebachite. Sorry if this is not close to what you had in mind.

Code: [Select]

RATES
    Calcite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_cc = SR("Calcite")
30 if (M<=0 and SR_cc<1) then goto 180
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
70 REM PARM(4)=k(precipitation-neutral)
80 REM PARM(5)=k(dissolution-acid)
90 REM PARM(6)=k(precipitation-acid)
100 REM PARM(7)=k(dissolution-alkali)
110 REM PARM(8)=k(precipitation-alkali)
120 IF SR_cc <= 1 AND act("H+") >= 1e-6 THEN moles = PARM(5) * PARM(2)* PARM(3)*(M/M0)^0.67 * act("H+") * (1-SR_cc) * TIME # pH component
130 IF SR_cc <= 1 AND act("H+") < 1e-6 AND act("H+") >= 1e-8 THEN moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_cc) * TIME # pH component
140 IF SR_cc <= 1 AND act("H+") < 1e-8 THEN moles = PARM(7) * PARM(2)* PARM(3)*(M/M0)^0.67 * -1*act("H+") * (1-SR_cc) * TIME # pH component
150 IF SR_cc > 1 AND act("H+") >= 1e-6 THEN moles = PARM(6) * PARM(2)* PARM(3)*(M/M0)^0.67 * act("H+") * (1-SR_cc) * TIME # pH component 
160 IF SR_cc > 1 AND act("H+") < 1e-6 AND act("H+") >= 1e-8 THEN moles = PARM(4) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_cc) * TIME # pH component
170 IF SR_cc > 1 AND act("H+") < 1e-8 THEN moles = PARM(8) * PARM(2)* PARM(3)*(M/M0)^0.67 * -1*act("H+") * (1-SR_cc) * TIME # pH component
180 save moles
-end
    Hannebachite
-start
10 REM use the simple dissolution rate from transition state theory
20 SR_han = SR("Hannebachite")
#30 if (M<=0 and SR_han<1) then goto 140
40 REM PARM(1)=k(dissolution-neutral)
50 REM PARM(2)=A/V
60 REM PARM(3)=correction factor
110 moles = PARM(1) * PARM(2)* PARM(3)*(M/M0)^0.67 * (1-SR_han) * TIME
140 SAVE moles
-end
PHASES
Brucite
    Mg(OH)2 + 2H+ = 2H2O + Mg+2
    log_k     16.298
    delta_h   -111.34 kJ
Hannebachite
    Ca(SO3):0.5(H2O) = Ca+2 + 0.5H2O + SO3-2
    log_k     6.51
Ettringite
    Ca6Al2(SO4)3(OH)12:26H2O + 12H+ = 2Al+3 + 6Ca+2 + 38H2O + 3SO4-2
    log_k     62.5362
    delta_h   -382.451 kJ
    -analytical_expression -1057 -0.11585 59580 385.85 1012.1 0
Magnetite
    Fe3O4 + 8H+ = Fe+2 + 2Fe+3 + 4H2O
    log_k     10.4724
    delta_h   -216.597 kJ
    -analytical_expression -305.1 -0.079919 18709 111.78 292.03 0
Sillimanite
    Al2SiO5 + 6H+ = 2Al+3 + H2O + H4SiO4
    log_k     16.308
    delta_h   -238.442 kJ
    -analytical_expression -71.61 -0.032196 12493 22.449 194.96 0
Portlandite
    Ca(OH)2 + 2H+ = Ca+2 + 2H2O
    log_k     22.5552
    delta_h   128.686 kJ
Magnesite
    MgCO3 +1.0000 H+  =  + 1.0000 HCO3- + 1.0000 Mg++
    log_k       2.2936
    -delta_H -44.4968 kJ
    -analytical_expression -1.6665e+002 -4.9469e-002 6.4344e+003 6.5506e+001 1.0045e+002
SOLUTION_MASTER_SPECIES
S(4) SO3-2 0 12 12
SOLUTION_SPECIES
SO4-2 = SO3-2 + 0.5O2
    log_k     -46.6244
    delta_h   267.985 kJ
    -analytical_expression 13.771 0.0006512 -13330 4.7164 -208 0
END
SOLUTION 0 AMD
    temp      25
    pH        3.18
    pe        7.61
    redox     pe
    units     mg/l
    density   1
    Ca        351.254
    Cl        11.7745
    K         3.5
    Mg        59.6634
    Na        4.4
    S(6)      1096.74
    Si        14.1943
    Fe        33.0487
    -water    0.2225
EQUILIBRIUM_PHASES 0
    CO2(g)    -3.5 1000
    O2(g)     -0.7 1000
SAVE solution 0
END
SOLUTION 1-6 Presaturated with AMD in Column
    temp      25
    pH        3.18
    pe        7.61
    redox     pe
    units     mg/l
    density   1
    Ca        351.254
    Cl        11.7745
    Fe        33.0487
    K         3.5
    Mg        59.6634
    Na        5.4
    S(6)      1096.74
    Si        14.19
    -water    0.2225
END
EQUILIBRIUM_PHASES 1-6
Gypsum 0 0
END
KINETICS 1-6
-cvode
Calcite
    -formula  CaCO3  1
    -m        0.0167 # increased from 0.02 to test pH stabilization
    -m0       0.0167
    -parms    1.55e-08 1360 0.1 2.17e-1 0.501 6.71e-9 3.31e-4 1.02e-5 # USGS
    -tol      1e-08
Hannebachite
    -formula  Ca(SO3):0.5(H2O)  1
    -m        0.3871
    -m0       0.3871
    -parms    2.1262e-7 1360 0.1
    -tol      1e-08
END
TRANSPORT
    -cells                 6
    -shifts                1 # 500
    -bc                    flux flux
    -time_step             14400 # seconds
    -lengths               6*0.05
    -dispersivities        0.005 # 6*0.03 # Runs faster for now
    -correct_disp          true
    -diffusion_coefficient 3e-09
    -thermal_diffusion     2   3e-09
    -punch_frequency       1
    -punch_cells           1-6
    -print_cells           1-6
USER_GRAPH 1 pH and Ca
    -headings               DIST Ca S(6) pH
    -axis_titles            "Liquid to Solid Ratio" "Concentration (mg/L)" "pH"
    -chart_title            "Gavin Fixated FGD Material"
    #-axis_scale x_axis      auto auto auto auto log
    -axis_scale y_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X DIST
20 GRAPH_Y TOT("Ca") * 40e3
30 GRAPH_Y TOT("S(6)") * 96e3
90 GRAPH_SY -LA("H+")
  -end
USER_GRAPH 2 pH and Minerals
    -headings               DIST KIN_Calc KIN_Hann EQUI_Gyp SI_Calc SI_Gyp
    -axis_titles            "Distance" "MOLES" "SATURATION INDEX"
    -chart_title            "Carmeuse LWM"
    #-axis_scale x_axis      auto auto auto auto log
    -axis_scale sy_axis      -10 2 auto auto
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 GRAPH_X DIST
 40 GRAPH_Y KIN("Calcite")
 50 GRAPH_Y KIN("Hannebachite")
 60 GRAPH_Y EQUI("Gypsum")
 70 GRAPH_SY SI("Calcite")
 90 GRAPH_SY SI("Gypsum")
  -end
    -active                 true


[wpowell9]

Hannebachite dissolving rapidly I can get to happen. Th issue is that the experimental results look like hannebachite is converting to gypsum (based on XRD data), which is part of what I am trying to simulate with this.

[dlparkhurst]

sit.dat has a CaSO3 mineral with log K -6.5; do you have a sign error on your Hannebachite log K definition?

To make gypsum, you need to increase calcium and (or) sulfate concentrations in the cells of your column. Your input solution could be the source, but I assume it is well defined. Calcite dissolution was included in the latest model, but did not produce enough extra calcium to precipitate gypsum. Lacking other sources, oxidation of sulfite seems the most likely reaction. However, it takes more than just oxygen saturation in the infilling solution; perhaps diffusion of oxygen into the column or contact with the atmosphere? I guess it is possible your Hannebachite is already oxidized to gypsum to some extent when you make the column.

[wpowell9]

Stupid question for using equilibrium phase calculations in conjunction with kinetic code:

Given:
EQUILIBRIUM_PHASES 0
    CO2(g)    -3.5 1000
    O2(g)     -0.7 1000
EQUILIBRIUM_PHASES 1-6
    Brucite    0.0 0.0343
    Gibbsite   0.0 0.00001

would the concentrations of brucite and gibbsite in cells 1-6 reset at each timestep to initial vales shown?

[dlparkhurst]

Assuming you are using TRANSPORT, no for cells 1-6; the initial conditions apply at time 0, but after each time step, the moles of equilibrium phases remaining in cell 1-6 are saved to start the next time step.

As for solution 0, at each time step, the original solution 0 is reacted with equilibrium_phases 0 to produce the water composition that is advected into the column.

[wpowell9]

That's what I thought, but the question came up.