kineticsSifit.ppi
|
|
|
# Fit log_k to synthetic data for the dissolution of quartz vs time based on the \kineticsSi\kineticsSi.ppi demo
# A random error of (RAND()-0.5)*0.005 has been added to each calculated value to produce a set of observations with a small random error
SPECIATION
calculationType fit
calculationMethod 1
fitMethod nlls
dataFile kineticsSiobs.dat # contains a list of observations (with errors) and independent variables
onePass TRUE # ALL points can be calculated in one PHREEQC simulation with step option providing it is a regular time series
selectedOutputLines auto # pick up all the lines output from the first simulation
dependentVariableColumnObs Si_obs # name of column in fit data file containing the observations (dep variable) (alias dependentVariableColumnObs)
dependentVariableColumnCalc Si_calc # name of column in selected output file containing the calcd values of the dep variable (alias dependentVariableColumnCalc)
fitWeightingMethod 0 # 0 = unit weights
fitFiniteDiffStepSize 1.0E-3 # initial step size for each adjustable parameter (here log_k) - MUST give a 'significant' change in the dependent variable
numberOfFitParameters 1
fitParameterNames log_k
fitLogParameters 0 # 0 = parameters on a linear scale, 1 = log10 parameters before fitting - log_k is ALREADY on a log scale so do not log again
fitAdjustableParameters 1 # 1 = adjustable
fitParameterValues -13 # initial estimate - needs to be quite close
fitMaxStepSize 2
# numericTags <log_k> = -13.7
PLOT
plotTitle "Estimate log_k for kinetics of quartz dissolution<br>(synthetic data, regular series)"
xtitle "Time (year)"
ytitle "Si (mmol/kgw)"
customXColumn time # plot x = time
lines Si_calc # plot y(line) = Si
pointcolor blue
points Si_obs residuals
extradat kineticsSiobs.dat # use this file to get the initial data for plotting
extradat kineticsSifit.pts
linecolor blue
linewidth 0.6
labelSize 0 # suppresses label
legendTextSize 0 # suppresses key
extraText extratextkineticsfit.dat # additional text (including some Greek symbols and goodness of fit)
extrasymbolslines extrasymbolslinesfit.dat # add dashed line at y=0
uselinecolordictionary 0
xaxislength 150
CHEMISTRY
# Kinetics of quartz dissolution
# from Appelo 'Get-going sheet #11'
#TRANSPORT
# -initial_time 0
RATES
#1 dQu/dt = -k * (1 - SR(Quartz). k = 10^-13.7 mol/m2/s (25 C)
#2 parm(1) = A (m2), parm(2) = V (dm3) recalculate to mol/dm3/s
Quartz # rate name
-start
10 moles = parm(1) / parm(2) * (m/m0)^0.67 * 10^<log_k> * (1 - SR("Quartz"))
20 save moles * time # integrate. save and time must be in rate definition
-end # moles count positive when added to solution
KINETICS # Sediment: 100% qu, grain size 0.1 mm, por 0.3, rho_qu 2.65 kg/dm3
Quartz # rate name
-formula SiO2
-m0 102.7 # initial moles of quartz
-parms 22.7 0.162 # parameters for rate eqn. Here:
# Quartz surface area (m2/kg sediment), water filled porosity (dm3/kg sediment)
-steps 1.5768e8 in 100 steps # 1.5768e8 seconds = 5 years
-tol 1e-8 # integration tolerance, default 1e-8 mol
INCREMENTAL_REACTIONS false # start integration from previous step
SOLUTION 1
SELECTED_OUTPUT
-reset FALSE
USER_PUNCH
-heading time Si_calc SIQtz # these are accumulated in the out file ready for plotting
-start
10 IF (step_no>0) THEN punch total_time/3.1536e7, tot("Si")*1e3, SI("Quartz")
-end
END
|
|