Calibrate salinity model to data — calibrate_salinity

Calibrate some or all salinity model parameters a, b, d, and C_d using salinity data.

calibrate_salinity_model(
  hydro_data,
  v,
  method = "Nelder-Mead",
  control = list(trace = T)
)

Arguments

hydro_data	data.frame containing `S_ppm`, `Q_cumec`, and `year` columns
v	logged parameters as a vector: log(a), log(b), log(d), log(C_d). Use NA for calibration params
method	"Nelder-Mead" or method supplied to `stats::optim` details below
control	control parameters supplied to `stats::optim` for the `SANN` optimization

Value

The function returns a vector of all (logged) parameters including calibrated ones

Details

This function calibrates the parameters in v set to NA using streamflow and salinity in hydro_data. The calibration is done by minimizing the sum of squared error of the modeled salinity with observed salinity. The minimization is performed using stats::optim function with simulated annealing.

Note that any values of hydro_data$S_ppm that are set to NA will be ignored in the calibration. In other words, it's possible to calibrate on a subset of dates within this data.frame.

The method input can be set to "auto", in which case "Nelder-Mead" is used for multi-variate optimization and for univariate optimization ("Brent" has given poor results). This can also be set to any of the method options allowed by stats::optim -- see ?optim for details.

Examples

library(deltasalinity)
hydro_data <- ganges_streamflow
hydro_data$S_ppm_orig <- sim_salin_annual(ganges_streamflow, ganges_params$param)

# Add random error to salinity output
set.seed(100)
hydro_data$S_ppm <- hydro_data$S_ppm_orig * runif(nrow(hydro_data), min = 0.9, max = 1.1)
hydro_data$year <- as.numeric(strftime(hydro_data$date,"%Y"))

# Calibrate parameter "a"
if (FALSE) {
# Nelder-Mead produces a warning with only 1 optimization parameter. Still seems to work though.
v <- ganges_params$param
v[1] <- NA
v_calibrated <- calibrate_salinity_model(hydro_data, v) # warning expected for Nelder-Mead

# Check the percent difference
(v_calibrated - ganges_params$param) / ganges_params$param * 100
}

# Calibrate parameters "b" and "d"
v <- ganges_params$param
v[c(2,3)] <- NA
v_calibrated <- calibrate_salinity_model(hydro_data, v)
#> Calibrating using Nelder-Mead method...
#>   Nelder-Mead direct search function minimizer
#> function value for initial parameters = 1825046197.139056
#>   Scaled convergence tolerance is 27.1953
#> Stepsize computed as 0.921034
#> BUILD              3 2226637971.075970 1812323319.196938
#> HI-REDUCTION       5 1825046197.139056 1412878712.686868
#> HI-REDUCTION       7 1812323319.196938 1217538363.578517
#> HI-REDUCTION       9 1412878712.686868 936072466.602010
#> LO-REDUCTION      11 1217538363.578517 936072466.602010
#> HI-REDUCTION      13 1014693954.791372 936072466.602010
#> HI-REDUCTION      15 960287175.507668 936072466.602010
#> HI-REDUCTION      17 956612833.986446 936072466.602010
#> EXTENSION         19 947792789.570164 899945211.639251
#> EXTENSION         21 936072466.602010 852289426.342894
#> EXTENSION         23 899945211.639251 729670267.997177
#> REFLECTION        25 852289426.342894 655988309.168991
#> HI-REDUCTION      27 753493539.733399 655988309.168991
#> HI-REDUCTION      29 729670267.997177 655988309.168991
#> EXTENSION         31 710859350.940728 575860139.221862
#> LO-REDUCTION      33 655988309.168991 575860139.221862
#> EXTENSION         35 590301634.712233 445367043.592501
#> LO-REDUCTION      37 575860139.221862 445367043.592501
#> REFLECTION        39 455879139.932054 334267334.274454
#> HI-REDUCTION      41 445367043.592501 334267334.274454
#> EXTENSION         43 410436430.812456 249978434.948027
#> REFLECTION        45 334267334.274454 215218656.331216
#> EXTENSION         47 249978434.948027 91257289.895338
#> HI-REDUCTION      49 215218656.331216 91257289.895338
#> EXTENSION         51 173424769.177592 51373184.363852
#> LO-REDUCTION      53 91257289.895338 36960705.892848
#> LO-REDUCTION      55 51373184.363852 25819167.193135
#> HI-REDUCTION      57 36960705.892848 25819167.193135
#> HI-REDUCTION      59 28512521.819791 21703028.228678
#> REFLECTION        61 25819167.193135 18427265.737055
#> HI-REDUCTION      63 21703028.228678 18427265.737055
#> LO-REDUCTION      65 19791038.087178 18427265.737055
#> HI-REDUCTION      67 18763606.802179 18427265.737055
#> HI-REDUCTION      69 18575889.764778 18295844.973268
#> HI-REDUCTION      71 18427265.737055 18226056.201577
#> HI-REDUCTION      73 18295844.973268 18156044.431435
#> LO-REDUCTION      75 18226056.201577 18156044.431435
#> HI-REDUCTION      77 18167181.112142 18156044.431435
#> REFLECTION        79 18160344.750118 18149369.487342
#> HI-REDUCTION      81 18156044.431435 18146861.589956
#> HI-REDUCTION      83 18149369.487342 18145887.134497
#> HI-REDUCTION      85 18146861.589956 18145661.938706
#> HI-REDUCTION      87 18145887.134497 18145160.297985
#> HI-REDUCTION      89 18145661.938706 18144871.169847
#> HI-REDUCTION      91 18145160.297985 18144871.169847
#> HI-REDUCTION      93 18145049.820170 18144871.169847
#> LO-REDUCTION      95 18144921.732866 18144871.169847
#> HI-REDUCTION      97 18144875.322017 18144835.420098
#> REFLECTION        99 18144871.169847 18144834.990372
#> Exiting from Nelder Mead minimizer
#>     101 function evaluations used
v_calibrated
#> [1] -12.560625  -4.176989  -6.174101  10.463103

# Check the percent difference
(v_calibrated - ganges_params$param) / ganges_params$param * 100
#> [1]  0.0000000 -0.7239763 -0.2124110  0.0000000