inverse model, silicate weathering, groundwater recharge

[Brendan]

Hello,

I am trying to create some inverse models to help understand pathways of mountain block recharge in volcanic rocks. I have plenty of analytical data to choose from. Here, I am hoping to see chemical evolution between a mountain spring (solution 1) and a monitoring well screen in dacite (solution 2). For the spring water I have temp, pH, pe, alkalinity, Al, Ba, Ca, Cl, F, Fe, Mg, Mn, PO4 as P, K, SiO2, Na and SO4. I have all that and more trace elements for the groundwater sample. Each of these has a small charge balance error when I run them as speciation models. So then I want to put them together in an inverse model, specifying that plagioclase compositions, along with some pyroxenes and hornblende should be subject to dissolution, based on the SI's and that they're present. Also I would specify CO2 dissolving. Also based on SI's, I would expect to see precipitation of maybe smectite and montmorillonite. I would not expect quartz or microcline to participate. I'm using llnl.dat.
First, I notice that the analytes I have that aren't in the phases are in the way - I have to list them as -balances, but then they're still often either not satisfying equality or inequality (errors), which I'm not sure I understand. I want those concentrations in my model for charge balance reasons, I think, but how do I go about scaling the uncertainties in -balances? I was able to get models by keeping a few such as Cl, K, Fe and Li, and putting a 1 after each of them, but I believe the models didn't make sense. They seemed to precipitate both kaolinite and CO2. I removed CO2(g) from the phases then. In each case, and in every model, the mass transfers are greater for albite than anorthite. Since albite reaches saturation somewhere along the flow path, I would want to see this reversed, so I think I need ion exchange there...? Sodium is kept in solution while calcium precipitates in clays?
I think perhaps I should use a mass balance equation for plagioclase, as in the example. The rocks have between An44 and An50. If I do, I'll have to adapt from H4SiO4 to SiO2, as dlparkhurst pointed out in one of these threads, because I'm using llnl.dat. I've played around with this to try to see what has what effect, here is where I was last.
Many thanks to anyone who can help me.
-Brendan

TITLE   Inverse model of regional recharge in Tschicoma Fm

SOLUTION 1   #Spring 5.29 water
   temp   4.8
   units   mg/L
   pH   5.6
   pe   6.05   #ORP is 333mV, using Nernst eqn
   Alkalinity   40.0   
   Al   2.07
   Ca   6.72
   Cl   1.6
   Fe(+3)   .703
   Mg   2.61
   K   1.49
   Si   34.4 as SiO2
   Na   3.76
SOLUTION 2   #R-37-2_S4_c
   temp   21.96
   units   mg/L
   pe   1.4   #ORP was 82.3 mV
   pH   7.64
   Alkalinity   49.3
   Cl   1.58
   Ca   9.1
   Mg   2.78
   Fe(+3)   .015
   K   1.26
   Na   10.6
   Al   .0029997
   Li   .0436835
   Si   61.3786 as SiO2
INVERSE_MODELING
   -solution 1 2
   -uncertainties 0.035 0.035
   -balances
      Cl   1
      Fe(+3)   1
      K   1
      Li   1
      C(4)   1
   -phases
      Kaolinite
      Montmor-Mg      precip
      Ferrosilite         dissolve
      Enstatite         dissolve
      Tremolite         dissolve
      Albite         dissolve
      Anorthite         dissolve
      Spodumene      dissolve
END

[dlparkhurst]

First, inverse modeling requires knowledge of the stoichiometry of the reactive minerals. In general, the stoichiometry of the trace elements is not well known, and trace elements can be affected by sorption to mineral surfaces. So I think inverse modeling works best for the major elements, including the major elements in the minerals, but skipping the trace elements.

If elements are trace, then removing them from the calculation should not affect the charge balance very much. Your waters have a positive imbalance, but you also lack sulfate analyses. If you have sulfate analyses, you should include it in the solution definition (could be an indicator of pyrite dissolution). Without the sulfate, you need to have a charge imbalance of about 15 percent for the inverse modeling, which is not great, but perhaps will allow some useful interpretations. Inverse modeling is far less quantitative than we like to think. It is really just testing that a set of reactions is plausible.

The following gives four minimal models that you can consider. I think you should probably not use albite, but a plag formula instead, or perhaps anorthite if you think only the Ca-rich parts of the feldspars are dissolved. I don't like having just Montmor-Mg in that the clay at least would have Ca, and probably some Na. However, if you add exchange for all the cations, to allow a variable clay composition, you may end up with way too many models to be useful.

TITLE   Inverse model of regional recharge in Tschicoma Fm

SOLUTION 1
    temp      4.8
    pH        5.6
    pe        6.05
    redox     pe
    units     mg/l
    density   1
    Al        2.07
    Alkalinity 40
    Ca        6.72
    Cl        1.6
    Fe(3)     0.703
    K         1.49
    Mg        2.61
    Na        3.76
    Si        34.4 as SiO2
    -water    1 # kg
SOLUTION 2
    temp      21.96
    pH        7.64
    pe        1.4
    redox     pe
    units     mg/l
    density   1
    Al        0.0029997
    Alkalinity 49.3
    Ca        9.1
    Cl        1.58
    Fe(3)     0.015
    K         1.26
    Li        0.0436835
    Mg        2.78
    Na        10.6
    Si        61.3786 as SiO2
    -water    1 # kg
INVERSE_MODELING 1
    -solutions      1        2
    -uncertainty    0.15
    -phases
        Kaolinite
        Montmor-Mg        pre
        Ferrosilite       dis
        Enstatite         dis
        Tremolite         dis
        Albite            dis
        Anorthite         dis
        #Spodumene      dissolve
     Goethite
    -tolerance         1e-10
    -mineral_water     true
    -minimal
END

[Brendan]

Dr. Parkhurst,

I want to thank you for your helpful and thoughtful response, and I certainly see what you are saying about the major vs. trace elements. Most of the trace elements don't really concern me, but I was trying to carry them through in case they're helpful later. Part of what I'm trying to do here is see where the Li in the perched intermediate zones comes from, for example.
I am hopeful, even though the aqueous concentrations are small, that iron and aluminum can be useful indicators in differentiating water-rock interactions by rock formation (I also have XRF data for the whole rock compositions). So while I am not, as you said, trying to use inverse modeling as a quantitative method, I am trying to use it to see the dominant interactions in a few different environments. There are a few things I keep coming across, which I wanted to ask about. I notice that you added goethite as a product, which certainly works. In every solution I have, it seems like the saturation indices for magnetite and hematite are even greater than goethite. Would one of these phases control the others, or would they precipitate together? Likewise, the SI for nontronite is always greater than those for the common clays I know a little about. That holds true even if I dial the Fe down to .001 mg/L - I get SI's of 8 for nontronite, 4 for montmorillonite, 3.5 for illite and -3.5 for smectite. I really don't know how to think about this, & I'm tempted to ignore it.
In the solutions I posted before, I had already omitted the sulfates, but with them the charge balances are well within 5%.
Thank you again for your input.

[dlparkhurst]

While Fe, Al, and Si are relatively small concentrations in solution, it is necessary to include mole balance for them because they are potentially major elements in the minerals.

Goethite vs Hematite is mostly a matter of time (or temperature and pressure). If you have some XRD information, you might be able to choose. However, for mole balance modeling it really does not matter which you pick because the stoichiometry really only differs by by a water of hydration.

I would look at what are the most likely secondary products rather than which mineral has the highest saturation index. Part of SI is strictly stoichiometry, for instance dolomite SI is often twice as big as calcite, just because there are twice as many ions in the dolomite formula.

For sulfate, atmospheric, gypsum (anhydrite), or pyrite (or other iron sulfide) are usually the likely sources for sulfate. Atmospheric can only account for relatively small sulfate concentrations. Gypsum is often not likely to be present, and pyrite can only dissolve to the extent that oxidants are present. So to get large concentrations of sulfate from pyrite, you need to react it in the unsaturated zone where there is a lot of oxygen (or have a lot of nitrate available to react.)

[Brendan]

The sulfate concentrations I am working with are generally low - below 5 mg/L in the dacite, but more like 8-10 in springs. But one spring I'm looking at is closer to 18 mg/L. Is that low enough to be from atmospheric sulfur?
I have been playing around with the secondary minerals, mostly smectite because I see this come up in thin section descriptions and well completion reports (I have a lot of data, it's Los Alamos). But another thing that makes smectite attractive or likely is that, in llnl.dat, it is a mixed-cation composition, where related clays like beidellite or montmorillonite are single-cation.
I'm trying to kind of calibrate my approach in a batch model, scaling dissolution inputs to hit target concentrations in a final solution. In particular, albite is the one, since it should be the only real source of Na. Unless I have zeolites (I don't think), there shouldn't really be a sink for Na either, aside for a small amount in smectite. My other inputs are scaled to albite.
With just plagioclase and pyroxene, I can hit the targets reasonably well, so I'm getting the idea of how much I need to dissolve. But I would also like to see the pH go from 5.6 to 8.4. Only by dissolving a lot of glass can I get close to the pH target, but then I end up with too much silica in solution. Could there also be, I don't know, chalcedony precipitation?
Then, I'm saving my final solution and using EQUILIBRIUM_PHASES to see what smectite precipitation will do. Instead of going all the way to eq., I'm using target saturation indices from my final solution (the well in the dacite). When I run the model, I notice that some of the SI's are where I put them, but others will get pushed to undersaturation. Is that because a single-cation clay removes all of the Mg, for instance, until something else that needs it is undersaturated?
Maybe this is like doing inverse modeling backwards, but I think I'm doing it to better structure my input for the inverse model later.
I see how clay-type partitioning for soils is based on climate... in the subsurface I might imagine it can't really be predicted. Would you agree with that?
Again, thank you very much. 

[dlparkhurst]

You can probably get a few mg/L of sulfate from the atmosphere if you have sufficient evaporation. In Oklahoma we figured the evaporation increased concentrations by a factor of 10 or 20 based on water balance calculations. But a few mg/L can also be obtained from pyrite oxidation. If rainwater has 0.25 mmol of O2, you can make ~0.1 mmol of SO4 from pyrite, or ~10 mg/L. Oxidation in the soil zone could produce more.

You should look at the SIs for speciation calculations for you endmember solutions. They give an idea of what is reasonable to dissolve and precipitate, although I never quite trust the SIs for alumino silicates. Ususally Al is too small to measure, and the thermodynamic data is not well constrained.

I find that trying to forward model is like whack-a-mole, if one element is right, another is wrong. I think inverse modeling is more efficient because it considers all elements simultaneously. You still have the choice of what you think is dissolving and precipitating, but it should also match the pHs.

[Brendan]

Thanks for all of your help previously.
I'm seeing how useful it can be to add phases, for inverse modeling. Not playing around with pure phase end-members seems like a huge plus, and from what I'm reading (Bowser & Jones), there is a lot better chance of coming up with meaningful mass-balance results with well-characterized phases to begin with. Especially when dealing with silicate dissolution, they say. I'm seeing it as a way to put XRF data to work, and hopefully get solutions that agree, more, with saturation indices.
I am having success with implementing these, but not with getting the kind of results I'm hoping to get. I can dissolve orthopyroxene and get smectite, for instance, but I'm expecting plagioclase to participate in this somehow, and it won't.
I also get a model that just precipitates a clay and an oxide, where I would like it to...do more?
I added my dacite_glass phase to this, and have the uncertainty set high, but still don't get very satisfying results.
Maybe I'm stumbling on something obvious.
Thanks again.

TITLE   Inverse model of regional recharge in Tschicoma Fm

SOLUTION 1   #avg spring water
   units   mg/L
   temp   6.33
   pH   6.92
   pe   6.05
   redox   pe
   Alkalinity   33.06
   Al   .93
   Ca   7.43
   Cl   1.4
   F   .083
   Fe   0.31
   Mg   2.77
   Mn   .00575
   N(5)   .0745   as N
   K   1.55
   Si   32.73 as SiO2
   Na   3.92
   S(6)   6.12 as SO4   
   -water   1 # kg

SOLUTION 2   #avg_R_37_2_S3
   units   mg/L
   temp   19.73
   pH   7.81
   pe   1.4
   redox   pe
   Alkalinity   57.14
   Al   .03031
   Ba   .01149
   B   .0118
   Br   .03
   Ca   9.86
   Cl   1.88
   F   .24
   Fe   .07956
   Li   .02297
   Mg   2.97
   Mn   .01238
   Ni   .00114
   N(5)   1.8   as NO3
   K   1.35
   Si   64.26 as SiO2
   Na   11.29
   S(6)   2.01 as SO4
   -water   1 # kg
INVERSE_MODELING 1
    -solutions   1   2
    -uncertainty   0.75
    -phases
      Kaolinite         pre
      Smectite-low-Fe-Mg   pre
        Orthopyroxene      dis
      Plag_An50         dis
   Goethite         pre
   Dacite_glass      dis
    -tolerance   1e-10
    -mineral_water   true
    -minimal
PHASES
   
Plag_An50
   Na0.5Ca0.5Al1.5Si2.5O8 + 6.0 H+ = 0.5 Na+ + 0.5 Ca+2 + 1.5 Al+3 + 2.5 SiO2 + 3.0 H2O
      log_k   0.0   #no log_k, Inverse modeling only

Orthopyroxene
   Mg0.913Fe1.087Si2O6 + 4 H+  =  + 2 H2O   + 0.913 Mg+2 + 1.087 Fe+2 + 2 SiO2
      log_k   0.0   #no log_k, Inverse modeling only
Dacite_glass
   SiAl0.26O2(OH)0.78 + 0.78H+ + 1.22H2O = SiO2 + 2H2O + 0.26Al+++
         log_k      0.10   
      delta_h   0.0
END

[dlparkhurst]

A couple of things. First, your two waters are pretty dilute, and some constituents increase and some decrease from one to the other.  So, these two waters are not that different, and I think you will have a hard time with mole balance, especially for relatively conservative constituents. Ideally, you would have several waters that demonstrate a consistent trend of increase and decreases.

I suggest for an exercise, you start with pure water and try to make your final water with reasonable reactants.

Second, you may need a source or sink for each major element. INVERSE_MODELing will determine whether there are simpler models, but in general you should probably start with a way to increase or decrease each element individually. In particular, you need sources for K, SO4, C, and perhaps others. The following gives a couple models, you can work from there.

TITLE   Inverse model of regional recharge in Tschicoma Fm
SOLUTION 3
pH 7 charge
SOLUTION 1   #avg spring water
   units   mg/L
   temp   6.33
   pH   6.92
   pe   6.05
   redox   pe
   Alkalinity   33.06
   Al   .93
   Ca   7.43
   Cl   1.4
   F   .083
   Fe   0.31
   Mg   2.77
   Mn   .00575
   N(5)   .0745   as N
   K   1.55
   Si   32.73 as SiO2
   Na   3.92
   S(6)   6.12 as SO4   
   -water   1 # kg

SOLUTION 2   #avg_R_37_2_S3
   units   mg/L
   temp   19.73
   pH   7.81
   pe   1.4
   redox   pe
   Alkalinity   57.14
   Al   .03031
   Ba   .01149
   B   .0118
   Br   .03
   Ca   9.86
   Cl   1.88
   F   .24
   Fe   .07956
   Li   .02297
   Mg   2.97
   Mn   .01238
   Ni   .00114
   N(5)   1.8   as NO3
   K   1.35
   Si   64.26 as SiO2
   Na   11.29
   S(6)   2.01 as SO4
   -water   1 # kg
INVERSE_MODELING 1
    -solutions   3   2
    -uncertainty   0.1
    -phases
      Kaolinite            pre
      Smectite-low-Fe-Mg   pre
      Orthopyroxene        dis
      Plag_An50            dis
      Goethite             pre
      Dacite_glass         dis
      chalcedony           pre
      O2(g)                dis
      CO2(g)
      K-Feldspar
      Albite
      Pyrite
    -tolerance   1e-10
    -mineral_water   true
    -minimal
PHASES
Plag_An50
   Na0.5Ca0.5Al1.5Si2.5O8 + 6.0 H+ = 0.5 Na+ + 0.5 Ca+2 + 1.5 Al+3 + 2.5 SiO2 + 3.0 H2O
      log_k   0.0   #no log_k, Inverse modeling only
Orthopyroxene
   Mg0.913Fe1.087Si2O6 + 4 H+  =  + 2 H2O   + 0.913 Mg+2 + 1.087 Fe+2 + 2 SiO2
      log_k   0.0   #no log_k, Inverse modeling only
Dacite_glass
   SiAl0.26O2(OH)0.78 + 0.78H+ + 1.22H2O = SiO2 + 2H2O + 0.26Al+++
         log_k      0.10   
      delta_h   0.0
END

[Brendan]

Dr. Parkhurst, thank you again for your help.
I see how it is necessary to have everything you're asking the model to rearrange in your phases, so here it makes sense to add K-Feldspar even if it isn't exactly present. If it isn't, and the saturation index for K-feldspar is positive (reflecting, sort of, its recalcitrance to weathering - I'm thinking of Goldich's sequence, and perhaps the idea that this melt crystallized at a lower temperature than K-feldspar would require), yet there is 2-3 wt% K in the rock, I am inclined to add it to the reaction for hydrated volcanic glass, using a molar ratio for K to Si. It seems like this fictive phase thing frees me up to do so...
In your view, is this a sensible or justifiable way to keep potassium in the picture?
Thank you again.

[dlparkhurst]

No Dr.

You should add a K phase that you think is reasonable; the K must come from somewhere. The only thing you must keep in mind is that you may need to be able to adjust each element concentration independently. As a simple example, you could not adjust both K and Mg independently if you have only one mineral that is the source for both K and Mg.  It might be possible, if the K/Mg ratio were the same in the solution and the mineral, but in general this won't be the case.

[Brendan]

Mr. Parkhurst,
Thank you again.
I think I've been successful in making a series of inverse models, based on your feedback and the minerals and conditions of my study. In general I am trying to minimize the uncertainty and still get between 2 and 6 or so solutions. I found that by specifying a greater uncertainty for sulfate with -balances, the main uncertainty can be between .05 and 0.1 usually. I'm trying to use barite as a phase, to deal with the sulfate, which seems to work sometimes but not always. I'm thinking that when the model won't really run below 0.2 uncertainty, perhaps that means mixing is important, so I am trying to test that out some.
I have set up some mixing models, using Cl as a starting point to set mixing fractions. For instance, here I am mixing groundwater with vadose zone infiltration through a wet canyon vadose zone. Before the mixing step, I'm using REACTION, and the mole phase transfers from an inverse model that I prefer, as relative moles. I did it in 4 steps, then mixed, and hit the target Cl concentration I was going for. The cations also seem reasonable, but the pH goes up much higher than I anticipated. Calcite seems to become supersaturated, so I set it to equilibrate at the end. This brought the pH down a little, but not very much. In one model it ended up around 8.8, where I was hoping to see something like 7.2; in the other it was above 9.3, where I would have liked to see about 7.5. In the final solution, where it gave C, the first model was very close when I expressed that number as HCO3, but the other not so close. Each had a few things that were quite close to the water chemistry I'm trying to recreate.
However, from looking at my silicate weathering reactions, I'm removing a lot of H+ from solution, so making it basic... but I wouldn't have expected the pH to change so much. I feel like I'm overlooking something very obvious, maybe?
Thank you again. Below is one of my inputs.

TITLE   Mixing model of intermediate recharge in Cerro Toledo to Otowi with canyon alluvial water


SOLUTION 1   #avg_R_26_S1
   units   mg/L
   temp   16.92
   pH   7.554
   Alkalinity   48.19 as HCO3
   Al   .003966
   Ba   .00861
   B   .01419
   Br   .05
   Ca   7.67
   Cl   1.3
   F   .15
   Fe   0.0239
   Mg   2.93
   Mn   .005346
   Ni   .002086
   N(5)   0.35 as NO3
   K   2.22
   Si   58.69 as SiO2
   Na   8.42
   S(6)   1.3585 as SO4
   -water   1 # kg

REACTION 1   Reacting minerals in the CT and Otowi
   Plag_An15         0.37
   rhyolite_glass_K      1.51
   OPX_CT-Otowi      1.2
   Smectite-low-Fe-Mg   -0.58
   Magnetite         -0.51
   O2(g)            0.28
   Chalcedony         -2.81
   1.0 millimoles in 4 steps
SAVE solution   3


EXCHANGE_SPECIES
        0.75Ca+2 + 0.25Mg+2 + 2X- = Ca.75Mg.25X2
        log_k   0.0
PHASES
Plag_An15
   Na0.85Ca0.15Al1.15Si2.85O8 + 4.6 H+ = 0.85 Na+ + 0.15 Ca++ + 1.15 Al+++ + 2.85 SiO2 + 2.3 H2O
Plag_An50
      Na0.5Ca0.5Al1.5Si2.5O8 + 6.0 H+ = 0.5 Na+ + 0.5 Ca+2 + 1.5 Al+3 + 2.5 SiO2 + 3.0 H2O
         log_k   0.0   #no log_k, Inverse modeling only
Plag_An52
   Na0.48Ca0.52Al1.52Si2.48O8 + 6.08 H+ = 0.48 Na+ + 0.52 Ca++ + 1.52 Al+++ + 2.48 SiO2 + 3.04 H2O
Orthopyroxene
      Mg0.913Fe1.087Si2O6 + 4 H+  =  + 2 H2O   + 0.913 Mg+2 + 1.087 Fe+2 + 2 SiO2
         log_k   0.0   #no log_k, Inverse modeling only
OPX_CT-Otowi   
   Mg0.5Fe1.5Si2O6 + 4 H+ = 2 H2O +  0.5 Mg++ +  1.5 Fe++ +  2 SiO2
Hornblende
   Ca2Mg2.28Fe2.72Si8O22(OH)2 + 14.0 H+  =  2.0 Ca++ + 2.28 Mg++ + 2.72 Fe++ + 8.0 H2O + 8.0 SiO2
Dacite_glass
      SiAl0.26O2(OH)0.78 + 0.78H+ + 1.22H2O = SiO2 + 2H2O + 0.26Al+++
      log_k      0.10   
         delta_h   0.0
dacite_glass_K
      K0.06SiAl0.27O2(OH)0.87 + 0.87 H+ = 0.87 H2O + SiO2 + 0.27 Al+++ + 0.06 K+
      log_k      0.1
      delta_h   0.0
rhyolite_glass_K
   K0.08SiAl0.2O2(OH)0.68 + 0.68 H+ =  0.2 Al+++ + SiO2 + 0.68 H2O + 0.08 K+
END

SOLUTION 30   #TDS is 216 avg   #aggregate alluvial water
units   mg/L
pH   7.0
Alkalinity   92.0 as HCO3
Al   3.455
B   .0328
Br   0.1519
Ca   19.86
Cl   23.29
F   0.24
Fe   2.165
Mg   5.8
Mn   .181
N(5)   .41
K   4.1
Si   39.26 as SiO2
Na   29.80
S(6)   14.99 as SO4
P   0.42

MIX 1
   3   0.753
   30   0.247
SAVE   solution   7
END
EQUILIBRIUM_PHASES 1
   Calcite   0.0
USE solution 7
END

 

[dlparkhurst]

So did you find a successful inverse model including mixing? If so, it gives you the mixing and mole transfers that you need for forward modeling.

It will only make sense to include equilibrium with minerals if the saturation state for the minerals are near equilibrium in a speciation calculation of the final water composition.