Ejemplo - Modelo hidrologico.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Antes de ver este ejemplo, deberias haber visto los Notebooks [Ejemplo - Onda de Sonido](Ejemplo%20-%20Onda%20de%20Sonido.ipynb) y [Ejemplo - Modelo deposito agua](Ejemplo%20-%20Modelo%20deposito%20agua.ipynb)\n",
    "\n",
    "# Modelo hidrológico\n",
    "En este Notebook vamos a crear un modelo hidrologico y vamos a utilizar widgets interactivos para intentar entender mejor como funciona el modelo, como influencian los parametros del modelo en el resultado y como podemos calibrar el modelo. \n",
    "\n",
    "Primero de todo, ¿que es un modelo hidrologico? Imagina que queremos predecir el hidrograma de un rio que recoje el agua de una cuenca cuando llueve. Para este proposito utilizamos un modelo hidrologico, que no es mas que un conjunto de ecuaciones que describen de una manera simplificada los procesos hidrologicos que ocurren en la cuenca. Estas ecuaciones incluyen diferentes parametros que describen algunas de las propiedades de la cuenca, por ejemplo las caracterisiticas del suelo.\n",
    "\n",
    "![diagrama cuenca](util/diagrama_cuenca.gif)\n",
    "\n",
    "En este ejemplo, vamos a utilizar un modelo simple (una adaptacion del modelo [HyMOD]:(https://doi.org/10.1002/9781118665671.ch14)) que tiene 5 parametros:\n",
    "\n",
    "- **Soil storage capacity** (mm): capacidad que tiene el suelo de retener agua de lluvia\n",
    "- **Evaporation ratio**: ratio de evaporacion o proporcion de lluvia que se evapora\n",
    "- **Infiltration ratio**: ratio de infiltracion proporcion del agua de lluvia efectiva (que no se evapora) que se infiltra en el suelo.\n",
    "- **Travel time - surface flow** (days): tiempo de concentracion del agua superficial o tiempo que tarda el agua superficial en llegar a la desembocadura de la cuenca.\n",
    "- **Travel time - underground flow** (days): tiempo de concentracion del agua subterranea o tiempo que tarda el agua subterranea en llegar a la desembocadura de la cuenca.\n",
    "\n",
    "En la imagen de abajo podemos ver como el modelo representa los procesos hidrologicos para asi obtener finalmente la prediccion del hidrograma del rio para los proximos meses. Como puedes ver el suelo es presentado como un deposito con una cierta capacidad que si es sobrepasada desborda generando flujo superficial y que tambien desaloja agua por su parte inferior generando agua subterranea.\n",
    "\n",
    "![diagrama modelo hidrologico](util/diagrama_modelo_hidrologico.gif)\n",
    "\n",
    "Vamos a crear el modelo y probarlo para comprender mejor como funciona.\n",
    "\n",
    "Lo primero como siempre es importar las librerias de funciones que vamos a necesitar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Periodo de simulacion\n",
    "Vamos a definir el periodo de simulacion, en este caso 120 dias (`T`), desde el 1 de Marzo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 120 # dias del periodo de simulacion\n",
    "dates = pd.date_range(start = '2022-03-01', periods = T)\n",
    "dates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Datos de entrada (inputs)\n",
    "Los inputs de nuestro modelo son los datos de lluvia (`prec`) y evapotranspiracion potencial (`etp`). Deben ser series temporales (diarias) y en formato array. Como solo vamos a hacer una prueba del modelo vamos a generar datos aleatorios con una distribucion uniforme pero podriamos haber usado datos reales, importandolos por ejemplo de un archivo Excel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prec = np.random.uniform(low = 0, high = 30, size = [T,1])\n",
    "etp = np.random.uniform(low = 0, high = 5, size = [T,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parametros del modelo\n",
    "Definimos los 5 paramtetros del modelo y le damos un valor inicial.\n",
    "\n",
    "![parametros_modelo_hidrologico](util/parametros_modelo_hidrologico.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sm_max = 10  # Maxima capacidad de retencion del del suelo (mm) [10 - 90]\n",
    "gamma  = 0.5 # Ratio de evaporacion () [0 - 1]\n",
    "alpha  = 0.5 # Ratio de infiltracion () [0 - 1]\n",
    "t_sup  = 1   # Tiempo concentracion del flujo superficial [0.8 - 2]\n",
    "t_sub  = 5   # Tiempo concentracion del flujo subterraneo [2 - 10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Declaracion de las variables del modelo\n",
    "Ahora declaramos (o inicalizamos) las variables del modelo. No hace falta que presteis mucha atencion a estas variables, simplemente quedaos con la idea de que van a formar parte de las ecuaciones del modelo que representan de manera simplificada los procesos hidrologicos de la cuenca. Como en cada paso (time-step) de la simulacion vamos a actualizar sus valores, necesitamos inicializarlas como un array de `T` ceros"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inf   = np.zeros((T,1))   # Tasa de infiltracion [mm/t]\n",
    "et    = np.zeros((T,1))   # Tasa de evapotranspiration [mm/t]\n",
    "sm    = np.zeros((T+1,1)) # Contenido de humdedad en el suelo [mm] (suponemos que el suelo esta seco inicialmente)\n",
    "sL    = np.zeros((T+1,1)) # Slow reservoir moisture [mm]\n",
    "sF    = np.zeros((T+1,1)) # Fast reservoir moisture [mm]\n",
    "Q_sub = np.zeros((T,1))   # Flujo subterraneo [mm/t]\n",
    "Q_sup = np.zeros((T,1))   # Flujo superficial [mm/t]\n",
    "Q_tot = np.zeros((T,1))   # Flujo total [mm/t]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulacion\n",
    "Para la simulacion vamos a escribir las ecuaciones del modelo dentro de un bucle (loop) `for` y asi realizar los calculos automaticamente para todos los dias (`T` = 120) desde el 1 de Marzo. Las ecuaciones son sencillas y son en su mayoria muy parecidas a la ecuacion de simulacion de un deposito (que vimos en el Notebook \"Ejemplo - Modelo deposito de agua\"). No tienes que entender las ecuaciones, de nuevo solo tener en mente que representan los procesos hidrologicos de la cuenca de una manera simplificada."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in range(T):\n",
    "\n",
    "    ##### Humedad del suelo #####          \n",
    "    sm_temp = max(min(sm[t] + prec[t] - inf[t], sm_max), 0)\n",
    "\n",
    "    ##### Infiltracion #####\n",
    "    inf[t] = inf[t] + max(sm[t] + prec[t] - inf[t] - sm_max, 0) + min(sm[t] + prec[t] - inf[t], 0)\n",
    "\n",
    "    ##### Evapotranspiracion #####\n",
    "    W = min(np.abs(sm[t]/sm_max)*gamma, 1) # Factor de correccion de la evapotranspiracion\n",
    "    et[t] = W * etp[t] # Calculo de la evapotranspiracion\n",
    "\n",
    "    ##### Humedad del suelo (t+1) #####\n",
    "    sm[t+1] = max(min(sm_temp - et[t], sm_max), 0) # Compute the water balance\n",
    "\n",
    "    ##### Flujo subterraneo ######\n",
    "    Q_sub[t] = 1/t_sub * sL[t]\n",
    "    sL[t+1] = sL[t] + (1-alpha)*inf[t] - Q_sub[t]\n",
    "\n",
    "    ##### Flujo superficial #####\n",
    "    sF[t+1] = sF[t] +  alpha*inf[t] - 1/t_sup*sF[t]\n",
    "    Q_sup[t] = 1/t_sup * sF[t]\n",
    "\n",
    "    ##### Flujo total #####\n",
    "    Q_tot[t] = Q_sup[t] + Q_sub[t]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizacion de los resultados\n",
    "Ahora podemos ver el resultado de la simulacion, en este caso representamos el hidrograma del flujo superficial, el flujo subterraneao y el hidrograma total (suma de los dos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,7))\n",
    "plt.plot(dates,Q_tot, label = 'total', color = 'blue')\n",
    "plt.plot(dates,Q_sup, linestyle = '--', color = 'orange',   label = 'superficial')\n",
    "plt.plot(dates,Q_sub, linestyle = ':',  color = 'green', label = 'subterraneo')\n",
    "plt.ylim(0,25) # para fijar los limites del eje y, en este caso 0 y 25\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si ahora quisieramos repetir la simulacion pero con distintos valores de los parametros para comprobar como cambia el resultado necesitariamos volver a la celda \"Parametros del modelo\" cambiar manualmente los valores y volver a ejecutar las celdas posteriores hasta generar la grafica del hidrograma... No suena muy eficiente. Vamos a hacer algo diferente.\n",
    "\n",
    "### Modelo hidrologico interactivo\n",
    "Para poder entender la influencia de los parametros del modelo en el resultado vamos a hacer nuestro modelo mas interactivo usando widgets, y mas concretamente deslizadores (sliders). Ahora vamos a crear una funcion `modelo_hidrologico` que va a contener nuestro modelo hidrologico y por medio de la funcion `interact` de la libreria Ipwidgets vamos a poder interactuar con los parametros del modelo mas facilmente."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(sm_max = (10, 90, 1), ratio_evap = (0.01, 0.99, 0.01), ratio_inf = (0.01, 0.99, 0.01), \n",
    "          t_sup = (0.8,2,0.1), t_sub = (2,10,0.1))\n",
    "\n",
    "def modelo_hidrologico(sm_max=10, ratio_evap=0.5, ratio_inf=0.5, t_sup=1, t_sub=5): # valores iniciales de los parametros\n",
    "    \"\"\"\n",
    "    Este modelo hidrologico es una adaptacion del modelo HyMOD. Tiene 5 parametros:\n",
    "\n",
    "    sm_max: Maxima capacidad de retencion del del suelo (mm) [10-90]\n",
    "    ratio_evap:  Ratio de evapotranspiracion () [0-1]\n",
    "    ratio_inf:  Ratio de infiltracion () [0-1]\n",
    "    t_sup:  Tiempo concentracion del flujo superficial [0.8 - 2]\n",
    "    t_sub:  Tiempo concentracion del flujo subterraneo [2 - 10]\n",
    "    \"\"\"\n",
    "    \n",
    "    #######################################################################\n",
    "    # Inicializar variables\n",
    "    #######################################################################\n",
    "    inf   = np.zeros((T,1))   # Lluvia infiltrada [mm/t]\n",
    "    sup   = np.zeros((T,1))   # Lluvia acumulada en superficie [mm/t]\n",
    "    et    = np.zeros((T,1))   # Tasa de evapotranspiration [mm/t]\n",
    "    sm    = np.zeros((T+1,1)) # Contendio de humdedad en el suelo [mm] (suponemos que el suelo esta seco inicialmente)\n",
    "    sL    = np.zeros((T+1,1)) # Slow reservoir moisture [mm]\n",
    "    sF    = np.zeros((T+1,1)) # Fast reservoir moisture [mm]\n",
    "    Q_sub = np.zeros((T,1))   # Flujo subterraneo [mm/t]\n",
    "    Q_sup = np.zeros((T,1))   # Flujo superficial [mm/t]\n",
    "    Q_tot = np.zeros((T,1))   # Flujo total [mm/t]\n",
    "\n",
    "    #######################################################################\n",
    "    # Simulacion\n",
    "    #######################################################################\n",
    "    for t in range(T):\n",
    "\n",
    "        ##### Humedad del suelo #####          \n",
    "        sm_temp = max(min(sm[t] + prec[t], sm_max), 0)\n",
    "\n",
    "        ##### Lluvia infiltrada #####\n",
    "        inf[t] = max(sm[t] + prec[t] - sm_max, 0) * ratio_inf\n",
    "        \n",
    "        ##### Lluvia acumulada en superficie #####\n",
    "        sup[t] = max(sm[t] + prec[t] - sm_max, 0) * (1- ratio_inf)\n",
    "\n",
    "        ##### Evapotranspiracion #####\n",
    "        W = min(np.abs(sm[t]/sm_max)*ratio_evap, 1) # Factor de correccion de la evapotranspiracion\n",
    "        et[t] = W * etp[t] # Calculo de la evapotranspiracion\n",
    "\n",
    "        ##### Humedad del suelo (t+1) #####\n",
    "        sm[t+1] = max(min(sm_temp - et[t], sm_max), 0)\n",
    "\n",
    "        ##### Flujo subterraneo ######\n",
    "        Q_sub[t] = 1/t_sub * sL[t]\n",
    "        sL[t+1] = sL[t] + inf[t] - Q_sub[t]\n",
    "\n",
    "        ##### Flujo superficial #####\n",
    "        Q_sup[t] = 1/t_sup * sF[t]\n",
    "        sF[t+1] = sF[t] +  sup[t] - Q_sup[t]\n",
    "\n",
    "        ##### Flujo total #####\n",
    "        Q_tot[t] = Q_sup[t] + Q_sub[t]\n",
    "\n",
    "    #######################################################################\n",
    "    # Visualizacion de los resultados\n",
    "    #######################################################################\n",
    "    # Dibujar la grafica\n",
    "    plt.figure(figsize=(20,7))\n",
    "    plt.plot(dates,Q_tot, label = 'hidrograma (flujo superficial + flujo subterraneo)', color = 'blue')\n",
    "    plt.plot(dates,Q_sup, linestyle = '--', color = 'orange',   label = 'flujo superficial')\n",
    "    plt.plot(dates,Q_sub, linestyle = ':',  color = 'green', label = 'flujo subterraneo')\n",
    "    plt.ylim(0,25) # para fijar los limites del eje y, en este caso 0 y 25\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calibracion del modelo hidrologico\n",
    "Un modelo hidrologico suele tener un gran número de parámetros. El usuario es quien decide el valor de estos parámetros para una aplicación particular en función de la información y datos que tengamos sobre los parámetros. Normalmente no vamos a tener mediciones directas de los valores de dichos parametros pero si es probable que tengamos datos climaticos historicos y del caudal del rio (hidrograma historico del rio). Con estos datos podemos inferir los valores de los parámetros encontrando los valores que hacen que el resultado del modelo se ajuste mejor al hidrograma historico del rio, esto se denomina **calibración** del modelo.\n",
    "\n",
    "La forma más sencilla de hacer esto es cambiando los valores de los parámetros de uno en uno y observando cómo esto cambia el resultado del modelo y como de bien se ajusta al hidrograma historico. Esto es similar a lo que hicimos en el Notebook [Ejemplo - Onda de Sonido](Ejemplo%20-%20Onda%20de%20Sonido.ipynb), solo que ahora no simulamos una onda de sonido sino el caudal de un rio y ahora no tenemos 3 parametros (amplitud, fase, frecuencia) sino 5 parametros.\n",
    "\n",
    "Intenta ahora, cambiando los valores de los 5 parametros, calibrar el modelo, es decir ajustar el hidrograma que simula el modelo al hidrograma historico."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(sm_max = (10, 90, 1), ratio_evap = (0.01, 0.99, 0.01), ratio_inf = (0.01, 0.99, 0.01), \n",
    "          t_sup = (0.8,2,0.1), t_sub = (2,10,0.1))\n",
    "\n",
    "def modelo_hidrologico(sm_max=10, ratio_evap=0.5, ratio_inf=0.5, t_sup=1, t_sub=5): # valores iniciales de los parametros\n",
    "    \"\"\"\n",
    "    Este modelo hidrologico es una adaptacion del modelo HyMOD. Tiene 5 parametros:\n",
    "\n",
    "    sm_max: Maxima capacidad de retencion del del suelo (mm) [10-90]\n",
    "    ratio_evap:  Ratio de evapotranspiracion () [0-1]\n",
    "    ratio_inf:  Ratio de infiltracion () [0-1]\n",
    "    t_sup:  Tiempo concentracion del flujo superficial [0.8 - 2]\n",
    "    t_sub:  Tiempo concentracion del flujo subterraneo [2 - 10]\n",
    "    \"\"\"\n",
    "    \n",
    "    #######################################################################\n",
    "    # Inicializar variables\n",
    "    #######################################################################\n",
    "    inf   = np.zeros((T,1))   # Lluvia infiltrada [mm/t]\n",
    "    sup   = np.zeros((T,1))   # Lluvia acumulada en superficie [mm/t]\n",
    "    et    = np.zeros((T,1))   # Tasa de evapotranspiration [mm/t]\n",
    "    sm    = np.zeros((T+1,1)) # Contendio de humdedad en el suelo [mm] (suponemos que el suelo esta seco inicialmente)\n",
    "    sL    = np.zeros((T+1,1)) # Slow reservoir moisture [mm]\n",
    "    sF    = np.zeros((T+1,1)) # Fast reservoir moisture [mm]\n",
    "    Q_sub = np.zeros((T,1))   # Flujo subterraneo [mm/t]\n",
    "    Q_sup = np.zeros((T,1))   # Flujo superficial [mm/t]\n",
    "    Q_tot = np.zeros((T,1))   # Flujo total [mm/t]\n",
    "\n",
    "    #######################################################################\n",
    "    # Simulacion\n",
    "    #######################################################################\n",
    "    for t in range(T):\n",
    "\n",
    "        ##### Humedad del suelo #####          \n",
    "        sm_temp = max(min(sm[t] + prec[t], sm_max), 0)\n",
    "\n",
    "        ##### Lluvia infiltrada #####\n",
    "        inf[t] = max(sm[t] + prec[t] - sm_max, 0) * ratio_inf\n",
    "        \n",
    "        ##### Lluvia acumulada en superficie #####\n",
    "        sup[t] = max(sm[t] + prec[t] - sm_max, 0) * (1- ratio_inf)\n",
    "\n",
    "        ##### Evapotranspiracion #####\n",
    "        W = min(np.abs(sm[t]/sm_max)*ratio_evap, 1) # Factor de correccion de la evapotranspiracion\n",
    "        et[t] = W * etp[t] # Calculo de la evapotranspiracion\n",
    "\n",
    "        ##### Humedad del suelo (t+1) #####\n",
    "        sm[t+1] = max(min(sm_temp - et[t], sm_max), 0)\n",
    "\n",
    "        ##### Flujo subterraneo ######\n",
    "        Q_sub[t] = 1/t_sub * sL[t]\n",
    "        sL[t+1] = sL[t] + inf[t] - Q_sub[t]\n",
    "\n",
    "        ##### Flujo superficial #####\n",
    "        Q_sup[t] = 1/t_sup * sF[t]\n",
    "        sF[t+1] = sF[t] +  sup[t] - Q_sup[t]\n",
    "\n",
    "        ##### Flujo total #####\n",
    "        Q_tot[t] = Q_sup[t] + Q_sub[t]\n",
    "\n",
    "    #######################################################################\n",
    "    # Visualizacion de los resultados\n",
    "    #######################################################################\n",
    "    # Cargar los datos historicos del caudal del rio (hidrograma historico)\n",
    "    from modulos.hidrograma_historico import hidrograma_historico\n",
    "    Q_tot_hist = hidrograma_historico(T,prec,etp)\n",
    "    # Dibujar la grafica\n",
    "    plt.figure(figsize=(20,7))\n",
    "    plt.plot(dates,Q_tot_hist, label = 'hidrograma historico', color = 'black')\n",
    "    plt.plot(dates,Q_tot, label = 'hidrograma simulado', color = 'blue')\n",
    "    plt.ylim(0,25) # para fijar los limites del eje y, en este caso 0 y 25\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Entender el concepto de calibracion de un modelo y para los que ya sabeis que es la calibracion, como ayuda a comprender el modelo en sí. Que significa cada parametro, que parametros influencian mas el resultado y comprobar que el model se comporta de una manera logica. Y con elementos interactivos como los widgets facilitamos esta tarea y por tanto la comprension del modelo (por ejemplo hidrologico) o de una ecuacion (por ejemplo, la de una onda de sonido)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}