# 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

  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

  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

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

#  -initial_time 0

#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 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


  -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")