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
Please note that some U.S. Geological Survey (USGS) information accessed through this page may be preliminary in nature and presented prior to final review and approval by the Director of the USGS. This information is provided with the understanding that it is not guaranteed to be correct or complete and conclusions drawn from such information are the sole responsibility of the user.
Any use of trade, product, or firm names in this publication is for descriptive purposes only and does not imply endorsement by the U.S. Government.
The URL of this page is:
Last modified: $Date: 2005-09-13 21:04:21 -0600 (Tue, 13 Sep 2005) $
Visitor number 2743 since Jan 22, 1998.