Reproducing hydrogeochemistry of carbonate and phyllosilicate assemblages

[mmelwani]

Hi all,

I'm trying to run a model by van Berk, W., & Fu, Y. (2011). (Reproducing hydrogeochemical conditions triggering the formation of carbonate and phyllosilicate alteration mineral assemblages on Mars (Nili Fossae region). Journal of Geophysical Research: Planets, 116(E10), E10006. https://doi.org/10.1029/2011JE003886).

The model in question is a simple batch reaction between water + CO2 + mafic rock in the first cell, followed by one-dimensional diffusive transport of the fluid into cells 2 - 50, which contain pure water in their pores. Their goal was to examine the resulting alteration mineralogy, especially the changing composition of carbonates, and quantify the total amounts of CO2 sequestered. Nothing particularly out of the ordinary, I should think.

Their whole input file is in their supplementary data (also here: https://agupubs.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1029%2F2011JE003886&file=jgre2976-sup-0002-ts01.doc).

However, when I create the PHREEQC input file and run it, I come across convergence issues. The first warnings are not consequential (Strengite is specified as a phase, but phosphorus is not included in the equilibrium phases or solution) and the model continues without a problem. I am confused about what might be causing convergence issues down the column, in the transport calculations. Any suggestions about what to try to get the model to run to completion?

Also, any thoughts about having pure water in equilibrium with the mafic rock in cells 2-50? Shouldn't the rock determine the composition of the fluid in the pores? Would it not be more appropriate to equilibrate pure water with the host rock (batch reaction), and then use the resulting fluid as the composition of the pore fluids? That would mean that the host rock does not consist of mafic minerals only, but also alteration minerals. I realized the model already equilibrates the pore fluids with the host rock for each of the cells.

Below is the code I used, translated to the best of my ability from the paper:

Code: [Select]

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\wateq4f.dat
TITLE 28MgO_pCO2=1.0_w/r=500_T=5
#------------------------------------------------------------- definition of additional phases (not included in wateq4f.dat)
PHASES
Enstatite
# G°f = -1459.0 kJ mol-1
# data from http://www.osti.gov/energycitations/product.biblio.jsp?osti_id=5274849
MgSiO3 + 2H+ + H2O = Mg+2 + H4SiO4
log_k           13.19

Fayalite
# G°f = -1379.4 kJ mol-1
# data from Stumm, W. & Morgan, J. J. Aquatic Chemistry, John Wiley, New York, 2nd edn., 1981.
Fe2SiO4 + 4H+ = 2Fe+2 + H4SiO4
log_k           16.647

Ferrosilite
# G°f = -1117,82 kJ mol-1
# data from Boyd T. D. & Scott S. D. Microbial and hydrothermal aspects of ferric oxyhydroxides
# and ferrosic hydroxides: the example of Franklin Seamount, Western Woodlark Basin, Papua New
# Guinea. Geochem. Trans. 7 (2001). DOI: 10.1039/b105277m.
FeSiO3 + 2H+ + H2O = Fe+2 + H4SiO4
log_k           7.11

Wollastonite
# log_k = 12.996
# from minteq.dat
# Parkhurst, D. L. & Appelo, C. A. J. User's guide to PHREEQC (Version 2) - A computer program
# for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. #U.S. Geological Survey Water-Resources Investigations Report 99-4259 (1999).
CaSiO3 + 2H+ + H2O = Ca+2 + H4SiO4
log_k           12.996
delta_h -19.498 kcal
#------------------------------------------------------------- definition of additional phases (not included in wateq4f.dat)

#------------------------------------------------------ pure water filling pore space of reactor no. 1 at starting conditions
SOLUTION 1                                                         
-pH 7.366
-density 1.000
-temp 5.0
-units mmol/kgw
#------------------------------------------------------ pure water filling pore space of reactor no. 1 at starting conditions

#------------------------------------------ unaltered rock composition [mol kgw-1] (reactor no. 1) at starting conditions
EQUILIBRIUM_PHASES 1
Enstatite         0.0 0.0378
Ferrosilite     0.0 0.0109
Wollastonite     0.0 0.0015
Forsterite     0.0 0.0364
Fayalite     0.0 0.0139
Albite     0.0 0.0020
Anorthite     0.0 0.0079
Pyrite     0.0 0.0005
CO2(g)                    0.0       # atmospheric pCO2 = 1.0 atm; amount of CO2 is 10.0 mol
#------------------------------------------ unaltered rock composition [mol kgw-1] (reactor no. 1) at starting conditions
#-------------------------------------------------------------------------- potential secondary mineral phases (reactor no. 1)
Al(OH)3(a) 0.0 0.0
Analcime 0.0 0.0
Anhydrite 0.0 0.0
Aragonite 0.0 0.0
Artinite 0.0 0.0
Basaluminite 0.0 0.0
Boehmite 0.0 0.0
Brucite 0.0 0.0
Calcite 0.0 0.0
#Chalcedony 0.0 0.0 # excluded
Chlorite14A 0.0 0.0
Chlorite7A 0.0 0.0
Chrysotile 0.0 0.0
Clinoenstatite 0.0 0.0
#Cristobalite 0.0 0.0 # excluded
Diaspore 0.0 0.0
#Diopside 0.0 0.0 # excluded
Dolomite 0.0 0.0
Dolomite(d) 0.0 0.0
Epsomite 0.0 0.0
Fe(OH)3(a) 0.0 0.0
Fe3(OH)8 0.0 0.0
FeS(ppt) 0.0 0.0
Gibbsite 0.0 0.0
Goethite 0.0 0.0
Greenalite 0.0 0.0
Greigite 0.0 0.0
Gypsum 0.0 0.0
Halloysite 0.0 0.0
Hematite 0.0 0.0
Huntite 0.0 0.0
Hydromagnesite 0.0 0.0
Jarosite-Na 0.0 0.0
JarositeH 0.0 0.0
Jurbanite 0.0 0.0
Kaolinite 0.0 0.0
Laumontite 0.0 0.0
Leonhardite 0.0 0.0
Mackinawite 0.0 0.0
Magadiite 0.0 0.0
Maghemite 0.0 0.0
Magnetite 0.0 0.0
Magnesite 0.0 0.0
Melanterite 0.0 0.0
Mirabilite 0.0 0.0
Montmorillonite-Ca 0.0 0.0
Nahcolite 0.0 0.0
Natron 0.0 0.0
Nesquehonite 0.0 0.0
Portlandite 0.0 0.0
Prehnite 0.0 0.0
#Pyrophyllite 0.0 0.0 # excluded
#Quartz 0.0 0.0 # excluded
Sepiolite 0.0 0.0
Sepiolite(d) 0.0 0.0
Siderite 0.0 0.0
Siderite(d)(3) 0.0 0.0
#Silicagel 0.0 0.0 # excluded
SiO2(a) 0.0 0.0
Strengite 0.0 0.0
Sulfur 0.0 0.0
Talc 0.0 0.0
Thenardite 0.0 0.0
Thermonatrite 0.0 0.0
#Tremolite 0.0 0.0 # excluded
Trona 0.0 0.0
Wairakite 0.0 0.0

SOLID_SOLUTIONS 1
TernarySolidSolution # not present at starting conditions
-comp Calcite 0.0 # not present at starting conditions
-comp Magnesite 0.0 # not present at starting conditions
-comp Siderite 0.0 # not present at starting conditions
#-------------------------------------------------------------------------- potential secondary mineral phases (reactor no. 1)
#------------------------------------------ unaltered rock composition [mol kgw-1] (reactor no. 1) at starting conditions

#---------------------------------------------- pure water filling pore space of reactors no. 2 to 50 at starting conditions
SOLUTION 2-50                                                       
-pH 7.366
-density 1.000
-temp 5.0
-units mmol/kgw
#------------------------------------- pure water filling pore space of reactors no. 2 to 50 at starting conditions

#---------------------------------------- unaltered rock composition (reactors no. 2 to 50) at starting conditions
EQUILIBRIUM_PHASES 2-50
Enstatite         0.0 0.0378
Ferrosilite     0.0 0.0109
Wollastonite     0.0 0.0015
Forsterite     0.0 0.0364
Fayalite     0.0 0.0139
Albite     0.0 0.0020
Anorthite     0.0 0.0079
Pyrite     0.0 0.0005
#------------------------------------------------- unaltered rock composition (reactors no. 2 to 50) at starting conditions
#------------------------------------------------------------------ potential secondary mineral phases (reactors no. 2 to 50)
Al(OH)3(a) 0.0 0.0
Analcime 0.0 0.0
Anhydrite 0.0 0.0
Aragonite 0.0 0.0
Artinite 0.0 0.0
Basaluminite 0.0 0.0
Boehmite 0.0 0.0
Brucite 0.0 0.0
Calcite 0.0 0.0
#Chalcedony 0.0 0.0 # excluded
Chlorite14A 0.0 0.0
Chlorite7A 0.0 0.0
Chrysotile 0.0 0.0
Clinoenstatite 0.0 0.0
#Cristobalite 0.0 0.0 # excluded
Diaspore 0.0 0.0
#Diopside 0.0 0.0 # excluded
Dolomite 0.0 0.0
Dolomite(d) 0.0 0.0
Epsomite 0.0 0.0
Fe(OH)3(a) 0.0 0.0
Fe3(OH)8 0.0 0.0
FeS(ppt) 0.0 0.0
Gibbsite 0.0 0.0
Goethite 0.0 0.0
Greenalite 0.0 0.0
Greigite 0.0 0.0
Gypsum 0.0 0.0
Halloysite 0.0 0.0
Hematite 0.0 0.0
Huntite 0.0 0.0
Hydromagnesite 0.0 0.0
Jarosite-Na 0.0 0.0
JarositeH 0.0 0.0
Jurbanite 0.0 0.0
Kaolinite 0.0 0.0
Laumontite 0.0 0.0
Leonhardite 0.0 0.0
Mackinawite 0.0 0.0
Magadiite 0.0 0.0
Maghemite 0.0 0.0
Magnetite 0.0 0.0
Magnesite 0.0 0.0
Melanterite 0.0 0.0
Mirabilite 0.0 0.0
Montmorillonite-Ca 0.0 0.0
Nahcolite 0.0 0.0
Natron 0.0 0.0
Nesquehonite 0.0 0.0
Portlandite 0.0 0.0
Prehnite 0.0 0.0
#Pyrophyllite 0.0 0.0 # excluded
#Quartz 0.0 0.0 # excluded
Sepiolite 0.0 0.0
Sepiolite(d) 0.0 0.0
Siderite 0.0 0.0
Siderite(d)(3) 0.0 0.0
#Silicagel 0.0 0.0 # excluded
SiO2(a) 0.0 0.0
Strengite 0.0 0.0
Sulfur 0.0 0.0
Talc 0.0 0.0
Thenardite 0.0 0.0
Thermonatrite 0.0 0.0
#Tremolite 0.0 0.0 # excluded
Trona 0.0 0.0
Wairakite 0.0 0.0

SOLID_SOLUTIONS 2-50
TernarySolidSolution # not present at starting conditions
-comp Calcite 0.0 # not present at starting conditions
-comp Magnesite 0.0 # not present at starting conditions
-comp Siderite 0.0 # not present at starting conditions
#---------------------------------------------------------------------- potential secondary mineral phases (reactor no. 2-50)

#----------------------------------------------------------------------------- one dimensional diffusive transport parameters
TRANSPORT
-cells 50
-lengths 50*2.0 # 50 * 2.0 meter
-shifts 100 # number of time steps
-time_step 3.1536e10 # seconds per shift; equivalent to 1000 a
-flow_direction diffusion_only
-diffusion_coefficient 1.0e-9 # effective diffusion coefficient [m2 s-1]
-boundary_conditions closed flux
-dump 1.dmp
-dump_frequency 1
-dump_restart 1
#----------------------------------------------------------------------------- one dimensional diffusive transport parameters
END


[dlparkhurst]

I think the solid-solutions calculation can be numerically difficult. It also does not help to have many, many phases to sort through for the equilibrium phases.

I was able to make the run by taking smaller steps in the Newton-Raphson iterations.

Code: [Select]

KNOBS
-step 10
-pe 5


[mmelwani]

Wonderful! Thank you, reducing the steps did indeed allow convergence. For future problems, I'll probably be satisfied with simple carbonate endmember compositions, rather than solid solutions.