Reactive Transport Modeling of minerals

[chrisnguen38]

Hi,

I would like to model the reactions of shale with CO2-saturated (sodium chloride) brine following the work described in the paper "A siltstone reaction front related to CO2- and sulfur-bearing fluids: Integrating quantitative elemental mapping with reactive transport modeling" by Maskell et al. (2018).
I want to plot mineral concentrations vs distance for different minerals (calcite, clay, quartz, and felspar). I have tried plotting the latter, but I have been struggling to do so. I only obtained element concentration vs. distance for different cells and the results were not great.
I have started with the code below:
TITLE Example. --
SOLUTION 1  PURE WATER
        pH      7.0
        temp    100
   -water      0.055784   # kg
EQUILIBRIUM_PHASES 1 Define amounts of phases in assemblage.
        Halite          0.0     0.055886379
SAVE SOLUTION 1
END
SOLUTION 0 Acidic water
USE solution 1
EQUILIBRIUM_PHASES 1 Equilibrate rock with carbon dioxide
CO2(g)      2.08807475135
SAVE solution 0
END
SOLUTION 1-20
USE solution 1
EQUILIBRIUM_PHASES 1 Equilibrate rock with pore water
   Kaolinite   0.0    0.000316662
   Quartz      0.0    0.001036349
        Calcite     0.0    0.027175368
   Dolomite    0.0    0.000143937
   K-Feldspar  0.0    0.000259087   
END
PRINT
        -selected_out true
KINETICS 1-20 Kinetic reaction for all cells # 1-20
Calcite
      -m   2.713e-02
      -parms   5.0      0.3
      -tol     1.e-8
TRANSPORT
        -cells                20               
        -length               0.005               
        -shifts               336               
        -time_step            3600             
        -flow_direction       forward
        -boundary_condition   flux flux
        -dispersivity         .05
        -correct_disp         true
        -diffusion_coef       0.0e-9
        -punch_cells          20               
        -punch_frequency      1                 
        -print_cells          20               
        -print_frequency      5                 
END
PRINT
    -selected_out true
USER_GRAPH 1 Plots Ca concentration profile
      -axis_titles "Distance, in m"
              "CO2 concentration, in moles per kilogram water"
       -axis_scale x_axis auto
       -axis_scale y_axis auto
           -connect_simulations true
           -plot_concentration_vs  x
END
My question is: how do I plot mineral concentrations vs. distance for different cells?

Thank you for your time. Your help would be much appreciated!

[dlparkhurst]

Look at the section of the manual on the Basic Interpreter for Basic functions that return the amounts of equilibrium phases and kinetic reactants.

Code: [Select]

TRANSPORT
        -cells                20               
        -length               0.005               
        -shifts               5               
        -time_step            3600             
        -flow_direction       forward
        -boundary_condition   flux flux
        -dispersivity         .05
        -correct_disp         true
        -diffusion_coef       0.0e-9
        -punch_cells          1-20               
        -punch_frequency      5                 
        -print_cells          20               
        -print_frequency      5                 
USER_GRAPH 1 Plots Ca concentration profile
    -axis_titles            "Distance, in m" "CO2 concentration, in moles per kilogram water" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X DIST
20 GRAPH_Y MOL("CO2")
  -end
    -active                 true
END


[Yongqiang]

Dear Parkhurst,

If we want to plot the concentration of calcite left in the cells, can we use the MOL("Calcite") in the USER_GRAPH section?

Regards,
Yongqiang

[Yongqiang]

I tried the MOL("Calcite"), but it gives a value of 0.

Regards,
Yongqiang

[dlparkhurst]

There is a very nice section of the manual describing Basic functions. Many of the items could also be defined by SELECTED_OUTPUT options, also described in the manual.

EQUI("Calcite") Moles of a phase in the equilibrium-phase assemblage.
 
MOL("HCO3-") Molality of an aqueous, exchange, or surface species.

See also SELECTED_OUTPUT; -equilibrium_phases.
 

[Yongqiang]

Yes, thank you, Professor PARKHURST. I also see the EQUI_DELTA("Calcite"), which is very useful. Thank you very much!!

[sanem]

Hello,

I have been trying to simulate CO2 injection into a reservoir using Reactive Transport Modelling. First of all, I tried to equilibrate the injected CO2 with a Solution 0. The mol value of CO2 is calculated from the CO2 solubility for Solution 1 at specific T and P values. Then, I added the reservoir mineralogy as “Equilibrium_Phases” and ran the code; however, when I try to examine the change in the amount of minerals with the increasing distance, there seems to be no change at all (except for the first reaction).

Code: [Select]


SOLUTION 0
    temp      25
    pH        7
    density   1
    -water    1 # kg
EQUILIBRIUM_PHASES 0
    CO2(g)    2.01 1.179

USE solution 0
USE equilibrium_phase 0
SAVE solution 0
END

SOLUTION 1-10
    temp      62
    pH        7.76
    pe        4
    redox     pe
    units     mg/l
    density   1
    Al        0.066
    C(4)      54
    Ca        25
    Cl(-1)    48
    Fe        0.058
    K         3.4
    Mg        0.2
    N(5)      0.18
    Na        193
    S(6)      360
    Si        32
    -water    1 # kg

EQUILIBRIUM_PHASES 1-10 EDR-1032
    Albite    0 49.664
    Anorthite 0 33.807
    Calcite   0 5.22
    Diopside  0 4.826
    Epidote   0 0
    K-feldspar 0 11.264
    Muscovite 0 0
    Phlogopite 0 11.216
    Quartz    0 260.902
    Tremolite 0 2.573
    Magnesite 0 0

TRANSPORT
    -cells                 10
    -shifts                10
    -time_step             3154000000 # seconds 100yrs
    -dispersivities        10*1
    -correct_disp          true

I also tried using “Kinetics” data block for the reservoir mineralogy, but then the simulation is stuck at Transport step 1 (my input file for the simulation with Kinetics is attached).   

I would like to understand whether I am using a wrong or insufficient input code for such simulations.

I would appreciate any help and/or advice.
Thank you in advance!

[dlparkhurst]

I did not look carefully, but I suspect when you used EQUILIBRIUM_PHASES, the solution reached equilibrium in the first step, and then simply remained at equilibrium, so there was no additional change in composition. Note that the calculation finds the stable phase assemblage such that some of your minerals dissolve entirely and the components precipitate as other minerals. You are approaching the phase rule, where the number of minerals that can coexist for a given set of components is limited.

With kinetics, the simulation is getting hung up on very small concentrations of Fe. If you remove Epidote, the calculation would proceed farther. That is not really a good answer, but I think in general you should start with a simpler set of reactions. Start with calcite, it is most like ly to react quickly. If some of the reactants approach equilibrium, it is a simpler calculation to set some of the minerals as equilibrium conditions and some as kinetic. For example, for longer simulations, assume calcite is in equilibrium.

Finally, relying on the relative rates of reaction of 10 minerals results in a very uncertain calculation. Unless you can simplify or put some constraints on the calculation, it is hard to put much faith in the result.

[sanem]

Thank you very much for the response! I have tried to simplify my calculations (moved some of the minerals from kinetics to equilibrium phase), and I was able to get results. However, the model is working as if CO2 does not move through cells and only reacting in the first cell. Is it because a data block/step is missing from my input code, or is it simply because of the long time step that I defined?

I would like to ask another question, if it is not too much trouble. What is the relationship between the time step in the “Kinetics” data block and the time step in the “Transport” data block when used in the same simulation? I defined a time step in my Kinetics data block previously, and when I deleted that part (thinking that I might be doing wrong), the program gave me the following error.

ERROR: Bad RK steps > 500. Please decrease (time)step or increase -bad_step_max.

I managed to get results with 30 minutes defined as time step in Transport (it was 100 years previously) and these results show changes in every cell (although the changes are really minor). Would it be also better to use shorter time steps for such simulations?

Best regards.

[dlparkhurst]

You probably need to understand the sequence of calculations in PHREEQC. If you define

SOLUTION
...
KINETICS
...
TRANSPORT
...
END

PHREEQC will perform 3 calculations. Looking at the output file, there will be an initial solution calculation, a reaction calculation of solution with kinetics, and a transport calculation. The time step for kinetics in the reaction calculation will be the one defined in KINETICS or 1 second by default. The time step for kinetics in TRANSPORT will be the one defined in TRANSPORT.

If you define the following

SOLUTION
END
KINETICS
END
TRANSPORT
END

there will be only two calculations, an initial solution calculation and a transport calculation. The time step used in the transport calculation is always the one defined in the TRANSPORT data block.

As for CO2 not moving through the system, perhaps the rate is fast (or other equilibrium phases consume the CO2) and most reaction occurs in the first cell.

[sanem]

Thank you so much David! I really appreciate your helpful comments.

Regards.