Estimar parámetros de la curva de calibración de caudales o descarga con R

#################################################################################
#
#
# Estimación de los parámetros de l curva de descarga
# Ajusta y plotea la grafica de curva descarga siguiendo una funcion exponencial
#
#  Walter Arnoldo Bardales Espinoza
#
#################################################################################

# remover todos los objetos de almacenamiento
rm(list = ls())

# cargar libreria
require(Hmisc)
library(minpack.lm)

# leer datos de aforos y niveles
dato.aforos = read.table("D:/Documents/Curva_calibracion_de_caudales/est_Matucuy.txt", header = TRUE, sep = "\t",na.strings = "-99.9")


plot(Caudal~Nivel)

# Para usar los nombres de los campos de las columnas
attach(dato.aforos)


# Estima los valores de los parametros iniciales de la curva
# Estableciendo h0 a un valor por debajo del minimo observado
rango=range(Nivel)
h.min=rango[1]
h.max=rango[2]
ho.inicial =  h.min - 0.1*(h.max - h.min)

# modelo lineal para estimar los valores inicial de k y n
modelo.log = lm( log10(Caudal) ~ log10(Nivel - ho.inicial) )
a.inicial = 10^modelo.log$coefficients[[1]]
b.inicial = modelo.log$coefficients[[2]]

# Utilizar los coeficientes calculados inicialmente para estimar los parametros de la regresion no lineal
curva.descarga = nlsLM( Caudal ~ k*(Nivel-ho)^n, start = list(k=a.inicial,n=b.inicial, ho=ho.inicial), control = list(maxiter=500))

# Resumen de los estadisticos de la curva de descarga
summary(curva.descarga)

# Grafico de la curva de descarga o calibracion
xmin = min(Nivel)
xmax = max(Nivel)
dx = 0.01*(xmax - xmin)
x = seq(xmin,xmax,dx)
coef.model = coef(curva.descarga)
qest = coef.model[1]*(x - coef.model[3])^coef.model[2]


plot(Nivel, Caudal,
     ylim = c(0, max(Caudal)+0.2*(max(Caudal)-min(Caudal))),
     xlim = c(min(Nivel)-0.1*(rango[2]-rango[1]), max(Nivel)+0.2*(rango[2]-rango[1])),
     xlab = "Nivel (m)",
     ylab = expression("Q ("*m^3*s^{-1}*")"),
     pch = 21, bg = "black",
     main = "Curva de calibracion de caudales de la estacion Matucuy, Periodo 2002 - 2015"
)
lines(x, qest, col = "red")
text(min(Nivel)+0.1*(max(Nivel)-min(Nivel)),max(Caudal), paste0("Q = ",round(coef.model[1],3),"*(H-(",round(coef.model[3],3),"))^",round(coef.model[2],3)))
text(min(Nivel)+0.1*(max(Nivel)-min(Nivel)),0.95*max(Caudal),paste0("RMS =", round(mean((Caudal-(coef.model[1]*(Nivel - coef.model[3])^coef.model[2]))^2))))

#################################################################### Fin código

Descargar datos y scripts R

Prueba de infiltración

Prueba de infiltración y ajuste con el modelo de Horton o Kostiakov

Infiltración (I):

La teoría de infiltración está basada en un análisis del movimiento del agua en el suelo bajo condiciones de no saturación. Durante la infiltración, la fase líquida y la fase gaseosa coexisten en la masa de suelo, con excepción de la zona de contacto entre el suelo y el agua en la superficie del terreno.

La infiltración, velocidad de infiltración o intensidad de entrada en el suelo, se puede definir como la velocidad de penetración del agua en el perfil del suelo, cuando la superficie del terreno está cubierta por una capa de agua poco profunda.

La infiltración tiene dimensión de velocidad (L T^-1), como la lámina de agua (L) admitida por el suelo en una unidad de tiempo (T) o como la cantidad de agua absorbida por la unidad de superficie del terreno en la unidad de tiempo (L³ T^-1 L^-2), respectivamente. Si las mismas unidades se usan en ambos casos, las expresiones son dimensionalmente equivalentes (L T^-1). En la primera forma la expresión común de velocidad de entrada es cm h^-1 o cm min^-1.  En la segunda forma, generalmente se expresa como m³ min^-1 m^-2. Es común, por varios autores, señalar que la infiltración obedece a un fenómeno de movimiento vertical del agua en el perfil de suelo.  A la disminución de la velocidad de infiltración después de cierto período de tiempo alcanzando una velocidad que tiende a ser constante y se denomina infiltración básica (I_b), la cual se considera cercana a la conductividad hidráulica o permeabilidad (K).

Factores que afectan la infiltración

·       Textura del suelo: Proporción de arenas, limos y arcilla

·       Estructura del suelo: Es la forma en la que se encuentra posicionadas las partículas del suelo o pedón (granular, migajosa, bloques, prismática o columnar o laminar)

·       Contenido de humedad del suelo (Capacidad de campo, punto de marchitez permanente y saturación)

Métodos para medir la velocidad de infiltración:

Las técnicas más utilizadas para determinar tasas de infiltración y conductividad hidráulica saturada, son: Lagunas de infiltración, permeámetros, infiltrómetros cilíndricos y el método de Porchet.

Porchet

Para la estimación de la tasa de infiltración en terreno se puede utilizar el método de Porchet, el cual consiste en excavar un cilindro de radio (R) y se llenarlo con agua hasta una altura (h).

La superficie por la cual pasa el agua al infiltrarse es:

S=2π*R*h+ π*R^2

S=πR(2h+R)

La velocidad de infiltración (f) que pasa por la superficie (S) del agujero cilíndrico es proporcional a la variación del volumen de agua del agujero cilíndrico para un intervalo de tiempo infinitamente pequeño (dt):

S*f =  dV/dt

Como el radio permanece constante, la expresión se puede simplificar de la siguiente forma:

f=  (πR^2)/S dh/dt

Sustituir la superficie del agujero cilíndrico

f=  (πR^2)/(πR(2h+R))  dh/dt

Al simplificar

f=  R/(2h+R)  dh/dt

Se parar los diferenciales

f dt=  R/((2h+R))  dh

Al integrar y valuar en los intervalos de tiempo (t1, t2) y el nivel de agua (h1 y h2), siendo t12

y h21

f (t2-t1)=  R/2*(Ln(2*h1+R)- Ln(2*h2+R))

Al simplificar

f=  R/(2(t2-t1))*Ln((2*h1+R)/(2*h2+R))

Métodos de ajuste de la ecuación de infiltración

Horton:

          Horton estudio la capacidad de infiltración, en su estudio determinó que la capacidad de infiltración tiene la siguiente forma.

fp= fc + (fo-fc )*e^(-kt)  

Donde:

fp= Capacidad de infiltración (mm/hr)

fc= Capacidad final o equilibrio (mm/hr)
fo= Capacidad inicial (mm/hr)

K= Coeficiente de decrecimiento de la capacidad de infiltración en el tiempo

Kostiakov

          Este investigador propuso una ecuación empírica, basada en un modelo potencial

f=a*t^b

Donde:

f= velocidad de infiltración (mm/min)

a y b = parámetros de ajuste

t= Tiempo (min)

Modelo HBV en excel

Introducción:

El modelo HBV fue desarrollado originalmente por Integrated Hidrological Modelling Sistem (IHMS) a principios de los años 70 para ayudar a las operaciones de energía hidroeléctrica. El enfoque del modelo HBV ha demostrado ser flexible y robusto en la solución de los problemas de los recursos hídricos y aplicaciones que abarcan una amplia gama hasta ahora. El modelo es la herramienta estándar de predicción en Suecia, donde alrededor de 75 cuencas, principalmente ríos pequeños y no regulados son calibrados por el servicio nacional de alerta. Adicionalmente se realiza una predicción para las empresas de energía hidroeléctrica en unas 80 cuencas. Por otra parte, aplicaciones científicas u operacionales del modelo HBV han sido reportadas alrededor de 50 países de todo el mundo.

Estructura del modelo:

El modelo HBV puede ser descrito de mejor forma como un modelo conceptual semi-distribuido. Con el paso de los años sean realizado pequeños cambios en la estructura del modelo base desde que fue creado. Los datos de entrada se han mantenido lo más sencillamente posible, normalmente sólo los valores medios diarios de precipitación se utilizan en este algoritmo. A pesar de su simplicidad, el rendimiento de la simulación es admirable, y el uso original para la predicción hidrológica se ha extendido a aplicaciones tales como llenar los vacíos en series de datos en el tiempo, así como la simulación de caudales en los ríos sin datos de aforo.

El algoritmo esta diseñado y adaptado a partir de la estructura básica del software IHMS/HBV creado por los suecos, este algoritmo se adapto a las condiciones climáticas de la región centro americana, motivo por el cual no se incluyó la temperatura dentro de los datos de entrada.

Rutina del suelo

La cuantificación de humedad del suelo es la parte principal de la rutina que controla la formación de escorrentía. Esta rutina es basada en los tres parámetros β, lp, y fc, como se observa en la figura.

β Controla la contribución de la función respuesta (DQ/DP) o el incremento de la humedad del suelo en el almacenamiento (1- DQ/DP) por cada milímetro de lluvia. La relación DQ/DP también es llamada coeficiente de escorrentía, y DQ es llamada precipitación efectiva. DQ/DP puede ser expresada como R/IN.

lp es el valor de la humedad del suelo sobre el cual la evapotranspiración alcanza su valor potencial, y fc es el valor máximo que alcanza la humedad del suelo en el almacenamiento (en mm) del modelo, también llamado capacidad de campo. El parámetro lp esta dado como una fracción de fc.

El efecto de la rutina del suelo en la contribución de la escorrentía por la lluvia es pequeña cuando el suelo está seco (valores bajos de humedad del suelo) y grande cuando el suelo esta húmedo. La evapotranspiración potencial disminuye a medida que el suelo se seca.

Valores medios de series largas son usadas para estimar la evaporación potencial en una determinada época del año. Por lo que asume que la variación interanual de la evapotranspiración real depende más de las condiciones de humedad del suelo que de la variación interanual de la evaporación potencial.

Rutina de respuesta

La rutina de generación de escorrentía es la función de respuesta que transforma el exceso de agua de la humedad del suelo en escorrentía. También incluye el efecto de la precipitación directa y evapotranspiración en una parte que representa los lagos, ríos y otras zonas húmedas. La función consta de un depósito superior no lineal y uno inferior lineal. Estos son los que influyen en que el escurrimiento sea rápido y lento dentro de los componentes del hidrograma.

El rendimiento de la humedad del suelo, por ejemplo, la precipitación efectiva, se añadirá al almacenamiento del depósito superior. Mientras haya agua en el depósito superior, el agua se percolara al depósito inferior según el parámetro de perc. La percolación no es suficiente para vaciar el depósito superior, y la descarga generada tendrá una contribución directa desde el depósito superior, que representa el drenaje a través de más canales superficiales. El depósito inferior, por el contrario, representa el almacenamiento de las aguas subterráneos de la cuenca que contribuyen al flujo base.

El caudal de salida del depósito superior se describe por una función que corresponde a un coeficiente de recesión en continuo aumento.

Q_o = K_o * UZ^{(1 + α)}

Donde:

Q_o = Caudal de salida del depósito superior

K_o = Coeficiente de recesión

UZ = Contenido del depósito superior

α = Coeficiente de no linealidad

El parámetro alfa es un coeficiente de no-linealidad, por lo general es de orden 1.

Una estimación del coeficiente de recesión se puede realizar con el hidrograma de los datos observados y utilizado como primera aproximación del valor k.

El caudal de salida del depósito inferior se describe por:

Q4 = K₄ * LZ

Donde:

Q4 = Caudal de salida del depósito inferior

K₄ = Coeficiente de recesión del depósito inferior

LZ = Contenido del depósito inferior

Función de transformación

La escorrentía generada a partir de la rutina de respuesta se convierte a través de una función de transformación, con el fin de obtener la forma correcta del hidrograma en la salida de la subcuenca. La función de transformación es una técnica de filtro simple, con una distribución triangular de los pesos.

Bajar el archivo de excel