Mixing/Equillibrating gas with solution

[Mgill9460]

Hi there, I am new to PHREEQC and have a rather simple question.

I have a set of different basic solutions. In these solutions I have experimentally measured the rate at which CO2 is dissolved into them. I am hoping to replicate the experimental results by equillibraitng the solutions with 1L of gas similar to the atmospheric composition of gas.

When I try to run the script it says there are input errors but it doesn't specify which error I have. I would really appreciate if anyone could help me out or point in the direction of some resources. TYIA.

Input:

Title: Gas-solution equillibrium
Site: 1
Date: September 2023

SOLUTION 1
   units mmol/kg
   Temp 25.0
   pH 11.8
   density 1.0
   Li 0.000541931
   Mg 0.000925735
   Si 0.000178031
   Ca 0.322832976

GAS_PHASE 1
    -fixed_pressure
    -pressure 1
    -volume 1
    -temperature 25
    CO2(g)    0.00042
    N2(g)     0.78
    O2(g)     0.2

USE SOLUTION 1
USE GAS_PHASE 1

END

[dlparkhurst]

The error is that you need either kgs or kgw in the units definition for kilogram solution or kilogram water. Should be an error message, but apparently not.

Code: [Select]

units mmol/kgw

You should use phreeqc.dat, Amm.dat, or pitzer.dat for the database. These databases have Ntg(g) instead of N2(g), which is a non-reactive nitrogen N2. N2 will react with oxygen to make nitrate, which is probably not what you want. Ntg(g) will be inert and simply dissolve in water.

To dissolve CO2 kinetically, I think you will need to make some further definitions to distinquish CO2 in the gas phase from dissolved CO2. Check out the forum for some discussions. I think this thread may help https://phreeqcusers.org/index.php/topic,2260.msg8452.html#msg8452.

[Mgill9460]

Hi Dr. Parkhurst, thank you for taking the time to reply.

If I could explain a little furhter basically what my experiments do is measure CO2 gas concentrations in an isolated chamber (volume ~7L) over natural bodies of water which sequester CO2.

I am wondering what is a simple way to simulate this system using a finite volume of gas.

I specified the units and changed the nitrogen from N2(g) to Ntg(g) but now I get the error:
ERROR: Gas not found in PHASES data base, Ntg(g).
ERROR: Gas not found in PHASES data base, O2(g).
WARNING: Could not find element in database, Si.
   Concentration is set to zero.
ERROR: Calculations terminating due to input errors.
Stopping.

Here is my current script:

Code: [Select]

DATABASE pitzer.dat

Title: Gas-solution equillibrium
Site: 1
Date: September 2023

SOLUTION 1
   units mmol/kgw
   Temp 25.0
   pH 11.8
   density 1.0
   Li 0.000541931
   Mg 0.000925735
   Si 0.000178031
   Ca 0.322832976

GAS_PHASE 1
    -fixed_pressure
    -pressure 1
    -volume 1
    -temperature 25
    CO2(g)    0.00042
    Ntg(g)     0.78
    O2(g)     0.2

USE SOLUTION 1
USE GAS_PHASE 1

END
Cheers,

M

[dlparkhurst]

Not sure what database you got, I don't think it was pitzer.dat. Use phreeqc.dat, or if you use pitzer.dat, you will need to use Oxg(g) for O2(g), which is a non-reactive form for O2. pitzer.dat does not have redox reactions.

[Mgill9460]

Hi Dr. Parkhurst

my code now executes. Although I still got an error that both Ntg(g) and Oxg(g) could not be identified in the phreeqc.dat.

My code currently is:

Code: [Select]

DATABASE phreeqc.dat

#Title: Gas-solution equillibrium
#Site: 1
#Date: September 2023

SOLUTION 1
   units mmol/kgw
   Temp 25.0
   pH 11.8
   density 1.0
   Li 0.000541931
   Mg 0.000925735
   Si 0.000178031
   Ca 0.322832976

SOLUTION 2
   units mmol/kgw
   Temp 25.0
   pH 13.0
   density 1.0
   Li 0.000541931
   Mg 0.000925735
   Si 0.000178031
   Ca 0.600000000

GAS_PHASE 1
    -fixed_pressure
    -pressure 1
    -volume 1
    -temperature 25
    CO2(g)    0.00042
    N2(g)     0.78
    O2(g)     0.2

KINETICS
USE SOLUTION 1
USE GAS_PHASE 1

END

KINETICS
USE SOLUTION 2
USE GAS_PHASE 1

END

So basically would this be reacting 1L of water with 1L of gas in a beaker? Is there a way to reduce the reaction time to one hour?

Thanks,

M

[dlparkhurst]

You must have an old version of PHREEQC. Download the latest version at https://www.usgs.gov/software/phreeqc-version-3.

[Mgill9460]

Hi Dr. Parkhurst,

I checked and I had the latest version of PHREEQC but I reinstalled just to be safe.

I did a bit of reading in the PHREEQC documentation and was able to modify the code to this stage:

Code: [Select]

#Title: Gas-solution equillibrium
#Site: 1
#Date: September 2023

#Solution 1
SOLUTION 1
    temp      25
    pH        11.98
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Br        0.002045287
    C         0.05965918
    Ca        1.004478766
    Cl        1.561218618
    Fe        4.79004e-06
    K         0.036419437
    Li        0.001239193
    Mg        0.002262909
    Na        0.99369291
    Si        0.000623108
    Sr        2.81443e-05
    Zn        1.00642e-05
#solution 2
SOLUTION 2
    temp      25
    pH        12.76
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Br        0.011740589
    C         0.062580487
    Ca        1.710374071
    Cl        9.89778016
    Fe        9.60516e-05
    K         0.443734015
    Li        0.010970821
    Mg        0.002454125
    Na        15.55023923
    Si        0.018297668
    Sr        0.000219202
    Zn        0.001718171

GAS_PHASE 1 headspace volume of gas over solution
    -fixed_volume
    -pressure 1
    -volume 1
    -temperature 25
    CO2(g)    0.00042
    Ntg(g)    0.78
    Oxg(g)    0.21

KINETICS 1
-steps       3600 in 1 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
USE SOLUTION 1
USE GAS_PHASE 1

END

KINETICS 2
-steps       3600 in 1 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
USE SOLUTION 2
USE GAS_PHASE 1

END
I know from experimental data that both solution 1 and solution 2 dissolve different amounts of carbon dioxide, however, in this simulation both of the KINETICS reactions seem to dissolve all the CO2 from the gas phase. Is there a potential problem with my setup?

Thank you,

M

[dlparkhurst]

First your initial solutions are way out of charge balance. In this script, I am assuming you added NaOH to adjust the pH.

Here a REACTION is used to add CO2 to the solution incrementally. You can think of it as kinetics without time. You can pick the endpoint of the reaction. If you have an infinite reservoir of air, it would be at P(CO2) 10^3.4. A fixed volume of gas would be lower.

You can work on converting REACTION to KINETICS if that is your intent.

Code: [Select]

#Title: Gas-solution equillibrium
#Site: 1
#Date: September 2023

#Solution 1
SOLUTION 1
    temp      25
    pH        11.98
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Br        0.002045287
    C         0.05965918
    Ca        1.004478766
    Cl        1.561218618
    Fe        4.79004e-06
    K         0.036419437
    Li        0.001239193
    Mg        0.002262909
    Na        0.99369291 charge
    Si        0.000623108
    Sr        2.81443e-05
    Zn        1.00642e-05
END
KNOBS
-iter 1000
USE solution 1
REACTION 1
CO2 1
0.02 mol in 50 steps
USER_GRAPH 1
    -headings               rxn P(CO2) pH
    -axis_titles            "CO2 added, moles" "P(CO2)" "pH"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X RXN
20 GRAPH_Y SI("CO2(g)")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true
END


[Mgill9460]

Hi Dr. Parkhurst,

Is it possible the reduction potential of the solutions could be affecting the charge balance? I have measured Eh values I could use but I was trying to keep the solutions as simple as possible for the initial setup.

When I run the code you provided I get a very large change in pH which is much larger than what we see experimentally.

I have played around a bit with the code but I cannot seem to replace the incremental reaction with a kinetic reaction in steps. Would it be something like the following though?

Code: [Select]

#Title: Gas-solution equillibrium
#Site: 1
#Date: September 2023

#Solution 1
SOLUTION 1
    temp      25
    pH        11.98
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Br        0.002045287
    C         0.05965918
    Ca        1.004478766
    Cl        1.561218618
    Fe        4.79004e-06
    K         0.036419437
    Li        0.001239193
    Mg        0.002262909
    Na        0.99369291 charge
    Si        0.000623108
    Sr        2.81443e-05
    Zn        1.00642e-05
END

GAS_PHASE 1 headspace volume of gas over solution
    -fixed_volume
    -pressure 1
    -volume 1
    -temperature 25
    CO2(g)    0.00042
    Ntg(g)    0.78
    Oxg(g)    0.2
END
REACTION 1
CO2 1
0.02 mol in 50 steps
END
KNOBS
-iter 1000
USE solution 1
KINETICS 1
-steps       3600 in 60 steps # seconds
-step_divide 60
-runge_kutta 3
-bad_step_max 500
USE SOLUTION 1
USE GAS_PHASE 1
USE REACTION 1
USER_GRAPH 1
    -headings               time P(CO2) pH
    -axis_titles            "Time, seconds" "P(CO2)" "pH"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X time
20 GRAPH_Y SI("CO2(g)")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true
END
Best,

M

[dlparkhurst]

The choice has a minimal effect on your charge balance. It only affects the distribution of Fe between Fe(2) and Fe(3) and C between C(4) and C(-4) in your solution. You are missing cations in your analysis or your anions are too large. Note that OH- is important at your pH.

Here is a script that kinetically approaches atmospheric CO2 in your solution. The rate is assumed proportional to the difference between the CO2(aq) concentration in equilibrium with the atmosphere and the CO2(aq) concentration in the solution at any point in time.

You can adjust the rate with the rate constant k in the RATES definition, or you can use a different rate equation entirely. I have arbitrarily set the rate to equilibrate in an hour or less.

(Note I added a definition in SOLUTION_SPECIES that makes the calculation more robust in this particular case.)

Code: [Select]

#Title: Gas-solution equillibrium
#Site: 1
#Date: September 2023

#Solution 1
SOLUTION_SPECIES
H2O + 0.01e- = H2O-0.01
log_k -8
SOLUTION 1
    temp      25
    pH        11.98
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Br        0.002045287
    C         0.05965918
    Ca        1.004478766
    Cl        1.561218618
    Fe        4.79004e-06
    K         0.036419437
    Li        0.001239193
    Mg        0.002262909
    Na        0.99369291 charge
    Si        0.000623108
    Sr        2.81443e-05
    Zn        1.00642e-05
END
RATES
CO2_ingassing
10 target_P = 10^-3.4
20 k = 1  # k is 1/s
30 target_m = target_P * 10^LK_PHASE("CO2(g)") / GAMMA("CO2")
40 rate = k*(target_m - MOL("CO2"))
50 moles = rate * TIME
60 SAVE moles
END
USE solution 1
KINETICS
CO2_ingassing
-formula CO2 1
-steps  3600 in 36
INCREMENTAL_REACTIONS
USER_GRAPH 1
    -headings               time P(CO2) pH
    -axis_titles            "Time, hours" "P(CO2)" "pH"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME /3600
20 GRAPH_Y SI("CO2(g)")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true
END


[Mgill9460]

Hi Dr. Parkhurst. I have been thinking about this and identified a rate at which CO2 will dissolve into solution.

I have added it to the code but now get an error that the mols of CO2 in solution are negative.

Also I was wondering if in addition to graphing the SI of CO2 if I also have the right code to graph the mols of CO2 in solution.

Code: [Select]

#Title: Gas-solution equillibrium
#Site: 1
#Date: September 2023

#Solution 1
SOLUTION_SPECIES
H2O + 0.01e- = H2O-0.01
log_k -8
SOLUTION 1
    temp      25
    pH        11.98
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Br        0.002045287
    C         0.05965918
    Ca        1.004478766
    Cl        1.561218618
    Fe        4.79004e-06
    K         0.036419437
    Li        0.001239193
    Mg        0.002262909
    Na        0.99369291 charge
    Si        0.000623108
    Sr        2.81443e-05
    Zn        1.00642e-05
    -water    1 # kg
END
RATES
CO2_ingassing
10 pH = 11.98
20 CO2 = (27.992*(pH^2))-(481.03*(pH))+(2116.1) #equation to define amount of CO2 in water
30 moles = CO2/60
40 SAVE moles
END
USE solution 1
KINETICS
CO2_ingassing
-formula CO2 1
-steps  3600 in 60 #1-hour in 60 second steps
INCREMENTAL_REACTIONS
USER_GRAPH 1
    -headings               time P(CO2) [CO2]
    -axis_titles            "Time, hours" "P(CO2)" "CO2 mols/L"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME /3600
20 GRAPH_Y SI("CO2(g)")
30 GRAPH_SY CO2
  -end
    -active                 true
END
Cheers,

[dlparkhurst]

First, I don't think your units are correct. You are adding 60 moles of CO2 at every invocation of the RATES, and rates is invoked several times per time step. So reconsider your units.

Second, you seem to be trying to approach a given concentration, so you should add a fraction of the difference between the desired concentration and the current concentration.