Convergence problems in coupled reactive transport using user-defined species.

[johan_s12345]

Dear phreeqcusers,

I run reactive transport simulation of a soil filter used to treat wastewater using Iphreeqc coupled to a transport simulator.
Iphreeqc solves a couple of kinetic reactions (RATES, KINETICS) of user-defined SPECIES such as:

example:

SOLUTION_MASTER_SPECIES
    Do   Do   0   1   1
    Sf           Sf     0   1   1

SOLUTION_SPECIES
Do = Do
    log_k     0
Sf = Sf
    log_k     0

The solution bloc when calling IPhreeqc looks like:

SOLUTION 1
units   mol/kgw
temp  20.0                  # temp
density   1
pH       7.0000     # comp 1
pe      4.0            # comp 2
Ogas          1.0            # comp 3
Do       9.0e-3            # comp 4
Sf        0.5e-3             # comp 5
Sa      1.03-3   # comp 6
Sin      0.5e-3   # comp 7
Snh      0.5e-3   # comp 8
Sno      40.0e-3   # comp 9
Sulphide   10.00e-3   # comp 10
Sso      10.0e-3   # comp 11
Xs_m      0.5e-3   # comp 12
Xs_im      0.5e-3   # comp 13
Xi_m      0.1e-3   # comp 14
Xi_im      0.1e-3   # comp 15
Xh        1.36e-4         # comp 16
Xa        1.36e-4         # comp 17
Xfb        1.36e-4        # comp 18
Xamb        1.36e-4     # comp 19
Xasrb      1.36e-4  # comp 20
Xsob      1.36e-4    # comp 21
-water       1    # kg      
#ende

However, IPhreeqc chrashes when the concentrations of a species get higher than 55mol/kgw with the error messages:

---
ERROR:                 H(1) Charge balance has not converged.    Residual: -9.208858e-16

ERROR:               A(H2O) Activity of water has not converged.    Residual: 6.845678e-05

ERROR: Model failed to converge for initial solution.
---

I think the molalities just get too high to allow convergence.

I tried to add a dummy species as described in the manual:

SOLUTION_SPECIES
H2O + 0.01e- = H2O-0.01     
log_k   -9.0

but without success.

Do you have any suggestions of how to improve/allow the convergence for charge and activity?

In fact, I do not need any equilibrium reactions computed by phreeqc, I just need the capability to solve user-defined reactions in RATES and KINETICS.

Best regards,

[dlparkhurst]

It does not look like your concentrations are too high.  It looks like you are defining immobile and sorbed species in the SOLUTION; normally EXCHANGE or SURFACE are used for the solid phases and only dissolved concentrations are defined in the SOLUTION

I would need your database file and an input file that fails to tell what the problem is.

[johan_s12345]

Dear Mr. Parkhurst,

thanks for your reply.

You are right, I specified dissolved and immobile species in the SOLUTION bloc.
I think that EXCHANGE and PHASES do not really serve my needs as the immobile species are not minerals that dissolve/precipitate accoring to a chemical equilibrium.

My model consists of different types of bacteria (Xh, Xa, Xfb, Xamb, Xasrb, Xsob) that grow by degrading different substrate material (Sf, Sa, Snh, Sno, Sulphie, Sso, Xs). The bacterial growth us modeled using kinetic reactions, where each reaction is a KINETIC species.
The bacteria are considered as immobile in the transport simulator.
Additionally, there are two components (Xs, Xi) that exist in dissolved (Xs_m, Xi_m) and immobilized form (Xi_im, Xs_im). These are solid components that can get transported in the water phase and attach to the porous media.
When running long simulations the attached components (Xi_im) gets pretty high (>50g/L) and the simulation crashes with the error message I posted before.

The input file for the transport simulator is: inputfile_Iphreeqc.pqc
The corresponding databasefile: phreeqc.dat  (user-defined stuff in SOLUTION_SPECIES, SOLUTION_MASTER_SPECIES, CALCULATE_VALUES).

A replication in phreeqc only (batch reaction) is: mod_phreeqc.pqi   using same database file.

The model mod_phreeqc.dat crashes immediately when defining a high concentration (Xi_im of 50000). I think this happens in the original coupled transport simulatin as well. Iphreeqc crashes when called if the concentration of Xi_im higher than a certain value.

Best regards,

[dlparkhurst]

I don't think PHREEQC can handle the high concentrations. The activity of water is calculates as 1 - sum(solutes), so at some point the activity of water will be negative. This may be the limit on your concentration, but it also will affect pH and potentially all of the aqueous equilibria. In general, the ion association model will have trouble with concentrations of 10 molal or more.

Perhaps you could scale everything by a factor of 1e-3 so the concentrations are not so large.

Also, 1 molal O2 is very large.

[johan_s12345]

I already scaled to g/L as, beside the immobilized species, mg/L is  more appropriate unit. I am a little afraid with respect to computation accuracy when scaling another 1e-3, however, I may try and see what the error margin looks like.

Would it be possible to pass by the activity computation in SOLUTION by defining all my transport simulator components as KINETIC reactions instead of SOLUTION species?

[dlparkhurst]

You cannot bypass the speciation (activity) calculation, but if you use KINETICS, you will not have the unreasonable concentrations of dissolved species. With some work, you could use your transport module to transport the moles of a KINETIC reactant, in addition to the dissolved constituents, when the reactant is assumed to be mobile.