Re: PHREEQC problem

["David L Parkhurst" <dlpark@xxxxxxxx>]

> Sorry I didn't explain the process clearly. I assume the sorption of Cu
by
organicmatter site (Org_x) as follow (see Table 6.6, Werner Stumm, 1995,
the
book: Aquatic Chemistry):
Cu+2  + Org_xH2  =  Org_xCu  + 2H+


> To the formula of Org_xh,  I mean H is 0 (zero), not H O. Maybe as you
said,
it's not necessary to mention this kinetic expression as it has the same
rate
like Org_xcu.

Sorry, my bad. I am a little slow, but I understand your input file now. I
don't have the latest Stumm and Morgan, but I did some investigation into
your problem. I am assuming the ".0001(org_cuex1 + org_cuex2) =
kin(org_xcu)" represents equilibrium; if not, then I what follows would
have to be modified.

Here are a few comments:


(1) First you are using kinetics to represent equilibrium and PHREEQC is
not particularly good for this kind problem. PHREEQC is not good for
handling stiff differential equations, an implicit integration method would
work better than the explicit Runge-Kutta routine that is currently used.
However, I think PHREEQC can beat the problem to death, although slowly.
Better yet would be to include the sorption model as a type of equilibrium
for PHREEQC, which would probably be hundreds of times faster.

(2) I'm not sure of the database you are using, but if it has both Cu(1)
and Cu(2), you could be running into some redox problems. You should define
the Cu in solution as "Cu(2)", otherwise, you may be removing Cu(1) and
effectively replacing it with one H+ and one H(0). The H(0) will cause some
unwanted redox reactions. With wateq4f.dat, if you define copper as Cu in
your initial solution, a significant fraction is Cu(1), so be careful that
the Cu starts out in the right redox state [Cu(2)].

(3) Aside from the redox-state issue, I think the fundamental problem is
integrating the rate expression. You have a dissolved Cutot [TOT("Cu")] and
sorbed Cukin [KIN(org_cux]. You calculate the sorbed Cu from the Cutot (and
H+), compare it to the Cukin, and calculate a rate from the difference
divided by 10. The rate is sufficient to ensure equilibrium in about 10
seconds; however, phreeqc attempts to extrapolate the rate for times up to
the time step (30000 seconds). There tends to be overshoot and undershoot
if the rate is extrapolated for two reason: (1) the long extrapolation
relative to the time to equilibrium, and (2) poor estimate of equilibrium.
The time scale of 10 seconds for the rate probably limits the internal time
step for the integration to about 1 second. thus there are a lot of
internal iterations to integrate the rate expression in each cell.

As for the poor estimation of equilibrium, the original rate expression
does not account for the change in solution and sorbed concentrations to
reach equilibrium. Let equilibrium be represented by Cutot,eq = Cutot +
dM/TOT("water") and Cukin,eq = Cukin - dM for some dM, where dM is moles.
You have equilibrium as dM = f(CuTot, H+) - Cukin, where CuTot and CuKin
are fixed at the values at the beginning of the integration step. I think
this leads to the possibility of oscillation, where one integration step
overshoots and the next over-corrects.

I've written a little routine that estimates dM and then calculates the
rate as a fraction of dM. I've used dM/10000 as discussed below. Regardless
of the factor used, the integration appears to be robust, only more or less
slow.

(4) Your dispersivity is large relative to the cell size and column length.
This (and -correct_disp) requires 9 "mixruns" to simulate a time step. As
you know, a smaller dispersivity would run faster. With as much dispersion
as you have, it would be considerably faster and equally accurate (I think)
to run the problem with phast; it would not require the 9 "mixruns" to
account for the dispersion.

(5) Finally, a rate that is of the same scale as the time step will run
faster (less stiff). I've used a rate that should be close to equilibrium
in 10000 seconds. I experimented a little with faster rates and the
accuracy seemed sufficient with the slower (dM/10000) rate. With this rate,
I was able to run the entire simulation in about 6 minutes on an early
pentium IV.

David



DATABASE ../database/wateq4f.dat
Rates
Org_xcu
-start
10 Cu_conc0=Tot("Cu")
20 H_conc=ACT("H+")
30 sorbed0 = M
35 moles_cu = cu_conc0 * TOT("water") + sorbed0
#40 print "cell, cu, sorbed, moles_cu: ", cell_no, cu_conc0, sorbed0, moles_cu, time, total_time
50 if (moles_cu <= 1e-14) then goto 350
60 if (Cu_conc0 > 1e-14) then x = log10(Cu_conc0/2) else x = log10(.01*moles_cu/TOT("water"))
#80 print "10^x, 10^x_cu, 10^x_tot", 10^x, Cu_conc0/2, .01*moles_cu/TOT("water")
130 REM loop to calculate equilibrium
140 gosub 500
150 fx1 = fx
160 REM calculate df(x)/dx
165 dx = 1e-7
170 x = x + dx
180 gosub 500
185 fx2 = fx
190 dfx = (fx - fx1)/dx
#200 x = x - dx - fx1/dfx
210 if (ABS(fx1/dfx) < 1) then x = x - dx - fx1/dfx else x = x - dx - (fx1/dfx)/ABS(fx1/dfx)
#220 print "cu_est, sorbed_est, resid, delta, mass_err", 10^x, sorbed, fx1, ABS(fx1/dfx), (10^x*TOT("water") + sorbed) - moles_cu, dfx
230 if ABS(fx1) > 1e-14 then goto 130

#0.0001*Cu_conc/(Cu_conc+H_conc)*parm(1)
#35 rate1=-(0.0001*(Org_xcue1+Org_xcue2)-Kin("Org_xcu"))/10
#335 rate1 = -(sorbed-kin("Org_xcu"))/10     # diff from equilibrium / 10
335 rate1 = -(sorbed-kin("Org_xcu"))/10000     # diff from equilibrium / 10000
# 0.0001--based on the quantity of organic matter
337 put(rate1,1)
340 moles=rate1 * time
345 if ABS(moles) < 1e-14 then moles = 0
350 save moles
360 goto 1000

# calculates sorbed_Cu from cu_conc, H_conc
500 REM x is log10 of Cu concentration
510 Cu_conc = 10^x
520 sorbed = sorbed0 - TOT("water")*(Cu_conc - Cu_conc0)
530 Org_xcue1=4.722*(1.8197*Cu_conc)^0.53/((1.8197*Cu_conc)^0.53+(218.78*H_conc)^0.66) \
*((1.8197*Cu_conc)^0.53+(218.78*H_conc)^0.66)^0.59/(1+((1.8197*Cu_conc)^0.53+(218.78*H_conc)^0.66)^0.59)
540 Org_xcue2=0.8811*(181970086*Cu_conc)^0.36/((181970086*Cu_conc)^0.36+(398107171*H_conc)^0.76) \
*((181970086*Cu_conc)^0.36+(398107171*H_conc)^0.76)^0.7/(1+((181970086*Cu_conc)^0.36+(398107171*H_conc)^0.76)^0.7)
550 fx = .0001*(org_xcue1 + org_xcue2) - sorbed
560 return

1000 REM end of program
-end
Org_xh
-start
35 rate1=get(1)
40 moles=-rate1*time
50 save moles
-end
end


SOLUTION 0  CCA concentration
  units mol/kgw
  pH 7.0 charge
  pe 4
  -water 0.03375
   Na 1
   Cl 1



SOLUTION 1  CCA concentration
  units mol/kgw
  pH 4.0
  pe 4
  -water 0.03375

   Cu(2)  0.123562
   Na 1
   Cl 1.3 charge

SOLUTION 2-16  CCA concentration
  units mol/kgw
  pH 4.0 charge
  pe 4
  -water 0.03375
  Cu  0.00
Na 1
Cl 1


kinetics 1-16
 Org_xcu
 -runge_k 6
 -Formula  H -2.0 Cu +1
 -m0 0.00
 -m 0.00
 -tol 1e-10
#-Parms 0.1
Org_xh
 -Formula H 0
 -m0 0.1
 -Parms 0.1

TRANSPORT
  -cells 16
  -shifts  33
  -time_step 30857  #10.5d/28*86400s/d
  -flow_direction forward
  -lengths 16*5e-3
  -boundary_conditions flux flux
  -dispersivities 16*0.012
  -correct_disp true
  -diffusion_coefficient 0.0000
  -punch_cells 16
  -punch_frequency 1
  -warnings false

SELECTED_OUTPUT
  -file mao.7.sel
  -selected_out true
  -user_punch true
  -reset false

USER_PUNCH
-headings time    cell  Cu    pH    Org_xcu     Tot_Cu fx
  -start
10 Cu_conc=Tot("Cu(2)")
20 H_conc=ACT("H+")
30 Org_xcue1=4.722*(1.8197*Cu_conc)^0.53/((1.8197*Cu_conc)^0.53+(218.78*H_conc)^0.66) \
*((1.8197*Cu_conc)^0.53+(218.78*H_conc)^0.66)^0.59/(1+((1.8197*Cu_conc)^0.53+(218.78*H_conc)^0.66)^0.59)
40 Org_xcue2=0.8811*(181970086*Cu_conc)^0.36/((181970086*Cu_conc)^0.36+(398107171*H_conc)^0.76) \
*((181970086*Cu_conc)^0.36+(398107171*H_conc)^0.76)^0.7/(1+((181970086*Cu_conc)^0.36+(398107171*H_conc)^0.76)^0.7)
50 fx = .0001*(org_xcue1 + org_xcue2) - Kin("Org_xcu")
100  PUNCH total_time, cell_no, ToT("Cu"), -LM("H+"),Kin("Org_xcu"), Cu_conc*TOT("water") + KIN("Org_xcu"), fx
 -end

USER_Graph
  -start
10  PUNCH ToT("Cu"), -LM("H+"),Kin("Org_xcu"),Kin("Org_xh")
 -end

END

David Parkhurst (dlpark@xxxxxxxx)
U.S. Geological Survey
Box 25046, MS 413
Denver Federal Center
Denver, CO 80225

Project web page: https://wwwbrr.cr.usgs.gov/projects/GWC_coupled