MgSiO4 dissolution in column - transport model

[ms14218]

Dear contributor,
 
This is my first attempt at transport module in PHREEQC. I would appreciate any guidance. Note that i have tried to attach images to this post but it was not working. If required, I am happy to provide them through another means.
 
I am attempting to model transport of co2-saturated water into a 6cm length, column (0.6cm diam) filled with forsterite powder assuming porosity of 30%., for a period of 5 days at 70bar pressure, ROOM temp, 0.2ml/min flowrate. I use the llnl.dat.
I have produced the following scripts, using information from literature where relevant. Version 1 and 2 of the same script are shown below. I have also attached some of the graphical outputs from both scripts
 
1) In both versions of the script I am experiencing issues understanding the graphical outputs e.g. the conc vs time plot generated (i have tried to attach the image but it was not working), which shows the concentrations of desired elements returning to zero at the end of the simulation in script 1's case and at regular intervals before climbing up again in script 2's case. I am not sure which part of the code is causing this or rather what was coded incorrectly in my script. In short im not sure how to interpret these graphs.
 
2) I would like to simulate 5 days of concentration evolution in this model. which is why the first script transport block was set up in this manner (with shifts and time_step resulting in 5 days of total simulation time). However, i would like to incorporate a flowrate of 0.2ml/min or in my system 2.95e-05 m/s and given the length of one cell is set at 1e-03 I set a time_step of 33.9 s in the second version of the script.  This would require 12,743.4 shifts to result in 5 days (or 432000 seconds) of total simulation time in the model, the run time is therefore very long. Secondly the results still show the trend of resetting concentrations to zero at regular intervals.
 
3) I also would like to produce a graph of conc. vs pore volumes injected as shown in the second script. To do this I have estimated the total pore volumes injected in my experiment to be 707. Given that (shifts+0.5)/no.of cells is the PHREEQC formula for PV, and number of cells is 60, i would require 42000 shifts, and10 seconds per time_step, but this new configuration does not support the desired flowrate of the experiment.
 
Regarding points 2 and 3 i am unsure how to balance these inputs to create the simulation that best resembles my experiment.
 
Additional question: In this experiment, the outlet end of the column is connected to a back pressure regulator which maintains upstream system pressure @ 70 bars, I would like to ask is the boundary conditions shown are suited for this model to replicate such system where the outlet fluid is then collected for analysis.
 
 
SCRIPT 1:

Code: [Select]

SOLUTION 0
    temp      25
    pH        5.8
    pe        4
    redox     pe
    units     mmol/l
    density   1
    -water    1 # kg
 
EQUILIBRIUM_PHASES 0
    CO2(g) -3.398 # CO2 partial pressure at 70 bar (pre-equilibrated)
 
SOLUTION 1-60
    temp 25
    pH 3.0
    pressure 70
    units mol/kgw
    C(4) 0 charge  # CO2-saturated water
    water 1 # 1 kg of deionized water
 
SAVE solution 1-60
 
RATES
    Forsterite
-start
 10 REM acid solution parameters
 20 a1 = 0.54 # Rosso and Rimstidt (2000) rates data
 30 e1 = 42600
 40 n1 = 0.5
 50 REM neutral solution parameters
 60 a2 = 1614
 70 e2 = 79000
 80 REM basic solution parameters
 90 a3 = 0
100 e3 = 0
110 n2 = 0
120 REM rate=0 if no mineral mass or undersaturated
130 IF(M <= 0) THEN GOTO 240
140 IF(SI("Forsterite") >= 0) THEN GOTO 240
150 IF(M0 <= 0) THEN sa = PARM(1) # Surface area in m?
160 IF(sa <= 0) THEN sa = 1
170 r = 8.31451
180 rate1 = a1*ACT("H+")^n1*EXP(-e1/r/TK)
190 rate2 = a2*EXP(-e2/r/TK)
200 rate3 = a3*ACT("OH-")^n2*EXP(-e3/r/TK)
210 rate = (rate1 + rate2 + rate3) * (1 - SI("Forsterite")) * sa * PARM(2)
220 moles = rate * TIME
230 IF(moles > M) THEN moles = M
240 SAVE moles
-end
 
KINETICS 1
Forsterite
    -formula  Mg2SiO4  1
    -m        0.0284 # Moles of forsterite (we state 4g as an approx. quantity then do 4g / 140.69g/mol)
    -m0       0.0284 # Initial moles
    -parms    1.62 1 # Parm 1 is surface area (we take BET measurement m2/g then multiply by approx. mass in g)
    -tol      1e-08  # a default value
-steps       1
-step_divide 1
-runge_kutta 3
-bad_step_max 500
 
REACTION_TEMPERATURE 1
    25
 
TRANSPORT
    -cells         60    # 6 cm or 0.06 m(i.e. mesoscale dimension) divided into 60 cells
    -length        1e-03 # Each cell is 1 mm in length
    -time_step     86400 # 1 day per shift in seconds
    -shifts        5     # 5 days of simulation
    -flow_direction forward
    -punch_cells   1-60  # Record data for all cells
    -punch_frequency 5
    -boundary_conditions flux constant # Pressure-maintained system
    -porosity      0.3   # Porosity of 30%
    -dispersivity  0.01  # Longitudinal dispersivity default value used
    -correct_disp  true
    -diffusion coefficient 0.78e-09
    -thermal diffusion coefficient 0.78e-09
    -warnings              false
 
#EXCHANGE 1
    #X 0.001  # lower total exchange capacity for Mg2+
 
#EXCHANGE_SPECIES
    #MgX2 = Mg+2 + 2X-
        #log_k 0.5  # Selectivity for Mg2+
 
#USE exchange 1
 
KNOBS
    -iterations 500  # Allow more iterations per time step
    -tolerance 1e-8  # Loosen the tolerance for convergence
 
USER_GRAPH 1
    -headings               Time Mg Si
    -axis_titles            "Time (days)" "Concentration (mmol/L)" ""
    -chart_title            "Concentrations with reaction time"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 IF (STEP_NO = 1) THEN CUM_MG = 0 #maintaining the cumulative values across shifts, to correctly track concentration changes over time without resetting them each day.
20 IF (STEP_NO = 1) THEN CUM_SI = 0
30 CUM_MG = CUM_MG + TOT("Mg")  # Accumulate Mg over time
40 CUM_SI = CUM_SI + TOT("Si")  # Accumulate Si over time
50 graph_x TOTAL_TIME / 86400  # Convert seconds to days
60 graph_y CUM_MG
70 graph_y CUM_SI
#10 graph_x TOTAL_TIME/86400
#20 graph_y TOT("Mg")
#30 graph_y TOT("Si")
  -end
 
USER_GRAPH 2
    -headings               Distance Mg Si
    -axis_titles            "Distance (m)" "Concentration, mmol/L" ""
    -chart_title            "Concentrations with distance across beadpack"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 IF (STEP_NO = 1) THEN CUM_MG = 0 #maintaining the cumulative values across shifts, to correctly track concentration changes over time without resetting them each day.
20 IF (STEP_NO = 1) THEN CUM_SI = 0
30 CUM_MG = CUM_MG + TOT("Mg")  # Accumulate Mg over time
40 CUM_SI = CUM_SI + TOT("Si")  # Accumulate Si over time
50 graph_x dist 
60 graph_y CUM_MG
70 graph_y CUM_SI
#10 graph_x dist
#20 graph_y TOT("Mg")
#30 graph_y TOT("Si")
  -end
 
SELECTED_OUTPUT
    -file forsterite_beadpack_output.txt
    -ph
    -totals Mg Si
 
END
 
 
SCRIPT 2:

Code: [Select]

SOLUTION 0
    temp      25
    pH        5.8
    pe        4
    redox     pe
    units     mmol/l
    density   1
    -water    1 # kg
 
EQUILIBRIUM_PHASES 0
    CO2(g) -3.398 # CO2 partial pressure at 70 bar (pre-equilibrated)
 
SOLUTION 1-60
    temp 25
    pH 3.0
    pressure 70
    units mol/kgw
    C(4) 0 charge  # CO2-saturated water
    water 1 # 1 kg of deionized water
 
SAVE solution 1-60
 
RATES
    Forsterite
-start
 10 REM acid solution parameters
 20 a1 = 0.54 # Rosso and Rimstidt (2000) rates data
 30 e1 = 42600
 40 n1 = 0.5
 50 REM neutral solution parameters
 60 a2 = 1614
 70 e2 = 79000
 80 REM basic solution parameters
 90 a3 = 0
100 e3 = 0
110 n2 = 0
120 REM rate=0 if no mineral mass or undersaturated
130 IF(M <= 0) THEN GOTO 240
140 IF(SI("Forsterite") >= 0) THEN GOTO 240
150 IF(M0 <= 0) THEN sa = PARM(1) # Surface area in m?
160 IF(sa <= 0) THEN sa = 1
170 r = 8.31451
180 rate1 = a1*ACT("H+")^n1*EXP(-e1/r/TK)
190 rate2 = a2*EXP(-e2/r/TK)
200 rate3 = a3*ACT("OH-")^n2*EXP(-e3/r/TK)
210 rate = (rate1 + rate2 + rate3) * (1 - SI("Forsterite")) * sa * PARM(2)
220 moles = rate * TIME
230 IF(moles > M) THEN moles = M
240 SAVE moles
-end
 
KINETICS 1
Forsterite
    -formula  Mg2SiO4  1
    -m        0.0284 # Moles of forsterite (we state 4g as an approx. quantity then do 4g / 140.69g/mol)
    -m0       0.0284 # Initial moles
    -parms    1.62 1 # Parm 1 is surface area in m?(we take BET measurement m2/g then multiply by approx. mass in g)
    -tol      1e-08  # a default value
-steps       1
-step_divide 1
-runge_kutta 3
-bad_step_max 500
 
REACTION_TEMPERATURE 1
    25
 
TRANSPORT
    -cells                 60
    -shifts                12743
    -time_step             33.9
    -boundary_conditions   flux constant
    -lengths               0.001
    -dispersivities        0.01
    -correct_disp          true
    -diffusion_coefficient 7.8e-10
    -thermal_diffusion     2   7.8e-10
    -punch_frequency       60
    -warnings              false
 
USER_GRAPH 1
    -headings               Time Mg Si
    -axis_titles            "Time (days)" "Concentration (mmol/L)" ""
    -chart_title            "Concentrations with reaction time"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 IF (STEP_NO = 1) THEN CUM_MG = 0 #maintaining the cumulative values across shifts, to correctly track concentration changes over time without resetting them each day.
20 IF (STEP_NO = 1) THEN CUM_SI = 0
30 CUM_MG = CUM_MG + TOT("Mg")  # Accumulate Mg over time
40 CUM_SI = CUM_SI + TOT("Si")  # Accumulate Si over time
50 graph_x TOTAL_TIME / 86400  # Convert seconds to days
60 graph_y CUM_MG
70 graph_y CUM_SI
#10 graph_x TOTAL_TIME/86400
#20 graph_y TOT("Mg")
#30 graph_y TOT("Si")
  -end
 
USER_GRAPH 2
    -headings               Distance Mg Si
    -axis_titles            "Distance (m)" "Concentration, mmol/L" ""
    -chart_title            "Concentrations with distance across beadpack"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 IF (STEP_NO = 1) THEN CUM_MG = 0 #maintaining the cumulative values across shifts, to correctly track concentration changes over time without resetting them each day.
20 IF (STEP_NO = 1) THEN CUM_SI = 0
30 CUM_MG = CUM_MG + TOT("Mg")  # Accumulate Mg over time
40 CUM_SI = CUM_SI + TOT("Si")  # Accumulate Si over time
50 graph_x dist 
60 graph_y CUM_MG
70 graph_y CUM_SI
#10 graph_x dist
#20 graph_y TOT("Mg")
#30 graph_y TOT("Si")
  -end
 
USER_GRAPH 3
    -headings               PV Mg Si
    -axis_titles            "Pore volumes injected" "Concentration, mmol/L" ""
    -chart_title            "Concentrations against injected poe volumes"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 IF (STEP_NO = 1) THEN CUM_MG = 0 #maintaining the cumulative values across shifts, to correctly track concentration changes over time without resetting them each day.
20 IF (STEP_NO = 1) THEN CUM_SI = 0
30 CUM_MG = CUM_MG + TOT("Mg")  # Accumulate Mg over time
40 CUM_SI = CUM_SI + TOT("Si")  # Accumulate Si over time
50 graph_x (step_no + 0.5)/cell_no  #calculation for pore volume to be plotted
60 graph_y CUM_MG
70 graph_y CUM_SI
#10 graph_x (step_no + 0.5)/cell_no 
#20 graph_y TOT("Mg")
#30 graph_y TOT("Si")
  -end
 
SELECTED_OUTPUT 1
    -file                 forsterite_beadpack_output_v2.txt
    -simulation           true
    -solution             true
    -pH                   true
    -temperature          true
    -totals               Mg  Si
    -saturation_indices   SiO2(am)  Magnesite
    -kinetic_reactants    Forsterite
 
END

[dlparkhurst]

Study the following script. I initially considered 0.001 m cell lengths, and then adjusted to 0.01 m cell lengths. The model with bigger cells runs quickly, and you can investigate whether it resembles your experiment. (The larger the dispersivity, the slower will be the run times, so you may want a small dispersivity until you are happy with the simulation.)

I calculated a water volume of 0.85 mL in a 0.001 m x pi*.3^2 cell, resulting in 255 s to transit the cell at a flow rate of 0.2 mL/min.

The rates in the cells are faster initially, when there is CO2 in each cell. Rates slow down when the only source of CO2 is the inflowing solution.

Code: [Select]

RATES
    Forsterite
-start
 10 REM acid solution parameters
 20 a1 = 0.54 # Rosso and Rimstidt (2000) rates data
 30 e1 = 42600
 40 n1 = 0.5
 50 REM neutral solution parameters
 60 a2 = 1614
 70 e2 = 79000
 80 REM basic solution parameters
 90 a3 = 0
100 e3 = 0
110 n2 = 0
120 REM rate=0 if no mineral mass or undersaturated
130 IF(M <= 0) THEN GOTO 240
140 IF(SI("Forsterite") >= 0) THEN GOTO 240
150 IF(M0 <= 0) THEN sa = PARM(1) # Surface area in m?
160 IF(sa <= 0) THEN sa = 1
170 r = 8.31451
180 rate1 = a1*ACT("H+")^n1*EXP(-e1/r/TK)
190 rate2 = a2*EXP(-e2/r/TK)
200 rate3 = a3*ACT("OH-")^n2*EXP(-e3/r/TK)
210 rate = (rate1 + rate2 + rate3) * (1 - SI("Forsterite")) * sa * PARM(2)
220 moles = rate * TIME
230 IF(moles > M) THEN moles = M
240 SAVE moles
-end
END
SOLUTION 0
    temp      25
    pH        5.8 charge
    pe        4
    redox     pe
    units     mmol/l
    density   1 calc
    C(4) 1 CO2(g) -3.398
    -water    1 # kg
END

SOLUTION 1-60
    temp 25
    pH 3.0
    pressure 70
    units mol/kgw
    C(4) 1 CO2(g) -3.398  # CO2-saturated water
    Cl  1 charge
    water 1 # 1 kg of deionized water
END
 
KINETICS 1-60
Forsterite
    -formula  Mg2SiO4  1
    -m        0.0284 # Moles of forsterite (we state 4g as an approx. quantity then do 4g / 140.69g/mol)
    -m0       0.0284 # Initial moles
    -parms    1.62 1 # Parm 1 is surface area (we take BET measurement m2/g then multiply by approx. mass in g)
    -tol      1e-08  # a default value
-steps       1
-step_divide 1
-runge_kutta 3
-bad_step_max 500
END
REACTION_TEMPERATURE 1-60
    25
REACTION_PRESSURE 1-60
    70 # atm
END
TRANSPORT
    -cells         6 # 60    # 6 cm or 0.06 m(i.e. mesoscale dimension) divided into 60 cells
    -length        0.01 # 1e-03 # Each cell is 1 mm in length
    -time_step     2550  # 255   # s/cell advective velocity 0.01/2550 = 4e-6 m/s
    -shifts        170 # 60    # days of simulation
    -flow_direction forward
    -punch_cells   1-6 # 1-60  # Record data for all cells
    -punch_frequency 1
    -boundary_conditions flux flux # Pressure-maintained system
    -porosity      0.3   # Porosity of 30%
    -dispersivity  0.01  # Longitudinal dispersivity default value used
    -correct_disp  true
    -diffusion coefficient 0.78e-09
    -thermal diffusion coefficient 0.78e-09
    -warnings              false
 
USER_GRAPH 1
    -headings               Time Mg Si Fosterite_dis
    -axis_titles            "Distance, m" "Concentration (mol/kgw)" "Fosterite dissolved, moles"
    -chart_title            "Concentrations with distance"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
#5 if (step_no <> 60) then goto 100
5 if (step_no <> 170) then goto 100
10 graph_x dist
20 graph_y TOT("Mg")
30 graph_y TOT("Si")
40 graph_sy M0 - M
100 END
  -end

USER_GRAPH 2
    -headings               Time Cummulative_dissolved Mg Si
    -axis_titles            "Time, days" "Fosterite dissolved, moles" "mol/kgw, end of column"
    -chart_title            "Concentrations with reaction time"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
30 PUT(M0 - M, cell_no)
40 if (cell_no <> 60) then goto 200
40 if (cell_no <> 6) then goto 200
45 sum = 0
50 FOR i = 1 TO 60
60 sum = sum + GET(i)
70 NEXT i
80 graph_x TOTAL_TIME / 86400
90 graph_y sum
100 graph_sy TOT("Mg")
110 graph_sy TOT("Si")
200 end
  -end
 
END

[dlparkhurst]

I should have mentioned that I used the database phreeqc_rates.dat. Like phreeqc.dat, Amm.dat, and pitzer.dat, phreeqc_rates.dat has the molar volumes of aqueous species to do calculations at pressure. Note that many of the phases lack molar volumes (forsterite, for example), but you can add them if you know the molecular weight and density of the mineral (v = mw/rho).

Note that llnl.dat cannot vary pressure and temperature independently. Its equilibrium constants are based on 1 atm pressure up to 100C, and then the steam curve above 100C.

[ms14218]

Many thanks for this detailed response and for the corrections made to the script. I am understanding the transport module more now.

I would like to ask where the phreeqc_rates.dat can be found? When clicking the options tab -> Set default database, I do not find it in the available databases list.

Regarding the molar volume of forsterite. I have a value in cm^3 calculated from average density value for the mineral. However, am unsure how this can be added into the script.

Lastly, since forsterite (Mg2SiO4) is not in the phreeqc_rates.dat, would it be possible to instead calculate the molar volume of a forsterite mineral phase with stoichiometry matching my experiment sample i.e. with Fe content?

[dlparkhurst]

You need to get the latest version 3.8.2 of PHREEQC from https://github.com/usgs-coupled/. Phreeqc_rates.dat will be in the database directory, or it will be available if you use PhreeqcI in the run dialog.

A definition for Forsterite is available in phreeqc_rate.dat. I would include the following in your input file. You can make changes as you want.

Code: [Select]

PHASES
Forsterite
Mg2SiO4 + 4 H+ = 2 Mg+2 + H4SiO4
log_k 28.609
delta_h -217.115 #kJ/mol #Internal calculation
-analytic -1.0983766E3 -1.5385695E-1 7.321503E4 3.91599E2 -3.7061609E6
        -Vm ?? #cm3/mol = MW / rho

Finally, the Palandri&Kharaka parameters are available in the data block RATE_PARAMETERS_PK in phreeqc_rates.dat. Examples of using the Basic function RATE_PK are available in comments in phreeqc_rates.dat, and both the data block and the Basic function are documented in documentation in PhreeqcI and in the doc directory of the installation. So, check out the documentation.

[ms14218]

Thank you for your guidance. I have run successfully in v3.8.2. I have some minor follow up questions:

1) What do the values in the analytical row of the PHASES block correspond to?

2)Re. the plotting of the pore volumes injected as i have in my original scripts, would this still fit in the current version of the script i.e. is the code I have for plotting PV correct?

[dlparkhurst]

Check the manual for the analytical expression coefficients. It is the same for all log Ks (SOLUTION_SPECIES, PHASES, ...).

Put a tracer in to see if the breakthrough occurs when you think it should.

[ms14218]

Many thanks for your responses.

I am wondering if it is possible to run this simulation, allowing the SI to be calculated for different phases, and then to remove the CO2 content + depressurise the system (resulting in a new pH) and re-calculate the SI to be able to see how the SI of phases dissolved in solution evolves from the reaction stage to the depressurisation stage?

The issue i am having is how to maintain the output from the first simulation and essentially utilising it as an input for the next simulation where CO2 and pressure are removed (back to ambient conditions).

[dlparkhurst]

Sounds like you need the SAVE and USE keywords.

[ms14218]

Hi,
I have applied the following code after the SELECTED_OUPUT from my transport model:

Code: [Select]

SAVE solution 4  # Save the final solution from transport model after the first step

# modify the pressure and CO2 after depressurization:
SOLUTION 5
    temp 25
    pH 3.0
    pressure 1 # Set the pressure to 1 atm (atmospheric)
    units mol/kgw
    C(4) 0 # Remove CO2(g) from the system
    water 1 # 1 kg of deionized water
END

EQUILIBRIUM_PHASES 1
    CO2(g) 0 # Set CO2 partial pressure to zero after depressurization

END

REACTION_PRESSURE 1
    1 # Set the pressure to 1 atm after depressurization

USE solution 5  # Use the saved solution for the next calculation
USE equilibrium_phases 1  # Apply the new equilibrium phases

It corrects to atmospheric conditions, however the SI indices wrt to different mineral phases are not calculated at the end of this step as they were at the end of the reactive transport simulation.  That data is not retained it seems for the final re-calculation after solution is depressurised and reset to atm conditions.

[dlparkhurst]

Well, no, solution 5 only has C and no other solutes, so there would be no other saturation indices.

I think you want to USE solution 4 with a lower CO2 pressure in equilibrium phases. Note that the first parameter after CO2(g) is the log of the partial pressure, so your script specifies 1 atm. If, for example, you make that -3.4, it will remove CO2 to the point of the atmospheric partial pressure of CO2.

Code: [Select]

SAVE solution 4  # Save the final solution from transport model after the first step
END
USE solution 4
EQUILIBRIUM_PHASES 1
    CO2(g) -3.4 # Set CO2 partial pressure to zero after depressurization
REACTION_PRESSURE 1
    1 # Set the pressure to 1 atm after depressurization
END


[ms14218]

Thank you for your quick response. I had initially thought to USE solution 4 in a similar, albeit not the exact way which you have shown, however I incorrectly opted to create a new solution (5) after saving solution 4. Thank you for pointing this mistake out and for the explanation.