Generalidades

La calibración del modelo se hace de forma conjunta en términos de cantidad y calidad de agua utilizando una función objetivo que indica el grado de similitud entre los valores simulados y observados de caudal diario (m³/s) y la concentración diaria de fósforo total (PT). Para establecer las concentraciones diarias de PT se utilizaron las relaciones exponenciales que se presentan en la sección Sistema de Monitoreo:muestreo de fósforo_y_nitrato

Para ello se utilizaron los siguientes paquetes de R: hydroPSO y SWATplusR.

El código utilizado para la calibración se muestra en el Anexo.

Función objetivo

La funcion objetivo se define según la siguiente ecuacion:

FO = 0.7(kge|₂₃ + kge|₃₀ + kge|₃₆ + kge|₃₈ + kge|₄₁)/5 + 0.3(kge|PT₄₁ + kge|PT₂₃)/2

Períodos de calibración validación

La siguiente tabla muestra los periodos de validación.

t23 = c("2004-01-01","2005-12-31") 
t30 = c("2005-01-01","2006-12-31") 
t36 = c("2008-01-01","2009-12-31") 
t38 = c("2007-01-01","2008-12-31") 
t41 = c("2006-01-01","2007-12-31")
T.valida = rbind(t23,t30,t36,t38,t41)
kable(T.valida, caption = "Periodos de Validación") %>%
  kable_styling(c("striped", "bordered"))

La siguiente figura muestra los períodos de validación, el resto del la serie se utiliza para calibración.

#lectura de datos
qobs <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/qobs_5est_2004_2012.RDS")

# ajuste de caudal
qobs$q23 = 1.630*qobs$q23
qobs$q30 = 0.789*qobs$q30
qobs$q36 = 1.184*qobs$q36
qobs$q38 = 0.724*qobs$q38
qobs$q41 = 1.120*qobs$q41

#conversion a mm
library(hydromad)
qobs$q23 = convertFlow(qobs$q23, from="cumecs", to="mm", area.km2 = 687)
qobs$q30 = convertFlow(qobs$q30, from="cumecs", to="mm", area.km2 = 1092)
qobs$q36 = convertFlow(qobs$q36, from="cumecs", to="mm", area.km2 = 2744)
qobs$q38 = convertFlow(qobs$q38, from="cumecs", to="mm", area.km2 = 3159)
qobs$q41 = convertFlow(qobs$q41, from="cumecs", to="mm", area.km2 = 4896)

# figura
fig = plotly_empty()
fig = fig %>% add_segments(x = ~t23[1], y=max(qobs$q23,na.rm = T), yend = max(qobs$q23,na.rm = T), 
                           xend = t23[2], name="Validación Paso Pache",
                           line = list(color = 'black', width = 2))
fig = fig %>% add_segments(x = ~t30[1], y=max(qobs$q30,na.rm = T), yend = max(qobs$q30,na.rm = T), 
                           xend = t30[2], name="Validación Paso Roldán",
                           line = list(color = 'blue', width = 2))
fig = fig %>% add_segments(x = ~t36[1], y=max(qobs$q36,na.rm = T), yend = max(qobs$q36,na.rm = T), 
                           xend = t36[2], name="Validación Fray Marcos",
                           line = list(color = 'red', width = 2))
fig = fig %>% add_segments(x = ~t38[1], y=max(qobs$q38,na.rm = T), yend = max(qobs$q38,na.rm = T),
                           xend = t38[2], name="Validación San Ramon",
                           line = list(color = 'orange', width = 2))
fig = fig %>% add_segments(x = ~t41[1], y=max(qobs$q41,na.rm = T), yend = max(qobs$q41,na.rm = T),
                           xend = t41[2], name="Validación Paso Pache",
                           line = list(color = 'green', width = 2))
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q23, type="scatter",mode = "lines",
               name="Paso de los Troncos", line = list(color = 'black', width = 2))
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q30, type="scatter",mode = "lines",
               name="Paso Roldán", line = list(color = 'blue', width = 2))
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q36, type="scatter",mode = "lines",
               name="Fray Marcos", line = list(color = 'red', width = 2))
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q38, type="scatter",mode = "lines",
               name="Paso San Ramón", line = list(color = 'orange', width = 2))
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q41, type="scatter",mode = "lines",
               name="Paso Pache", line = list(color = 'green', width = 2))
fig <- fig %>% layout(#legend = list(x = 1, y = 0.9),
                      yaxis = list(title="Runoff (mm)",
                                   showticklabels = TRUE, ticks="outside",
                                   showline=TRUE, showgrid=TRUE),
                      xaxis = list(title="",showticklabels = TRUE, 
                                   ticks="outside",showline=TRUE, 
                                   showgrid=TRUE, range=as.Date(c("2004-01-01", "2012-12-31"))),
                      height=300,
                      width=700)
fig

Parámetros a optimizar

Se consideraron 19 parámetros. El rango de variación de cada parámetro así como el tipo de cambio se muestra en la siguiente tabla:

par_all <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/par_all.RDS")
par_def = readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/par_definition.RDS")

par_upr = apply(par_all,2,max)
par_lwr = apply(par_all,2,min)
par_upr = par_upr[par_def$par_name]
par_lwr = par_lwr[par_def$par_name]

par_optim = data.frame(parametro=par_def$par_name,
                       archivo=par_def$file_name,
                       cambio=par_def$change,
                       L_inf = as.numeric(round(par_lwr,2)),
                       L_sup = as.numeric(round(par_upr,2)))


kable(par_optim, caption ="Tipo de cambio y rangos de los parámetros en la optimización") %>%
  kable_styling(c("striped", "bordered"))

Resultados de la optimización

#observados
qobs_cal <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/qobs_5est_2004_2012_para_calibracion.RDS")
qobs_cal$q23 = 1.630*qobs_cal$q23
qobs_cal$q30 = 0.789*qobs_cal$q30
qobs_cal$q36 = 1.184*qobs_cal$q36
qobs_cal$q38 = 0.724*qobs_cal$q38
qobs_cal$q41 = 1.120*qobs_cal$q41


qobs_val = readRDS("~/R_projects/SWAT_calibration_ETmodis/data/qobs_5est_2004_2012_para_validacion.RDS")

qobs_val$q23 = 1.630*qobs_val$q23
qobs_val$q30 = 0.789*qobs_val$q30
qobs_val$q36 = 1.184*qobs_val$q36
qobs_val$q38 = 0.724*qobs_val$q38
qobs_val$q41 = 1.120*qobs_val$q41



# simulados
flow23 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow23.RDS")
flow30 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow30.RDS")
flow36 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow36.RDS")
flow38 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow38.RDS")
flow41 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow41.RDS")

library(hydroGOF)

kge_stat = cbind(
kge23cal = apply(flow23[,-1], 2,function(x){KGE(x, qobs_cal$q23)}),
kge23val = apply(flow23[,-1], 2,function(x){KGE(x, qobs_val$q23)}),
kge30cal = apply(flow30[,-1], 2,function(x){KGE(x, qobs_cal$q30)}),
kge30val = apply(flow30[,-1], 2,function(x){KGE(x, qobs_val$q30)}),
kge36cal = apply(flow36[,-1], 2,function(x){KGE(x, qobs_cal$q36)}),
kge36val = apply(flow36[,-1], 2,function(x){KGE(x, qobs_val$q36)}),
kge38cal = apply(flow38[,-1], 2,function(x){KGE(x, qobs_cal$q38)}),
kge38val = apply(flow38[,-1], 2,function(x){KGE(x, qobs_val$q38)}),
kge41cal = apply(flow41[,-1], 2,function(x){KGE(x, qobs_cal$q41)}),
kge41val = apply(flow41[,-1], 2,function(x){KGE(x, qobs_val$q41)}))

pbias_stat = cbind(
  pbias23cal = apply(flow23[,-1], 2,function(x){pbias(x, qobs_cal$q23)}),
  pbias23val = apply(flow23[,-1], 2,function(x){pbias(x, qobs_val$q23)}),
  pbias30cal = apply(flow30[,-1], 2,function(x){pbias(x, qobs_cal$q30)}),
  pbias30val = apply(flow30[,-1], 2,function(x){pbias(x, qobs_val$q30)}),
  pbias36cal = apply(flow36[,-1], 2,function(x){pbias(x, qobs_cal$q36)}),
  pbias36val = apply(flow36[,-1], 2,function(x){pbias(x, qobs_val$q36)}),
  pbias38cal = apply(flow38[,-1], 2,function(x){pbias(x, qobs_cal$q38)}),
  pbias38val = apply(flow38[,-1], 2,function(x){pbias(x, qobs_val$q38)}),
  pbias41cal = apply(flow41[,-1], 2,function(x){pbias(x, qobs_cal$q41)}),
  pbias41val = apply(flow41[,-1], 2,function(x){pbias(x, qobs_val$q41)}))

rm(list = c("qobs_cal", "qobs_val"))

par(mfrow=c(2,1), mar=c(6,4,1,1))
boxplot(kge_stat, las=2, col=c(2,2,3,3,4,4,5,5,6,6), ylab="KGE", ylim=c(0,1))
boxplot(pbias_stat, las=2, col=c(2,2,3,3,4,4,5,5,6,6), ylim=c(-70,70), ylab="Bias (%)")
abline(h=c(-20,20), col=8,lty=2)

#lectura de datos
qobs <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/qobs_5est_2004_2012.RDS")

# ajuste de caudal
qobs$q23adj = 1.630*qobs$q23
qobs$q30adj = 0.789*qobs$q30
qobs$q36adj = 1.184*qobs$q36
qobs$q38adj = 0.724*qobs$q38
qobs$q41adj = 1.120*qobs$q41

# simulados
flow23 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow23.RDS")
flow30 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow30.RDS")
flow36 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow36.RDS")
flow38 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow38.RDS")
flow41 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow41.RDS")

est = "23"
##
Rdate = qobs$Rdate
qo = qobs[,paste0("q",est)]
qa = qobs[,paste0("q",est,"adj")]
qs = readRDS(paste0("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow",est,".RDS"))

qh = apply(qs[,-1], 1,function(x){max(x)})
ql = apply(qs[,-1], 1,function(x){min(x)})
qm = apply(qs[,-1], 1,function(x){median(x)})

fig <- plotly_empty()
fig = fig %>% add_polygons(x=c(Rdate,rev(Rdate)),
                           y=c(ql,rev(qh)), name="Qsim_range",
                           line=list(width=2,color="gray"),
                           fillcolor='gray', inherit = FALSE, opacity = 0.5)  

fig <- fig %>% add_lines(x = ~Rdate, y = ~qm, name="Qsim_median",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'orange', width = 2))

fig <- fig %>% add_lines(x = ~Rdate, y = ~qo, name="Qobs",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'black', width = 2))


fig <- fig %>% add_lines(x = ~Rdate, y = ~qa, name="Qobs_adjusted",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'red', width = 2))

fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9),
                      yaxis = list(title="Flow (m3/s)",showticklabels = TRUE, ticks="outside",showline=TRUE, showgrid=TRUE),
                      xaxis = list(title="",showticklabels = TRUE, ticks="outside",showline=TRUE, showgrid=TRUE))

fig

est = "30"

Rdate = qobs$Rdate
qo = qobs[,paste0("q",est)]
qa = qobs[,paste0("q",est,"adj")]
qs = readRDS(paste0("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow",est,".RDS"))

qh = apply(qs[,-1], 1,function(x){max(x)})
ql = apply(qs[,-1], 1,function(x){min(x)})
qm = apply(qs[,-1], 1,function(x){median(x)})

fig <- plotly_empty()
fig = fig %>% add_polygons(x=c(Rdate,rev(Rdate)),
                           y=c(ql,rev(qh)), name="Qsim_range",
                           line=list(width=2,color="gray"),
                           fillcolor='gray', inherit = FALSE, opacity = 0.5)  

fig <- fig %>% add_lines(x = ~Rdate, y = ~qm, name="Qsim_median",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'orange', width = 2))

fig <- fig %>% add_lines(x = ~Rdate, y = ~qo, name="Qobs",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'black', width = 2))


fig <- fig %>% add_lines(x = ~Rdate, y = ~qa, name="Qobs_adjusted",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'red', width = 2))

fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9),
                      yaxis = list(title="Flow (m3/s)",showticklabels = TRUE, ticks="outside",showline=TRUE, showgrid=TRUE),
                      xaxis = list(title="",showticklabels = TRUE, ticks="outside",showline=TRUE, showgrid=TRUE))

fig

est = "36"
##
Rdate = qobs$Rdate
qo = qobs[,paste0("q",est)]
qa = qobs[,paste0("q",est,"adj")]
qs = readRDS(paste0("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow",est,".RDS"))

qh = apply(qs[,-1], 1,function(x){max(x)})
ql = apply(qs[,-1], 1,function(x){min(x)})
qm = apply(qs[,-1], 1,function(x){median(x)})

fig <- plotly_empty()
fig = fig %>% add_polygons(x=c(Rdate,rev(Rdate)),
                           y=c(ql,rev(qh)), name="Qsim_range",
                           line=list(width=2,color="gray"),
                           fillcolor='gray', inherit = FALSE, opacity = 0.5)  

fig <- fig %>% add_lines(x = ~Rdate, y = ~qm, name="Qsim_median",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'orange', width = 2))

fig <- fig %>% add_lines(x = ~Rdate, y = ~qo, name="Qobs",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'black', width = 2))


fig <- fig %>% add_lines(x = ~Rdate, y = ~qa, name="Qobs_adjusted",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'red', width = 2))

fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9),
                      yaxis = list(title="Flow (m3/s)",showticklabels = TRUE, ticks="outside",showline=TRUE, showgrid=TRUE),
                      xaxis = list(title="",showticklabels = TRUE, ticks="outside",showline=TRUE, showgrid=TRUE))

fig

est = "38"
##
Rdate = qobs$Rdate
qo = qobs[,paste0("q",est)]
qa = qobs[,paste0("q",est,"adj")]
qs = readRDS(paste0("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow",est,".RDS"))

qh = apply(qs[,-1], 1,function(x){max(x)})
ql = apply(qs[,-1], 1,function(x){min(x)})
qm = apply(qs[,-1], 1,function(x){median(x)})

fig <- plotly_empty()
fig = fig %>% add_polygons(x=c(Rdate,rev(Rdate)),
                           y=c(ql,rev(qh)), name="Qsim_range",
                           line=list(width=2,color="gray"),
                           fillcolor='gray', inherit = FALSE, opacity = 0.5)  

fig <- fig %>% add_lines(x = ~Rdate, y = ~qm, name="Qsim_median",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'orange', width = 2))

fig <- fig %>% add_lines(x = ~Rdate, y = ~qo, name="Qobs",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'black', width = 2))


fig <- fig %>% add_lines(x = ~Rdate, y = ~qa, name="Qobs_adjusted",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'red', width = 2))

fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9),
                      yaxis = list(title="Flow (m3/s)",showticklabels = TRUE, ticks="outside",showline=TRUE, showgrid=TRUE),
                      xaxis = list(title="",showticklabels = TRUE, ticks="outside",showline=TRUE, showgrid=TRUE))

fig

est = "36"
##
Rdate = qobs$Rdate
qo = qobs[,paste0("q",est)]
qa = qobs[,paste0("q",est,"adj")]
qs = readRDS(paste0("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow",est,".RDS"))

qh = apply(qs[,-1], 1,function(x){max(x)})
ql = apply(qs[,-1], 1,function(x){min(x)})
qm = apply(qs[,-1], 1,function(x){median(x)})

fig <- plotly_empty()
fig = fig %>% add_polygons(x=c(Rdate,rev(Rdate)),
                           y=c(ql,rev(qh)), name="Qsim_range",
                           line=list(width=2,color="gray"),
                           fillcolor='gray', inherit = FALSE, opacity = 0.5)  

fig <- fig %>% add_lines(x = ~Rdate, y = ~qm, name="Qsim_median",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'orange', width = 2))

fig <- fig %>% add_lines(x = ~Rdate, y = ~qo, name="Qobs",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'black', width = 2))


fig <- fig %>% add_lines(x = ~Rdate, y = ~qa, name="Qobs_adjusted",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'red', width = 2))

fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9),
                      yaxis = list(title="Flow (m3/s)",showticklabels = TRUE, ticks="outside",showline=TRUE, showgrid=TRUE),
                      xaxis = list(title="",showticklabels = TRUE, ticks="outside",showline=TRUE, showgrid=TRUE))

fig

#lectura de datos
qobs <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/qobs_5est_2004_2012.RDS")

# ajuste de caudal
qobs$q23adj = 1.630*qobs$q23
qobs$q30adj = 0.789*qobs$q30
qobs$q36adj = 1.184*qobs$q36
qobs$q38adj = 0.724*qobs$q38
qobs$q41adj = 1.120*qobs$q41

qlt_41 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/rch41_calidad.RDS")
#qlt_41$qobs = 1.12*qlt_41$qobs
qlt_23 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/rch23_calidad.RDS")
#qlt_23$qobs = 1.66*qlt_23$qobs

# simulados
flow23 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow23.RDS")
flow41 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/flow41.RDS")

library(hydroGOF)

PT23 = readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/pt23.RDS")
PT41 = readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_27/pt23.RDS")

PTvsFlow41 = matrix(ncol=2,nrow = 0)
PTvsFlow23 = matrix(ncol=2,nrow = 0)

for(i in 2:ncol(flow23)){
  PTvsFlow41 = rbind(PTvsFlow41,
                     cbind(PT = 1000000*unlist(PT41[,i])/(24*60*60*4896),
                           Flow = 1000*unlist(flow41[,i])/4896))
  PTvsFlow23 = rbind(PTvsFlow23,
                     cbind(PT = 1000000*unlist(PT23[,i])/(24*60*60*687),
                           Flow = 1000*unlist(flow23[,i])/687))
}

PTvsFlow23[PTvsFlow23<=0] = NA
PTvsFlow23 = na.omit(PTvsFlow23)

PTvsFlow41[PTvsFlow41<=0] = NA
PTvsFlow41 = na.omit(PTvsFlow41)

par(mfrow=c(1,2), mar=c(4,4,1,1))

plot(PTvsFlow23[,"Flow"], PTvsFlow23[,"PT"], log="xy",col=8,pch="+", 
     xlab = "flow (lps/km2)", ylab="PT (mg/s/km2)", main="Paso de los Troncos",
          ylim=c(1e-10, 1e+3), xlim=c(1e-4,1e+3))
lines(1000*qobs$q23adj/687, 0.06*(1000*qobs$q23adj/687)^0.93, col=2, lwd=2)
points(1000*qlt_23$qobs/687, qlt_23$PT_mgxs/687, pch=19, col=2,cex=1.5)

fit = lm(log(PTvsFlow23[,"PT"])~log(PTvsFlow23[,"Flow"]))
fit = c(exp(fit$coefficients[1]), fit$coefficients[2])
fit = fit[1]*(PTvsFlow23[,"Flow"])^fit[2]
lines(PTvsFlow23[,"Flow"], fit, col=1,lwd=2)


plot(PTvsFlow41[,"Flow"], PTvsFlow41[,"PT"], log="xy",col=8,pch="+", 
     xlab = "flow (lps/km2)", ylab="PT (mg/s/km2)", main="Paso Pache", 
     ylim=c(1e-10, 1e+3), xlim=c(1e-4,1e+3))
lines(1000*qobs$q41adj/4896, 0.16*(1000*qobs$q41adj/4896)^1.02, col=2, pch=19)
points(1000*qlt_41$qobs/4896, qlt_41$PT_mgxs/4896, pch=19, col=2,cex=1.5)

fit = lm(log(PTvsFlow41[,"PT"])~log(PTvsFlow41[,"Flow"]))
fit = c(exp(fit$coefficients[1]), fit$coefficients[2])
fit = fit[1]*(PTvsFlow41[,"Flow"])^fit[2]
lines(PTvsFlow41[,"Flow"], fit, col=1,lwd=2)

Anexo

rm(list = ls())

# Librerias
library(SWATplusR)
library(stringr)
library(hydroGOF)
library(hydromad)
library(tibble)
library(hydroPSO)

# Caudales observados
obs.data <- readRDS("./data/qobs_5est_2004_2012_para_calibracion.RDS")

# ajuste de caudal
obs.data$q23 = 1.630*obs.data$q23
obs.data$q30 = 0.789*obs.data$q30
obs.data$q36 = 1.184*obs.data$q36
obs.data$q38 = 0.724*obs.data$q38
obs.data$q41 = 1.120*obs.data$q41


# SWAT project mensual
swat_project = "./SWAT_SL_20200401/"

# Funcion de optimización
swat_optim <- function(mypar, obs.data, project_dir, start_sim, end_sim,runrun_path, head_log) {
  
  
  names(mypar) = head_log[25:length(head_log)]
  path_run = paste0(runrun_path,"run_",Sys.getpid())
  
  par_set = tibble("CN2::CN2.mgt|change = pctchg" = mypar["CN2"]*100,
                   "ESCO::ESCO.hru|change = absval" = mypar["ESCO"],
                   "SOL_AWC::SOL_AWC.sol|change = pctchg " = mypar["SOL_AWC"]*100,
                   "ALPHA_BF::ALPHA_BF.gw|change = absval" = mypar["ALPHA_BF"],
                   "GWQMN::GWQMN.gw|change = absval" = mypar["GWQMN"],
                   "GW_DELAY::GW_DELAY.gw|change = absval |" = mypar["GW_DELAY"],
                   "REVAPMN::REVAPMN.gw|change = absval" =  mypar["REVAPMN"],
                   "GW_REVAP::GW_REVAP.gw|change = absval" = mypar["GW_REVAP"],
                   "OV_N::OV_N.hru|change = absval" = mypar["OV_N"],
                   "SLSUBBSN::SLSUBBSN.hru|change = pctchg" = mypar["SLSUBBSN"]*100,
                   
                   "HRU_SLP::HRU_SLP.hru|change= pctchg" = mypar["HRU_SLP"]*100,
                   "P_UPDIS::P_UPDIS.bsn|change= absval"= mypar["P_UPDIS"], 
                   "FRT_SURFACE::FRT_SURFACE.mgt|change=absval" = mypar["FRT_SURFACE"], 
                   "FRT_KG::FRT_KG.mgt|change=pctchg" = mypar["FRT_KG"]*100,
                   "PPERCO::PPERCO.bsn|change= absval"= mypar["PPERCO"],                    
                   "CMN::CMN.bsn|change= absval"= mypar["CMN"],                    
                   "PSP::PSP.bsn|change= absval"= mypar["PSP"],
                   "PHOSKD::PHOSKD.bsn|change= absval"= mypar["PHOSKD"],
                   "RSDCO::RSDCO.bsn|change= absval"= mypar["RSDCO"]
  ) 
  
  
  # run swat
  qsim = run_swat2012(project_path = project_dir,
                      output = list(FLOW = define_output(file = "rch",
                                                         variable = "FLOW_OUT",
                                                         unit = c(23,30,36,38,41)),
                                    PT = define_output(file = "rch",
                                                       variable = "TOT_Pkg",
                                                       unit = c(23,41)),
                                    NO3 = define_output(file = "rch",
                                                        variable = "NO3_OUT",
                                                        unit = c(23,41))),
                      parameter = par_set,
                      start_date = as.Date("2000-01-01"),
                      end_date = as.Date("2012-12-31"),
                      years_skip = 4,
                      run_path=path_run,
                      quiet = TRUE,
                      output_interval="d",
                      add_parameter = FALSE)
  
  # objective function
  nse23 = NSE(as.vector(qsim$FLOW_23), as.vector(obs.data$q23))
  nse30 = NSE(as.vector(qsim$FLOW_30), as.vector(obs.data$q30))
  nse36 = NSE(as.vector(qsim$FLOW_36), as.vector(obs.data$q36))
  nse38 = NSE(as.vector(qsim$FLOW_38), as.vector(obs.data$q38))
  nse41 = NSE(as.vector(qsim$FLOW_41), as.vector(obs.data$q41))
  
  kge23 = KGE(as.vector(qsim$FLOW_23), as.vector(obs.data$q23))
  kge30 = KGE(as.vector(qsim$FLOW_30), as.vector(obs.data$q30))
  kge36 = KGE(as.vector(qsim$FLOW_36), as.vector(obs.data$q36))
  kge38 = KGE(as.vector(qsim$FLOW_38), as.vector(obs.data$q38))
  kge41 = KGE(as.vector(qsim$FLOW_41), as.vector(obs.data$q41))
  
  pbias23 = pbias(as.vector(qsim$FLOW_23), as.vector(obs.data$q23))
  pbias30 = pbias(as.vector(qsim$FLOW_30), as.vector(obs.data$q30))
  pbias36 = pbias(as.vector(qsim$FLOW_36), as.vector(obs.data$q36))
  pbias38 = pbias(as.vector(qsim$FLOW_38), as.vector(obs.data$q38))
  pbias41 = pbias(as.vector(qsim$FLOW_41), as.vector(obs.data$q41))
  
  
  # objective quality
  qlt_41 = data.frame(flow_lps_km2=1000*qsim$FLOW_41/4896,
                      PT_mgs_km2 = 1000000*qsim$PT_41/(24*60*60*4896),
                      NO3_mgs_km2 = 1000000*qsim$NO3_41/(24*60*60*4896),
                      qobs_lps_km2 = 1000*obs.data$q41/4896,
                      PTobs_mgs_km2 = 0.16*(1000*obs.data$q41/4896)^1.02,
                      NO3obs_mgs_km2 = 0.1*(1000*obs.data$q41/4896)^1.41)
  
  qlt_41[qlt_41<=0] = NA
  qlt_41 = na.omit(qlt_41)
  
  qlt_23 = data.frame(flow_lps_km2=1000*qsim$FLOW_23/687,
                      PT_mgs_km2 = 1000000*qsim$PT_23/(24*60*60*687),
                      NO3_mgs_km2 = 1000000*qsim$NO3_23/(24*60*60*687),
                      qobs_lps_km2 = 1000*obs.data$q23/687,
                      PTobs_mgs_km2 = 0.06*(1000*obs.data$q23/687)^0.93,
                      NO3obs_mgs_km2 = 0.04*(1000*obs.data$q23/687)^1.41)
  
  qlt_23[qlt_23<=0] = NA
  qlt_23 = na.omit(qlt_23)
  
  fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt_41)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_41$PT_fit = a*(qlt_41$flow_lps_km2)^b
  
  fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt_41)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_41$NO3_fit = a*(qlt_41$flow_lps_km2)^b
  
  fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt_23)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_23$PT_fit = a*(qlt_23$flow_lps_km2)^b
  
  fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt_23)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_23$NO3_fit = a*(qlt_23$flow_lps_km2)^b
  
  kge_PT_41 = KGE(qlt_41$PT_fit, qlt_41$PTobs_mgs_km2)
  kge_NO3_41 = KGE(qlt_41$NO3_fit, qlt_41$NO3obs_mgs_km2)
  kge_PT_23 = KGE(qlt_23$PT_fit, qlt_23$PTobs_mgs_km2)
  kge_NO3_23 = KGE(qlt_23$NO3_fit, qlt_23$NO3obs_mgs_km2)
  
  pbias_PT_23 = pbias(qlt_23$PT_fit, qlt_23$PTobs_mgs_km2)
  pbias_PT_41 = pbias(qlt_41$PT_fit, qlt_41$PTobs_mgs_km2)
  
  pbias_NO3_23 = pbias(qlt_23$NO3_fit, qlt_23$NO3obs_mgs_km2)
  pbias_NO3_41 = pbias(qlt_41$NO3_fit, qlt_41$NO3obs_mgs_km2)
  
  FO = 0.7*(kge23 + kge30 + kge36 + kge38 + kge41)/5 + 0.3*(kge_PT_41 + kge_PT_23)/2 
  
  # logfile
  vout = data.frame(FO,kge23,kge30,kge36,kge38,kge41,
                    nse23,nse30,nse36,nse38,nse41,
                    pbias23,pbias30,pbias36,pbias38,pbias41,
                    kge_PT_23, kge_PT_41, kge_NO3_41, kge_NO3_23,
                    pbias_PT_23, pbias_PT_41, pbias_NO3_23, pbias_NO3_41,
                    t(mypar))
  colnames(vout) = head_log
  log.file = paste0(runrun_path,"log_",Sys.getpid(),".txt")
  
  if(!any(list.files(runrun_path)==paste0("log_",Sys.getpid(),".txt"))){
    write.table(t(head_log), log.file, col.names = FALSE, row.names = FALSE, quote = FALSE)
  }
  
  ff = read.table(log.file, col.names = head_log,colClasses = "character")
  ff = rbind(ff,vout)
  write.table(ff, log.file, col.names = FALSE, row.names = FALSE, quote = FALSE)
  
  #output
  return(FO)
}

CN2 = c(-0.2,0.2)
ESCO = c(0.1,0.45)
SOL_AWC = c(-0.2,.2)
ALPHA_BF = c(0.4,0.8)  
GWQMN = c(800,7000)
GW_DELAY = c(2,45)
REVAPMN = c(300,700)
GW_REVAP = c(0.02, 0.01)
OV_N = c(0.4, 0.6)
SLSUBBSN = c(-0.1, 0.3)

HRU_SLP = c(-0.1,0.1)
P_UPDIS = c(0, 150)
FRT_SURFACE = c(0.1,1)
FRT_KG = c(-0.3,0.5)
PPERCO = c(10,17.5) #bsn absval
CMN = c(0,0.03) #bsn absval
PSP = c(0.01,0.8)
PHOSKD = c(10,1000)
RSDCO = c(0.01, 0.1)

names_par = c("CN2","ESCO","SOL_AWC",
              "ALPHA_BF","GWQMN","GW_DELAY","REVAPMN","GW_REVAP","OV_N","SLSUBBSN","HRU_SLP",
              "P_UPDIS","FRT_SURFACE","FRT_KG","PPERCO","CMN","PSP", "PHOSKD", "RSDCO")

#           CN2   ESCO  SOL_AWC ALPHA_BF  GWQMN GW_DELAY  REVAPMN GW_REVAP  OV_N  SLSUBBSN
par_lwr = c(CN2[1],
            ESCO[1], SOL_AWC[1], ALPHA_BF[1], GWQMN[1], GW_DELAY[1], REVAPMN[1], GW_REVAP[1], OV_N[1],SLSUBBSN[1],
            HRU_SLP[1], P_UPDIS[1], FRT_SURFACE[1], FRT_KG[1], PPERCO[1], CMN[1],PSP[1], PHOSKD[1], RSDCO[1])

par_upr = c(CN2[2],
            ESCO[2], SOL_AWC[2], ALPHA_BF[2], GWQMN[2], GW_DELAY[2], REVAPMN[2], GW_REVAP[2], OV_N[2],SLSUBBSN[2],
            HRU_SLP[2], P_UPDIS[2], FRT_SURFACE[2], FRT_KG[2], PPERCO[2], CMN[2],PSP[2], PHOSKD[2], RSDCO[2])


par_init = (par_upr + par_lwr)/2


# Optimización
head_log = c("FO","KGE23","KGE30","KGE36", "KGE38","KGE41",
             "NSE23","NSE30","NSE36","NSE38", "NSE41",
             "pbias23","pbias30","pbias36","pbias38", "pbias41",
             "kge_PT_23", "kge_PT_41", "kge_NO3_41", "kge_NO3_23",
             "pbias_PT_23", "pbias_PT_41", "pbias_NO3_23", "pbias_NO3_41",
             names_par)


fr = list.files("./run_optim_mensual/", pattern = "*.txt", full.names=TRUE)
file.remove(fr)
opt_pso = hydroPSO(par=par_init, fn = swat_optim,lower=par_lwr, upper=par_upr, 
                   control = list(MinMax = "max",
                                  maxit=500,
                                  reltol=0.00001,
                                  parallel="parallelWin", 
                                  par.nnodes=12,
                                  par.pkgs = list("hydromad", 
                                                  "SWATplusR",
                                                  "stringr",
                                                  "hydroGOF",
                                                  "tibble"),
                                  write2disk=FALSE,
                                  REPORT=1,
                                  npart=36,
                                  normalise=TRUE),
                   obs.data = obs.data, 
                   project_dir = swat_project, 
                   runrun_path = "./run_optim_mensual/",
                   head_log = head_log
)