# 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

  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.0E-3                        # 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

  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

# Kinetics of quartz dissolution
# from Appelo 'Get-going sheet #11' (

#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
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 1e-8                                                            # integration tolerance, default 1e-8 mol

# INCREMENTAL_REACTIONS FALSE                                        # not necessary when only one step


  -reset FALSE

  -heading time Si_calc SIQtz                                        # these are accumulated in the out file ready for plotting
  10 IF (step_no>0) THEN punch total_time/3.1536e7, tot("Si")*1e3, SI("Quartz")