Error when running reactive transport modelling with PHREEQCRM

[Yongqiang]

Hi David,

I set up a simulation by coupling phreeqcrm and openfoam. The simulations can run well with a bigger time step. However, with a shorter time step, a few errors popped out after running a while. The following is the error message. Do you have any suggestion to fix this?

Regards,
Michael

Code: [Select]

WARNING: Maximum iterations exceeded, 100

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying smaller step size, pe step size 10, 5 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying reduced tolerance 1e-16 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying increased tolerance 1e-14 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying diagonal scaling ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying diagonal scaling and reduced tolerance 1e-16 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying scaling pure_phase columns 1e-10 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying scaling pure_phase columns and diagonal scale 1e-10 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying increased scaling 1e-09 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Skipping optimize equations for first 5 iterations ...

WARNING: Maximum iterations exceeded, 100

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Adding inequality to make concentrations greater than zero.

WARNING: Maximum iterations exceeded, 100

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying reduced tolerance 1e-17 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying reduced tolerance 1e-18 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: The program has failed to converge to a numerical solution.

The following equations were not satisfied:
ERROR:               A(H2O) Activity of water has not converged. Residual: 3.878177e-04

ERROR: Numerical method failed on all combinations of convergence parameters, cell/soln/mix 0
Stopping.
ERROR: ERROR:               A(H2O) Activity of water has not converged. Residual: 3.878177e-04

ERROR: Numerical method failed on all combinations of convergence parameters, cell/soln/mix 0

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

ERROR: Numerical method failed on all combinations of convergence parameters, cell/soln/mix 0

ERROR: PhreeqcRM failed.
ERROR: PhreeqcRM::RunCells


[dlparkhurst]

I'm not sure if I can provide a solution or not. Changing some of the convergence parameters might help.

What time step are you using and did it happen at the first time step?

It looks like the first cell of your numbering system (cell 0) failed. It would help to see the input file that defines the initial solutions and reactants. What is the solution volume of the first cell (rv, porosity, saturation)? Is it a cell with a boundary condition?

Sometimes it helps to add a species to poise redox in systems with very little redox buffering. (Makes sure this code is run by the workers):

Code: [Select]

SOLUTION_SPECIES
H2O + 0.01e- = H2O-0.01
log_k -8
END

[Yongqiang]

Many thanks for the comments, David.

The following script is the aqueous reactions, set in the initialization part:

Code: [Select]

SELECTED_OUTPUT
-reset false
-high_precision true
-molalities AT H2O H+
END

SOLUTION_MASTER_SPECIES
T T 0 1 1
A A 0 1 1

SOLUTION_SPECIES
T = T
A = A
A + T = AT
log_k 2

SOLUTION 0 formation water
     units   Mol/L
     T 0.00001
     A 0.00001
END

SELECTED_OUTPUT
-reset false
-high_precision true
-molalities AT H2O H+
END

The species array is as the following:

Code: [Select]

0  A
1  AT
2  H+
3  H2
4  H2O
5  O2
6  OH-
7  T

The species concentration after one step of transport is as the following:

Code: [Select]

DILUPBiCGStab:  Solving for T, Initial residual = 1, Final residual = 1.34983e-22, No Iterations 1
0  1.00201
1  1.108e-08
2  1.11063e-07
3  8.80961e-26
4  61.6451
5  0
6  1.07006e-07
7  1.00201
          Estimated efficiency of chemistry without communication: 58.3619
          Cells shifted between threads     0
          Time rebalancing load             1.16015e-05
0  0.0901333
1  0.812401
2  9.92048e-08
3  1.44626e-39
4  55.5253
5  2.94271e-15
6  9.5551e-08
7  0.0901333
The top species is the concentration after one step transport of the initialization species. The bottom concentration is the reacted species concentration for the transported species.

The boundary condition is set in the openfoam with a fixed value type. Do you have any suggestion to set it directly in Phreeqcrm? I tried to set the boundary condition directly in Phreeqcrm. However, it is not succeeded.

Thanks

[dlparkhurst]

I think you may have issues with units. Your SOLUTION 1 has A and T concentrations of 1e-5 mol/L. If you are using mg/L as your transport units, that would be about 0.010 mg/L (using gram formula weights of 1 for A and T). So, I am confused where your concentration of 1 for A and T in your top concentrations.

The 1 for A and T in the top concentrations appears to be mol/L. Look at the concentration of species AT. At 1e-5 mol/L of A and T, the concentration of AT is smaller than A and T, somewhat like the upper set of concentrations. At 1 mol/L of A and T, AT is the dominant species. In your lower concententration set, AT is dominant.

Something else is wrong with H2O. I don't see how it's concentration can change from 61 to 55 (while total A and T concentrations are the same in top and bottom sets) if the two sets of concentrations have the same units.

Code: [Select]

SOLUTION_MASTER_SPECIES
T T 0 1 1
A A 0 1 1

SOLUTION_SPECIES
T = T
A = A
A + T = AT
log_k 2

SOLUTION 0 formation water
     units   Mol/L
     T 1e-5
     A 1e-5
SOLUTION 1 formation water
     units   Mol/L
     T 1e0
     A 1e0
END

So, I am suspicious that you are not transferring the concentrations correctly. You should pick the transporting concentration unit with SetUnitsSolution--mol/L, mg/L, or mass fraction. You should transport in one of these units. The SetSpeciesConcentrations will convert your transporting units to PHREEQC solution definitions, so you should not have to do any conversions on your side.

If you haven't already, I would start with simple Na and Cl to make sure you are doing conservative transport correctly. Make sure you first transfer your initial conditions with InitialPhreeqc2Module. Then use GetSpeciesConcentrations to get the concentrations for transport, then SetSpeciesConcentrations after transport. (Unless you are doing multicomponent diffusion calculations, you can transport total concentrations [SetConcentrations and GetConcentrations] rather than species concentrations.)

PhreeqcRM does not care about boundary conditions, that is totally on the transport side. The strategy is to use InitialPhreeqc2Concentrations to get the concentrations you need to apply boundary conditions in the transport calculation.

[Yongqiang]

Hi David,

The code has unsteady velocity problem. The inlet velocity was 10 while the internal velocity was set to be 5. So this could be one of the reasons for the strange output.

When I run the following phreeqc code, the inlet concentration is 10 mol/l for Na+ and Cl-. However, the selected output for Na+ is 24 mol/l, which is higher than the inlet concentration. Could you please help find the reason for this high concentraton?

Code: [Select]


SOLUTION 0 formation water
     units   mol/l
     Na 10.00001
     Cl 10.00001
END

SOLUTION 1-100 formation water
     units   mol/l
     Na 10.00001
     Cl 10.00001
END

USER_PUNCH
    -head Na+ PH K_mmol Cl_mmol
    10 PUNCH MOL("Na+"), -LA("H+") #TOT("K")*1000, TOT("Cl")*1000
SELECTED_OUTPUT
-reset false
-high_precision true
#-molalities Na+ Cl- H2O
END

TRANSPORT           
   -cells   100
   -shifts  100
   -flow_direction  forward
   -boundary_conditions flux  flux
   -lengths 0.1
   -dispersivities        0.50     # No dispersion
   -diffusion_coefficient 0.3e-9     # No diffusion
   -time_step              1  # 317 years give 122 mixes
   -punch_cells 1
   -punch_frequency 1
USER_GRAPH 1 Example 11
  -chart_title "Transport data"
  -headings NaCl Na Cl
  -axis_titles "PORE VOLUMES" "MOLES PER KILOGRAM WATER"
  -axis_scale x_axis auto
  -axis_scale y_axis auto
  -plot_concentration_vs time
  -start
  10 x = (STEP_NO + 0.5) / cell_no#DIST
  #20 PLOT_XY x, TOT("Li")*1000, symbol = None
  #30 PLOT_XY x, -LA("H+"), symbol = None
  #30 PLOT_XY x, MOL("NaCl"), symbol = None
  40 PLOT_XY x, MOL("Na+"), symbol = None
  #50 PLOT_XY x, MOL("Cl-"), symbol = None
  -end
END

[dlparkhurst]

I'm not convinced that you are transporting correctly, so first, you should run conservative transport time steps without any calls to PhreeqcRM at all. For testing, you should use upstream-in-space and backward-in-time weighting, which should guarantee all concentrations calculated will be between the maximum and minimum concentrations; other weighting could result in oscillations.

Once you are satisfied that transport is working correctly, you should skip RunCells, simply SetConcentrations and GetConcentrations, which should produce the same concentrations before and after the calls.

Finally, you can add RunCells, and, assuming you do not define any reactions (EQUILIBRIUM_PHASES, EXCHANGE, and others), should also produce the same concentrations before and after.

At that point, you could start trying other weighting strategies, reactions, boundary conditions, etc.

[dlparkhurst]

During testing, I would also use concentrations of 1 mmol/kgw or less in your Phreeqc definitions to avoid major changes in density effects.

Note that when you use mol/L and are using phreeqc.dat or pitzer.dat, you should probably use -density 1 calc in the SOLUTION definition. Otherwise the mol/L calculated for transport will not match the input mol/L.

[Yongqiang]

Thanks, David!
This is a really good, really useful, smart programming method!
I followed the method to check the code step by step. Before Run_cells, the output and mass balance were all fine. However, after activating Run_cells, all the concentrations were normal except the AT (A+ T = AT). The inlet and outlet boundary values for A and T were set as 10 mol/L. During the running, the AT concentrations went up to 20 mol/L. I tried to fix it. However, the concentration still went up after running an enough long time. Could you please share some insights in this problem?

Many thanks.

[dlparkhurst]

The total concentrations of A and T should not have changed by RunCells (assuming only SOLUTIONS have been defined), except possibly for some density effects. The changes due to density should be small, not double. The concentration of an aqueous species could change by a factor of 2 in RunCells due to an aqueous equilibrium calculation (but total A and total T should not change much).

You don't say, but for getting started, you should use GetConcentrations and SetConcentrations. You should transport all of the components included in these two methods, which are known after a call to FindComponents.

If your concentrations are gradually changing after repeated calls to RunCells, it could be because of changes in density. Here is an excerpt from the documentation:

Code: [Select]

Two options are available for the volume and mass of solution that are used in converting to transport concentrations: (1) the volume and mass of solution are calculated by PHREEQC, or (2) the volume of solution is the product of saturation (SetSaturation), porosity (SetPorosity), and representative volume (SetRepresentativeVolume), and the mass of solution is volume times density as defined by SetDensity. UseSolutionDensityVolume determines which option is used. For option 1, the databases that have partial molar volume definitions needed to accurately calculate solution volume are phreeqc.dat, Amm.dat, and pitzer.dat.
The units conversion from mol/L (provided in SetConcentrations) to moles in a solution for PHREEQC calculations is done as follows:

Code: [Select]

To convert from mg/L to moles of element in the representative volume of a reaction cell, mg/L is converted to mol/L and multiplied by the solution volume, which is the product of porosity (SetPorosity), saturation (SetSaturation), and representative volume (SetRepresentativeVolume). To convert from mol/L to moles of element in the representative volume of a reaction cell, mol/L is multiplied by the solution volume.

To convert from moles in solution in PHREEQC to mol/L (which is returned in GetConcentrations):

Code: [Select]

To convert from moles of element in the representative volume of a reaction cell to mg/L, the number of moles of an element is divided by the solution volume resulting in mol/L, and then converted to mg/L. To convert from moles of element in a cell to mol/L, the number of moles of an element is divided by the solution volume resulting in mol/L. To convert from moles of element in a cell to mass fraction, the number of moles of an element is converted to kg and divided by the total mass of the solution. Two options are available for the volume and mass of solution that are used in converting to transport concentrations: (1) the volume and mass of solution are calculated by PHREEQC, or (2) the volume of solution is the product of porosity (SetPorosity), saturation (SetSaturation), and representative volume (SetRepresentativeVolume), and the mass of solution is volume times density as defined by SetDensity. Which option is used is determined by UseSolutionDensityVolume.

If you use UseSolutionDensityVolume(false), there should not be changes due to the calculated change in density.

I will send you a message through the forum, and we can communicate further offline.