Última actualización: 03 junio 2020



1 Escenarios

Se definieron 26 escenarios utilizando diferentes anchos de la zona buffer (5, 10, 20, 40, 80 metros). Los 26 escenarios surgen de la combinación por zona y por tipo de cauce. Se definen tres zonas en la cuenca en función de la capacidad que tiene cada región geográfica para producir sedimentos y nutrientes ver Zonas. Se consideran dos tipos de cauces según la norma (cauces principales y tributarios).

Los escenarios quedan definidos de la siguiente forma:

  • Tipo “ABC” variación uniforme en toda la cuenca (6 escenarios).
  • Tipo “BC” variación de la zona “BC”. La zona “A” se fija según la norma (5 escenarios).
  • Tipo “A”, “B”, o “C” variacion de las zonas “A”, “B” o “C”, en el resto de la cuenca se fija la norma (“BC”, “AC”, “AB”, respectivamente) (15 escenarios).
Tamaños de la zona buffer para cada escenario
escenario BP_a BP_b BP_c BT_a BT_b BT_c
ABC_1 80 80 80 40 40 40
ABC_2 40 40 40 20 20 20
ABC_3 20 20 20 10 10 10
ABC_4 10 10 10 5 5 5
ABC_5 40 40 40 10 10 10
ABC_6 40 40 40 5 5 5
A1 80 40 40 40 20 20
A3 20 40 40 10 20 20
A4 10 40 40 5 20 20
A5 40 40 40 10 20 20
A6 40 40 40 5 20 20
B1 40 80 40 20 40 20
B3 40 20 40 20 10 20
B4 40 10 40 20 5 20
B5 40 40 40 20 10 20
B6 40 40 40 20 5 20
C1 40 40 80 20 20 40
C3 40 40 20 20 20 10
C4 40 40 10 20 20 5
C5 40 40 40 20 20 10
C6 40 40 40 20 20 5
BC_1 40 80 80 20 40 40
BC_3 40 20 20 20 10 10
BC_4 40 10 10 20 5 5
BC_5 40 40 40 20 10 10
BC_6 40 40 40 20 5 5

2 Zonas

En función de los resultados de la simulación (producción de sedimentos y fósforo por subcuenca), se identificaron 3 zonas contrastantes. La zona A, tiene la menor producción de sedimentos, se ubica en la parte alta. La zona B, tiene la mayor producción de sedimentos, se ubica en la vertiente norte del cauce principal. La zona C, presenta una producción de sedimentos inferior a la zona B y significativamente mayor que la zona A.

3 Resultados

Se computaron los escenarios en el servidor de INIA. El tiempo de simulación es alrededor de 1 hora por cada escenario (12 procesadores). Cada escenario tiene 100 realizaciones del modelo. Para disminuir el volumen de datos se guarda únicamente las distribuciones de probabilidad acumulada (mínima, media y máxima) de la concentración de fósforo en cada uno de los 41 canales.

La siguiente figura muestra la probabilidad de no excedencia de la concentración de fósforo para el escenario base (situación actual) y los escenarios de zona buffer 40-20 y 10-05 (Principal-tributario en metros). Ambos escenarios de zona buffer producen resultados bastante similares. Esto será realmente así o hay algún proceso que no se esté representando bien por el modelo?

Anexo: código para correr los escenarios

rm(list=ls())

library(tibble)
library(SWATplusR)
library(stringi)
library(stringr)

# lectura de parámetros optimizados
par_set = readRDS("./input_buffer/best_par_tibble.RDS")

# Lectura de filter ratio para todos los escenarios
FR = readRDS("./input_buffer/03_filter_ratio.RDS")

#Lectura metadata HRU
hru_info <- readRDS("./input_buffer/hru_info.RDS")

model_path = "./SWAT_SL_20200401/"
runrun_path = "./run_buffer"

# escenarios a modelar
escenarios_buffer = readRDS("./input_buffer/escenarios_buffer.RDS")

# Zonificación
zonaA = c(1,7,8,10,14,15,16,17,18,19,20,22,24,28)
zonaB = c(4,5,25,31,3,27,2,23,29)
zonaC = c(41,12,21,30,9,13,32,38,36,34,6,11,37,35,39,40,26,33)



## Inicio de loop en los escenarios
########################################

for(n in 1:nrow(escenarios_buffer)){

# Cuidado aqui se elimina la carpeta runrun_path!!!!!!!!!!!!!!!!
unlink(runrun_path, recursive = TRUE)
dir.create(runrun_path)

# creacion de un nuevo directorio
R.utils::copyDirectory(from=model_path, to=runrun_path)

# Seleccion de HRU a aplicar cambios (FILTER_RATIO) 
hru2change = hru_info[hru_info$LANDUSE %in% c("GRAS", "AGRC", "AGRP", "LECH"),]

# Alteracion de parametros del Filter Strip
cc =  as.numeric(paste(FR$ZB_principal_km2))==escenarios_buffer[n,"BP_a"] & 
  as.numeric(paste(FR$ZB_tributarios_km2))==escenarios_buffer[n,"BT_a"] & FR$subassin%in%zonaA |
  as.numeric(paste(FR$ZB_principal_km2))==escenarios_buffer[n,"BP_b"] & 
  as.numeric(paste(FR$ZB_tributarios_km2))==escenarios_buffer[n,"BT_b"] & FR$subassin%in%zonaB |
  as.numeric(paste(FR$ZB_principal_km2))==escenarios_buffer[n,"BP_c"] & 
  as.numeric(paste(FR$ZB_tributarios_km2))==escenarios_buffer[n,"BT_c"] & FR$subassin%in%zonaC

FilterRatio2sub = FR[cc,]

hru2change$FilterRatio = FilterRatio2sub$FILTER_RATIO[match(hru2change$Subbasin, FilterRatio2sub$subassin)]

# creacion de archivos *.ops
for(i in 1:nrow(hru2change)){
  MONTH = 1
  DAY = 1
  IYEAR = 2009
  MGT_OP = 4
  FILTER_I = 1
  FILTER_RATIO = hru2change$FilterRatio[i]
  FILTER_CON = 0.5
  FILTER_CH = 0
  
  fs_line =  paste0(str_pad(MONTH, 3, pad=" "),
                    str_pad(DAY, 6-3, pad=" "),
                    str_pad(IYEAR, 15-6, pad=" "),
                    str_pad(MGT_OP,18-15, pad = " "),
                    str_pad(FILTER_I,23-18, pad = " "),
                    str_pad(round(FILTER_RATIO,2),34-23, pad = " "),
                    str_pad(round(FILTER_CON,5),47-34, pad = " "),
                    str_pad(round(FILTER_CH,2),54-47, pad = " "))
  
  ops_file = paste0(hru2change$HRUGIS[i],".ops")
  ops_out = rbind(".ops generado en R", fs_line)
  write.table(ops_out, file=paste0(runrun_path,"/", ops_file), quote = F, col.names = F,
              row.names = F)
}

#Modificacion de archivos *.sub
sub_files = list.files(runrun_path, pattern = ".sub", full.names = T)
for(i in 1:length(sub_files)){
  sub_file = readLines(sub_files[i])
  cc = substr(sub_file,1,9) %in% paste(hru2change$HRUGIS)    
  substr(sub_file[cc], 66, 78) <-  paste0(substr(sub_file[cc],1,10),"ops")
  write.table(sub_file, file=sub_files[i], quote = F, col.names = F,
              row.names = F)
}


# ejecucion del modelo
buffer_run = run_swat2012(project_path = runrun_path,
                          output = list(FLOW = define_output(file = "rch",
                                                             variable = "FLOW_OUT",
                                                             unit = 1:41),
                                        PT = define_output(file = "rch",
                                                           variable = "TOT_Pkg",
                                                           unit = 1:41),
                                        SYLD = define_output(file = "rch",
                                                             variable = "SED_OUT",
                                                             unit = 1:41)),
                          parameter = par_set,
                          start_date = as.Date("2005-01-01"),
                          end_date = as.Date("2015-12-31"),
                          years_skip = 1,
                          add_parameter = FALSE,
                          n_thread = 12,
                          output_interval = "m")

# output files:
# 1 archivo por variable y por escenario
# contiene 100 realizaciones en cada canal a escala mensual

flow = buffer_run[1:41]
pt = buffer_run[42:82]
sed = buffer_run[83:123]


saveRDS(flow,paste0("./output_buffer/flow_",escenarios_buffer[n,"escenario"],".RDS"))
saveRDS(pt,paste0("./output_buffer/pt_",escenarios_buffer[n,"escenario"],".RDS"))
saveRDS(sed,paste0("./output_buffer/sed_",escenarios_buffer[n,"escenario"],".RDS"))

}
 



A work by Rafael Navas