Re: PHREEQC problem

[xiaomin Mao <xmao@xxxxxxxx>]

David,

Thank you very much! It can work very well now.
Yes, the .0001(org_cuex1 + org_cuex2) = kin(org_xcu)" represents equilibrium. In fact it's only part of the Nica-Donnan model, as at first
I think it is too difficult to model it using phreeqc. Now maybe I should continue to explore how to perform the whole model.
I am further convinced that PHREEQC can do anything you want.  :-)

cheers,
xiaomin

David L Parkhurst wrote:

> > 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