convert meq/kgw to meq/L in alkalinity and convergence error

[Chatura]

Hi all,

I am in the process of calculating the saturation index and to validate my model I have to add alkalinity to the model. But the value is in meq/kgw and has to be input in meq/L. I converted it and it gives a convergence error. Is there any coefficient to convert meq/kgw to meq/L.

And also I want to change the water/Rock ratio since I am modelling with rock. How to change the W/R ratio in the model.

Really appreciate your support.

[dlparkhurst]

For most solutions, meq/kgw and meq/L are numerically nearly equal.

SOLUTION allows you to pick the units for the analytical data. It can be ppm (mg/kgs), mg/L, mmol/L, or mmol/kgw (among others). However, you must enter all of the data with the same denominator /kgs, /L, or /kgw. So you can enter elements as mg/L and mol/L in the same SOLUTION definition, but you cannot enter mg/L and mg/kgw in the same SOLUTION definition.

If not all the data have the same denominator, you are not far off calling all of them /L or /kgw.

[dlparkhurst]

A few problems with your SOLUTION. First, at 250C, pH 7.4, and pe 4.0, your solution is outside the stability field of water. The O2(aq) would be some large number. The results of the script below indicate that a pe of 2.8 gives one atm partial pressure of O2(g), conventionally, the edge of the stability field of water.

Second, Alkalinity in SOLUTION is the total alkalinity, including contributions from all alkalinity species. In your solution, Fe and Si contribute more than 0.6 meq/L, without any carbon in system. PHREEQC will try to find the amount of carbon in your system, but it fails if it needs a negative contribution to alkalinity. You can define C(4), which is total dissolved inorganic carbon. I've not done a lot of calculations at 250 C, but if most of the carbon is HCO3-, then you can just use the alkalinity value for C(4); still you are on your on at this temperature.

Third, you need to be careful with density. If you enter the data as per liter, then the density is used to convert to per kilogram solution, and then the sum of the solutes is subtracted to calculate the mass of water in a kilogram of solution. The final step is to use the mass of water to arrive at moles per kilogram water, which is the unit used for all calculations. If you add "calc" after your density estimate the program will iterate to find the density in SOLUTION that gives the same calculated density shown in the output. In this case, it is 0.8 kg/L.

Perhaps you want to run your SOLUTION at 25 C, which may be the temperature at which the analysis took place, and then use REACTION_TEMPERATURE to heat it back up to 250 C, possibly with some EQUILIBRIUM_PHASE constrants on gases?

Code: [Select]

PRINT
-alkalinity true
SOLUTION 1
    temp      250
    pH        7.41
    pe        2  O2(g) 0
    redox     pe
    units     mmol/l
    density   1.1 calc
    Al        0.02
    #Alkalinity 0.6 meq/L
    Ca        3.69
    Fe        0.21
    K         6.82
    Mg        0.64
    Na        934.78
    Si        3.05
    #-water    20 # kg

I generally scale the rock (minerals) to one L of water, which at low temperature is pretty close to 1 kgw. In your case, 1 L of solution differs quite a bit from 1 kgw. I would put in the amount of minerals and other reactants that would be present in 1 L or 1 kgw, your choice and adjust the solution accordingly. You can scale the solution to a different volume with -water, with the limitation that PHREEQC works best if the mass of water in the solution is within a couple orders of magnitude of 1.0.

[dlparkhurst]

The definition of KINETICS was corrupted in your file, but I have rearranged it to make it work.

I like to define RATES, SOLUTION_SPECIES, and other database-like data blocks first in the script to ensure that the definitions apply to all of the calculations.

I have limited REACTION_TEMPERATURE and KINETICS to one step just to simplify. Note that the 5 steps defined in REACTION_TEMPERATURE will apply to the first 5 steps defined in KINETICS, so the temperature will increase from 25 to 250 C in the first 5 steps, and remain at 250 C for the remaining 95 steps of KINETICS. I added INCREMENTAL_REACTIONS to speedup a multistep KINETICS integration. By default KINETICS will integrate from time zero to each incremental time, so the first increment will be integrated 100 times in your case. INCREMENTAL_REACTIONS causes the reaction to be integrated from 0 to T1, T1 to T2, ... until the total time has been integrated.

Code: [Select]

RATES
Quartz
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_qtz = SI("Quartz")
    30 if (M <= 0 and si_qtz < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_qtz > 0) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 10^(-13.40)*EXP(-90.90e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_qtz))
    190 moles = r * TIME
    200 SAVE moles
     -end
Anorthite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_anort = SI("Anorthite")
   30 if (M <= 0 and si_anort < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_anort > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-3.5)*EXP(-16.60e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.411)                 
   70 k_neut = 10^(-9.82)*EXP(-31.50e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_anort))
   190 moles = r * TIME
   200 SAVE moles
   -end
Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_alb = SI("Albite")
    30 if (M <= 0 and si_alb < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_alb > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.16)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("H+")^(0.457))                 
    70 k_neut = 10^(-12.56)*EXP(-69.80e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 10^(-15.60)*EXP(-71.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("OH-")^(-0.572))
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_alb))
    190 moles = r * TIME
    200 SAVE moles
    -end
K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_kfeld = SI("K-Feldspar")
    30 if (M <= 0 and si_kfeld < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_kfeld > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.06)*EXP(-51.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)                 
    70 k_neut = 10^(-12.41)*EXP(-38.00e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 10^(-21.20)*EXP(-94.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.823)
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_kfeld))
    190 moles = r * TIME
    200 SAVE moles
    -end
Diopside
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_diopsi = SI("Diopside")
    30 if (M <= 0 and si_diopsi < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_diopsi > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-6.36)*EXP(-96.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.710)                 
    70 k_neut = 10^(-11.11)*EXP(-40.60e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_diopsi))
    190 moles = r * TIME
    200 SAVE moles
    -end
Enstatite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_enst = SI("Enstatite")
   30 if (M <= 0 and si_enst < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_enst > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-9.02)*EXP(-80.00e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.600)                 
   70 k_neut = 10^(-12.72)*EXP(-80.00e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_enst))
   190 moles = r * TIME
   200 SAVE moles
   -end
Magnetite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_magnet = SI("Magnetite")
   30 if (M <= 0 and si_magnet < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_magnet > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-8.59)*EXP(-18.60e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.279)                 
   70 k_neut = 10^(-10.78)*EXP(-18.60e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_magnet))
   190 moles = r * TIME
   200 SAVE moles
   -end
Epidote
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_epidt = SI("Epidote")
    30 if (M <= 0 and si_epidt < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_epidt > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.60)*EXP(-71.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.338)                 
    70 k_neut = 10^(-11.99)*EXP(-70.70e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 10^(-17.33)*EXP(-79.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.556)
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_epidt))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_calc = SI("Calcite")
    30 if (M <= 0 and si_calc < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_calc > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-0.30)*EXP(-14.40e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)                 
    70 k_neut = 10^(-5.81)*EXP(-23.50e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_carb = 10^(-3.48)*EXP(-35.40e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(1.000)
    90 k_rateconst = k_acid + k_neut + k_carb
    100 r = k_rateconst * SA * (1-(10^si_calc))
    190 moles = r * TIME
    200 SAVE moles
    -end
Illite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_ill = SI("Illite")
    30 if (M <= 0 and si_ill < 0) then goto 200           
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_ill > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.12)*EXP(-58.00e+03/8.314*(1.0/TK-1.0/373.15))*ACT("H+")^(0.550)                 
    70 k_neut = 10^(-12.26)*EXP(-54.00e+03/8.314*(1.0/TK-1.0/373.15))
    80 k_base = 10^(-10.60)*EXP(-77.00e+03/8.314*(1.0/TK-1.0/373.15))*ACT("OH-")^(-0.350)
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_ill))
    190 moles = r * TIME
    200 SAVE moles
    -end
Chamosite-7A
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_cham = SI("Chamosite-7A")
    30 if (M <= 0 and si_cham < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_cham > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-11.11)*EXP(-88.00e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)                 
    70 k_neut = 10^(-12.52)*EXP(-88.00e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_cham))
    190 moles = r * TIME
    200 SAVE moles
    -end
Dolomite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_dolo = SI("Dolomite")
    30 if (M <= 0 and si_dolo < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_dolo > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-3.76)*EXP(-56.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)                 
    70 k_neut = 10^(-8.60)*EXP(-95.30e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_carb = 10^(-5.37)*EXP(-45.70e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(0.500)
    90 k_rateconst = k_acid + k_neut + k_carb
    100 r = k_rateconst * SA * (1-(10^si_dolo))
    190 moles = r * TIME
    200 SAVE moles
    -end
Anhydrite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_anhyd = SI("Anhydrite")
   30 if (M <= 0 and si_anhyd < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_anhyd > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = 10^(-2.79)*EXP(-0.00e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_anhyd))
   190 moles = r * TIME
   200 SAVE moles
   -end
END
SOLUTION 1
    temp      25
    pH        7.41
    pe        2
    redox     pe
    units     mmol/l
    density   1.1
    Al        0.02
    Alkalinity 0.6 meq/l      as HCO3-
    Ca        3.69
    Fe        0.21
    K         6.82
    Mg        0.64
    Na        934.78
    Si        3.05
    -water    20 # kg
INCREMENTAL_REACTIONS true
REACTION_TEMPERATURE 1
    25 #250 in 5 steps

KINETICS 1
-steps       3024000 #in 100 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
Quartz
    -formula  SiO2  1
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
Anorthite
    -formula  CaAl2(SiO4)2  1
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
Albite
    -formula  NaAlSi3O8  1
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
K-Feldspar
    -formula  KAlSi3O8  1
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
Diopside
    -formula  CaMgSi2O6  1
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
Enstatite
    -formula  MgSiO3  1
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
Magnetite
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
Epidote
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
Calcite
    -formula  CaCO3  1
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
Illite
    -formula  K0.6Mg0.25Al2.3Si3.5O10(OH)2  1
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
Chamosite-7A
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
Dolomite
    -formula  CaMg(CO3)2  1
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
Anhydrite
    -formula  CaSO4  1
    -m        0
    -m0       0
    -parms    0
    #-tol      1e-12
END

The "negative moles for Al" is a warning. As the rate integration proceeds, sometimes a bad step is taken and an element concentration becomes negative. The method then backs up to the last good step and repeats with a smaller time step to avoid the negative concentration. So the message is just a "warning", and, and if the calculation finishes, the results of the calculation should be correct. If more serious problems occur, the calculation will fail with an "ERROR" message, and the calculation will not be correct or completed.

[dlparkhurst]

I am skeptical that these rate equations are applicable at 250 C, but if you use -cvode in KINETICS, the integrations can be completed.

[dlparkhurst]

Are you asking why the minerals are not at saturation? You are only running the simulation for about 0.1 years. You probably need to run for many thousands of years to approach equilibrium with some of these minerals.

Check out the minerals that are supersaturated at the end of longer runs to see if there are other reasonable products, especially in the iron system.

To speed the calculation, you may want to assume equilibrium with some phases, particularly calcite. You also may want remove some of the minerals that have zero mass and never form (have very small number of moles) from the KINETICS definition (or put them in EQUILIBRIUM_PHASES). I suspect the simulation would run much faster.