# 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.
# About half of the original dataset was deleted to give an irregular (in time) series.
# Only one value is sent to the KINETIC block at a time which means that you do not need a regular time series of observations
# i.e. uses PhreePlot looping rather than PHREEQC's own internal (and faster) looping
SPECIATION
calculationType fit
calculationMethod 1
fitMethod nlls
dataFile kineticsSiobsirregular.dat # contains a list of observations (with errors) at irregular times
onePass FALSE # One PHREEQC run per observation rather than one run for the whole set of observations.
# Slower but enables irregular time series to be fitted without much effort.
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.0E3 # initial step size for each adjustable parameter (here log_k)  MUST give a 'significant' change in the dependent variable
blockRangeColumn sim # column header in fit data file for which PHREEQC simulation(s) to use for each observation  here only 1 simulation anyway
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
numericTags <secs> = <time>*3.1536e7 # <time> form fit data file, 3.1536e7 secs in one year
PLOT
plotTitle "Estimate log_k for kinetics of quartz dissolution<br>(synthetic data, irregular 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 kineticsSiobsirregular.dat # use this file to get the initial data for plotting
extradat kineticsSifit1.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 'Getgoing sheet #11' (www.xs4all.nl/~appt)
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 <secs> # substitute here, reading down fit data file
tol 1e8 # integration tolerance, default 1e8 mol
# INCREMENTAL_REACTIONS FALSE # not necessary when only one 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
