[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

*To*: xiaomin Mao <xmao@xxxxxxxx>*Subject*: Re: PHREEQC problem*From*: "David L Parkhurst" <dlpark@xxxxxxxx>*Date*: Mon, 2 Jun 2003 14:47:57 -0600*Bcc*: "David L Parkhurst" <dlpark@xxxxxxxx>*In-reply-to*: <3ED4A9C9.6D100711@ed.ac.uk>

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