INTRODUCCIÓN

Para predecir posibles crecientes, se analizan y aplican métodos hidrológicos que transforman lluvias de diseño en las predicciones buscadas, teniendo en cuenta las condiciones físicas actuales o futuras de las áreas o cuencas objeto de estudio. El proceso comienza con la obtención de las curvas Intensidad-Duración-Frecuencia, (IDF), las cuales representan las características relevantes de las tormentas que ocurren en la zona (^{Aparicio, 1997}; ^{Rodríguez, 2018}). Al respecto, ^{Carballo et al. (2013}); ^{Quispe (2018}), generaron modelos matemáticos IDF para estimar lluvias de diseño, determinando la intensidad de lluvias extremas para varios rangos de duración; ^{Pizarro Tapia et al., (2013}) establecieron una plataforma interfaz georeferenciada para procesar información en forma remota con modelos matemáticos para obtener curvas IDF, los datos se tomaron tanto de estaciones pluviográficas (las intensidades máximas se calculan para duraciones desde 15 minutos a 24 horas), como de estaciones pluviométricas que sólo contaban con datos cada 24 horas; ^{Mamun bin Salleh et al. (2018}) con datos de precipitación máxima anual, mediante análisis estadístico y técnicas de ajuste de gráficos logarítmicos desarrollaron una correlación entre las precipitaciones de corta duración y los valores de precipitación diaria, analizaron la información de la precipitación para tiempos de 15, 30 y 45 minutos y cuyos porcentajes de lluvia diaria son de 32.4 %, 47.1 % y 57.4 % respectivamente; ^{Manzano et al., (2015}) mediante el análisis de los registros diarios de lluvias, diseñaron para México un mapa que muestra la regionalización de la relación entre las intensidades máximas de precipitación en intervalos de 1 y 24 horas (parámetro K) y observan una alta variabilidad espacial de este parámetro en todo el país; y ^{Muñoz y Zamudio (2018}) aplican el método de las isolíneas para regionalizar curvas IDF en el departamento de Boyacá, según la relación entre los parámetros que describen una ecuación de estas curvas y su ubicación geoespacial.

Frecuentemente, para obtener estos modelos, se comprueba la homogeneidad, la no estacionariedad de los registros de precipitación correspondientes a diferentes periodos e intervalos de tiempo, se determinan y comparan las curvas IDF utilizando distribuciones de probabilidad y análisis de eventos extremos (López et al., 2018; ^{Acosta y Sierra, 2013}; Gutiérrez et al., 2017). Cuando no se dispone de registros suficientes de intensidad de lluvia, se aplican técnicas como: relacionar información por grupos (los máximos diarios anuales, los máximos mensuales y los satelitales de precipitaciones tropicales) (Ayman y ^{Awadallah 2013}); establecer relaciones generales entre la intensidad, duración y frecuencia debido a que las propiedades convectivas de las precipitaciones en periodos cortos son similares en diferentes regiones hidrológicas (^{Kothyari y Ramchandra 1992}). Los modelos IDF se aplican para: identificar las zonas más propensas a inundaciones según períodos de retorno de lluvias extremas, sintetizando los períodos críticos (^{Correa, 2016}; ^{Rodríguez, 2018}); diseñar la infraestructura necesaria en los sistemas de drenaje y protección contra crecientes (^{Campos, 2010}); analizar los efectos del cambio climático en el comportamiento futuro de las precipitaciones (^{Rodríguez et al., 2013}, ^{Arnbjerg, 2012}); solucionar algunos problemas de ingeniería hidráulica, relacionando los períodos de retorno de una tormenta, con la magnitud de las precipitaciones y la distribución tanto espacial como temporal de las mismas (Berríos, 2008). En esta dirección, para analizar la precipitación en un sitio particular, existen básicamente dos métodos con los que se pueden determinar las relaciones entre las variables intensidad I, duración d y período de retorno T. El primero es denominado intensidad-período de retorno, que relaciona estas dos variables por separado mediante alguna de las funciones de probabilidad usadas en hidrología y el segundo, que se utilizó en esta investigación, relaciona simultáneamente las tres variables en una familia de curvas según la ecuación (1), donde las constantes k, m, n y c se calculan mediante un análisis de correlación lineal múltiple entre las variables (^{Aparicio, 1997}; ^{Chen, 1983}).

Pese a la relevancia y a la cantidad de investigaciones en el tema, el problema hidrológico de describir el comportamiento de la precipitación del rio chulo que pasa por la ciudad de Tunja es único, por las particularidades climatológicas de la región. Al respecto, el clima en esta ciudad se clasifica como oceánico, húmedo, cálido y templado, CFB (McKnight Hess, 2000); con promedio de precipitación anual de 917 mm, el mes más seco es enero y el de mayor precipitación es octubre con promedios de 23 mm y 122 mm de lluvia respectivamente; la temperatura media de 12.8 grados centígrados, caracterizada por pequeñas o moderadas oscilaciones térmicas diurnas y anuales; marzo es el mes más cálido del año y junio el mes más frio con temperaturas de 13.7 y 11.9 grados centígrados en promedio; además con una humedad relativa elevada.

Por tanto, el objetivo de la investigación corresponde a generar un modelo para estimar lluvias extremas de diseño en la cuenca inicial de este rio. Los resultados son los valores de los parámetros de los modelos que relacionan la intensidad, duración y frecuencia de lluvia, correspondientes a cuatro estaciones climatológicas que inciden en esta cuenca. El estudio contribuye a caracterizar la cuenca del rio Chicamocha que es un factor de desarrollo para el sector agrícola, ganadero e industrial de la región central del departamento de Boyacá.

METODOLOGIA

Para el diseño, desarrollo y simulación del modelo, se realizaron las siguientes etapas: 1) Recolectar y analizar la información de las precipitaciones diarias, medidas por el pluviógrafo y registradas en la estación climatológica de la Facultad de Agronomía de la Universidad Pedagógica y Tecnológica de Colombia ECFA, para calcular los valores máximos mensuales de precipitación de cada año en el período de estudio, 1967-2016; 2) Calcular la intensidad en mm/h, a partir de los valores de cada serie dividida por su duración (en horas); 3) Evaluar funciones de distribución de probabilidad en los valores de intensidad de precipitación; 4) Determinar qué distribuciones de probabilidad describen mejor la serie histórica de los datos, aplicando pruebas de bondad de ajuste; 5) Obtener el modelo matemático que comprende: calcular las precipitaciones diarias máximas probables para diferentes periodos de retorno, con la función de distribución de probabilidad adoptada; establecer las precipitaciones máximas para tiempos de duración de lluvia de 1, 2, 3, 4, 5, 6, 8, 12, 18 y 24 horas, como un porcentaje de la precipitación máxima probable para 24 horas, para cada periodo de retorno; obtener los parámetros de la ecuación de intensidad en función de la duración y tiempo de retorno aplicando regresión potencial; 6) Representar y analizar las curvas IDF para diferentes períodos de retorno, obtenidas al simular el modelo matemático; 7) Aplicar el procedimiento para encontrar los valores de los parámetros del modelo en otras estaciones.

RESULTADOS

Realizadas las etapas descritas en la metodología, se evaluó el modelo matemático con los registros pluviométricos de las intensidades obtenidas de la ECFAT. Los parámetros de ajuste, se obtuvieron a partir de las precipitaciones máximas mensuales en 24 horas de los años de estudio. Con el modelo obtenido se representaron gráficamente las curvas IDF, que son de gran importancia porque permiten calcular la intensidad promedio para cierta probabilidad de excedencia y duración (^{Maidment, 1993}); además, las gráficas evidencian las características de las tormentas de la zona o región con respecto a las variables magnitud, duración y frecuencia (Campos,1998).

Funciones de probabilidad

Para seleccionar el modelo matemático que mejor ajuste los datos de precipitación medidos, se analizaron las siguientes funciones de distribución de probabilidad, frecuentemente usadas en hidrología (^{Abramowitz y Stegum, 1972}): Normal, Lognormal, Pearson III, y Gumbel.

Normal, representada por la ecuación (2) donde μ y σ son, respectivamente la media y la desviación estándar, (^{Kite, 1988}):

Lognormal, representada por la ecuación (3), con parámetros α y β dados respectivamente por los logaritmos naturales de la media y la desviación estándar de la variable aleatoria:

Tabla 1: Precipitaciones máximas mensuales de 1967 a 2016 (Datos tomados del IDEAM, 2018) 

Pearson III o Gamma de tres parámetros, ecuación (4) con parámetros α 1 , β 1 y σ 1 , (^{Kite, 1988}):

Gumbel, definida por la ecuación (5) con parámetros estadísticos α y β; σ y x son respectivamente la desviación y la media de los datos medidos (^{Gumbel, 2004}):

Pruebas de homogeneidad y ajuste

Prueba de homogeneidad Estándar (SNHT). Es una prueba paramétrica, que asume la hipótesis nula como la homogeneidad de la serie y la hipótesis alterna como la existencia de una fecha en la que hay un cambio en la media de la serie. La prueba estadística T(k), está definida por la ecuación (6) (^{Alexandersson, 1986}).

Donde z 1 y z 2 , definidas por la ecuación (7), son las medias de los primeros k años y de los últimos n-k años de la serie que son comparadas.

T(k) alcanza su valor máximo cuándo hay un punto de cambio ubicado en el año k. El estadístico de prueba T 0 se define en la ecuación (8).

Si T 0 es superior al valor crítico, la hipótesis nula se rechaza. Los valores críticos dependen del tamaño de la muestra cómo se relaciona en la Tabla 2.

Tabla 2: Valores críticos del estadístico T 0 de la prueba de homogeneidad normal estándar en función del tamaño de muestra (n). (Datos tomados de (^{Alexandersson y Moeberg, 1997})) 

Para la serie de precipitaciones máximas mensuales, n=50, el valor máximo de T(k) es T 0 = 2.01397 y se registró en el año 2014, k=48, lo cual sugiere un punto de cambio en la serie y según los valores de la Tabla 2, T 0 < T 0 crítico . Por tanto, el punto de cambio detectado en el año 2014 no es estadísticamente significativo y se deduce que la serie se considera homogénea.

Para evaluar el ajuste de las funciones de distribución de probabilidad a los datos, se sometió la información a las pruebas: Error cuadrático mínimo, Ji cuadrada χ 2 , Kolmogorov-Smirnov y la eficiencia de Nash-Sutcliffe.

Error cuadrático mínimo. Con los valores observados xo i y los estimados teóricamente xe i , para cada función de distribución, se procedió a calcular este error según la fórmula (9) los resultados están en la Tabla 3.

Tabla 3: Cálculos para obtener el error cuadrático entre las funciones de probabilidad 

Ji cuadrada χ 2 . Para esta prueba se dividieron los datos en k=7 intervalos de clase (Tabla 4) donde θ i es el número observado en el intervalo i.

Tabla 4: intervalos de clase 

En este caso, se calcula el valor del número esperado de eventos en el mismo intervalo ϵ i , según la ecuación (10), donde F S i y F L i es la función de distribución de probabilidad evaluada en el límite superior e inferior del intervalo i respectivamente y n es el número de datos.

El estadístico D c se calculó con la ecuación (11), obteniendo los siguientes valores para cada distribución: Normal, 131.48; Lognormal, 12.50; Pearson III, 50.76; y Gumbel, 6.67.

Calculados los parámetros D c , se determina el valor de una variable aleatoria con distribución χ 2 , para v=k-1-m grados de libertad y un nivel de significancia α, donde m es el número de parámetros estimados para cada distribución. Para aceptar una función de distribución de probabilidad, el parámetro D c correspondiente, debe cumplir la condición descrita por la expresión (12).

Kolmogorov-Smirnov. Compara el máximo valor absoluto de la diferencia D K entre la función de distribución observada F o ( x m ) y la estimada F( x m ), según la ecuación (13), con un valor crítico d que depende del número de datos y el nivel de significancia seleccionado (^{Benjamin y Cornell, 1970}; ^{Simard y L'Ecuyer, 2011}).

En esta prueba, si D K <d, se acepta la hipótesis nula, es decir, se acepta la función de distribución de probabilidad. Para esta prueba, la función de distribución de probabilidad observada se calcula mediante la ecuación (14), donde m es el orden del dato x m en la lista ordenada de mayor a menor y n es el número total de datos:

Al aplicar esta prueba, en la Tabla 5 se resumen los cálculos según la ecuación (10), con n=50 datos y un nivel de significancia del 5 % se obtiene un valor crítico d=0.1923, que al compararlo con los valores D K , números subrayados de la Tabla 5, se infiere que se aceptan la funciones de probabilidad Lognormal, Gumbel, Normal y se rechaza la función Pearson.

Tabla 5: Cálculo del estimador D K , para la prueba de Kolmogorov Smirnov 

Eficiencia de Nash-Sutcliffe. El criterio se define en la ecuación (15) y mide cuánto de la variabilidad de las observaciones es explicada por la función de probabilidad. Algunos valores sugeridos para la toma de decisiones son resumidos en la Tabla 6 (^{Nash y Sutcliffe, 1970}).

Tabla 6: Valores referenciales del Criterio de Nash-Sutcliffe (Datos tomados de Molnar, 2011 (Molnar, 2011)) 

Los valores del criterio de Nash- Sutcliffe (E) obtenidos para cada distribución, según la ecuación (15) , son los siguientes: Normal, 1.0099; Lognormal, 1.0060; Pearson III, 1.0053; y Gumbel, 1.0043.

Selección de la función de distribución

En la Tabla 7, se presenta el resumen de las calificaciones de las funciones de distribuciones de probabilidad según las pruebas aplicadas, asignando el valor de 1 a la “mejor” y 5 a la “peor”. De estos resultados se infiere que el orden en que mejor se ajustan las funciones de probabilidad a los datos corresponde a: 1) Gumbel, 2) Lognormal, 3) Pearson III y 4) Normal. Por esta razón, para calcular las curvas IDF, se toma la decisión de aplicar la distribución Gumbel.

Tabla 7: Pruebas de bondad del ajuste de los datos a las funciones de distribución de probabilidad 

Cálculo de las curvas IDF

La primera fase para hallar las curvas IDF, corresponde a calcular la precipitación máxima probable Pd, a partir de la precipitación máxima mensual por año medida por el pluviómetro, representada por x_i en la Tabla 1. Según estos datos, el promedio de la precipitación máxima mensual es x = 31.9380 mm, la desviación estándar es S=8.9889 mm y los valores de los coeficientes según la ecuación (5) corresponden a los relacionados en (16).

En la Tabla 8 se presentan las precipitaciones máximas probables para diferentes períodos de retorno Tr, la variable reducida está dada por la ecuación (17) la precipitación es XT’=µ+α·YT y la probabilidad de ocurrencia de la precipitación se obtiene al aplicar la función Gumbelecuación (5), con los parámetros XT’, µ, α. La columna 5, corresponde a la corrección del intervalo fijo XT=XT’·1.13, según el estudio de proporción de lluvia máxima verdadera de intervalo fijo (^{Weiss, 1964}), valor que se utilizó en la investigación.

Tabla 8: Cálculo de las precipitaciones diarias máximas probables para diferentes períodos de retorno 

Ecuación de intensidad

Las relaciones o cocientes de lluvia de 24 horas que se emplean para duraciones de varias horas, son los propuestos en ^{Campos (1978}) y que se relacionan en Tabla 9.

Tabla 9: Valores para las relaciones a la lluvia de duración 24 horas (Datos tomados de ^{Campos (1978})) 

El cálculo de la precipitación máxima para diferentes tiempos de duración, se presenta en la Tabla 10, que corresponde a los porcentajes de la precipitación máxima probable en 24 horas (establecidos en la Tabla 9), para cada período de retorno.

Tabla 10: Precipitaciones máximas para diferentes tiempos de lluvias 

Intensidades de lluvia a partir de Pd, según duración de precipitación y frecuencia de la misma

La intensidad es la tasa temporal de precipitación, la profundidad por unidad de tiempo (mm/h); comúnmente se utiliza la intensidad promedio, ecuación (18), donde P es la profundidad de lluvia en (mm) y Td es la duración en horas. La frecuencia se expresa en función del período de retorno T_r (ver Tabla 11).

Tabla 11: Cálculo de la intensidad de lluvia en mm/h en la ciudad de Tunja 

El modelo matemático descrito por la ecuación (19) relaciona simultáneamente: la intensidad máxima de precipitación i en mm/hr, la duración de la precipitación t en minutos, la frecuencia en función del período de retorno T y los parámetros K, m, n, que dependen de la zona de estudio.

Para calcular los valores de K, m, n, se aplica la regresión lineal múltiple a la ecuación i=d t -n , obtenida al reemplazar d=K T m en la ecuación (19), obteniendo los valores para d, n y el coeficiente de correlación R, para cada período de retorno, como se presenta en la Tabla 12 y en la Fig. 1.

Tabla 12: Regresión potencial por el método de distribución de Gumbel para un período de retorno de 2 años 

Fig. 1: Gráfica de la función intensidad de lluvia un período de retorno de 2 años 

Un proceso similar se realizó para los periodos de retorno T de 5, 10, 25, 50 y 100 años. Los valores d y n se resumen en la Tabla 13.

Tabla 13: Resumen de aplicación de regresión 

Con los valores de la Tabla 13, se realizó otra regresión potencial entre las columnas del período de retorno T y el término constante de regresión d, Tabla 14, Fig. 2.

Tabla 14: Regresión potencial, d=K T m , T período de retorno 

Fig. 2: Gráfica de la regresión d=KT^m 

Para d=K T m se encuentran los valores: K= 126.8174, m= 0.1508 y n= 0.6163. Por consiguiente, el modelo matemático que describe el comportamiento de las precipitaciones para la zona de influencia de la ECFA se representa por la ecuación (20).

En la Tabla 15, se presentan los valores de la intensidad I en mm/min para diferentes tiempos de duración t en minutos y correspondientes a períodos de retorno T en años, como resultado de aplicar el modelo descrito por la ecuación (14). En la Fig. 3, se presentan las curvas IDF, generadas a partir de este modelo, que se comportan como funciones decrecientes en la intensidad a medida que aumenta el tiempo de duración del evento de lluvia.

Tabla 15: Valores de Intensidad de precipitación según duración y frecuencia de precipitación 

Fig. 3 Curvas IDF, de la cuenca con datos de precipitación máxima en 24 horas 

Respecto a los resultados obtenidos en el modelo de la ecuación (20), por ejemplo para un período de retorno T de 2 años y una duración de 60 minutos, la intensidad de una lluvia es 11.29 mm/h. Esto se interpreta que en promedio un evento de lluvia con una intensidad de 11.29 mm/h mayor y de duración de 60 minutos se presenta en una vez cada 2 años.

La curva IDF para un período de retorno dado se interpreta como una curva masa de precipitación, por ejemplo para el caso particular de T=2 años, resulta de la ecuación (20) que i=140.7926 t -0.61638 y a partir de esta expresión se obtiene la altura de precipitación para este periodo, según la ecuación (21). La gráfica de la altura de precipitación en función del tiempo de duración se representa en la Fig. 4.

Los procedimientos anteriores se aplicaron a datos de las estaciones Chiquiza, el Encanto y Cómbita, que inciden en la cuenca del río chulo, que pasa por la ciudad de Tunja, obteniendo los valores de los parámetros K, m, n del modelo matemático según lo descrito por la ecuación (19), el coeficiente de determinación R² y se resume en la Tabla 16.

Fig. 4 Gráfica de altura de precipitación 

Tabla 16 Valores de los parámetros, K,m,n de otras estaciones climatológicas 

Tabla 16 (continuación)

DISCUSION

Respecto a las siguientes características de la precipitación desde el punto de vista hidrológico que inciden en el uso y control del agua en la cuenca del rio que pasa por la ciudad de Tunja: altura o intensidad, su distribución en el espacio y en el tiempo, y la frecuencia o probabilidad de ocurrencia; en esta investigación se relacionó la intensidad, la distribución en el tiempo y la frecuencia. Por consiguiente, se requiere continuar el estudio para establecer la relación entre la intensidad con la distribución en el espacio y la intensidad con la distribución en el espacio y el tiempo. Sobre la metodología para determinar las curvas IDF, en el estudio de (^{Acosta y Sierra, 2013}) aplicaron el método de intensidad-período de retorno, mientras que en esta investigación se relacionan simultáneamente las tres variables como se presenta en la Tabla 17. Estos dos estudios coinciden en que la función de distribución Gumbel, es la mejor que representa los datos y fue la utilizada para desarrollar estos trabajos.

Al respecto, el IDEAM para estimar las curvas IDF, los datos de intensidades son ajustados a la distribución de probabilidad Gumbel y el modelo adoptado es el propuesto en la ecuación (22), para estimar los parámetros utilizan el método L-momentos y son calculados para diferentes periodos de retorno y duraciones.

En particular, los parámetros obtenidos del modelo IDF para la estación UPTC código 2403513 en la ventana de observación de 1968 a 2010, son los relacionados en la Tabla 17, los resultados de las intensidades estimados en la investigación y los presentados por el IDEAM, son similares para los periodos de retorno y duración correspondientes con una diferencia entre estos dos de I IDEAM -I =137.8430.

Tabla 17 Valores de los parámetros del modelo (22) de la estación U.P.T.C (Datos tomados de curvas IDF estación UPTC IDEAM) 

Sin embargo, ajustando la corrección del intervalo fijo como XT=XT’·2.375, se encuentran los valores de los parámetros de la ecuación (23) mejorando una aproximación entre las intensidades de I IDEAM - I 2.375 =23.7628

Los resultados teóricos que genera el modelo (20) son aceptables y son los esperados, como se puede verificar al revisar los datos en las Tablas 11 y 15. Respecto a la comparación con otros modelos IDF, se observa por ejemplo: para un evento de lluvia de duración de una hora y para un periodo de retorno de 25 años, en (^{Correa, 2016}) la intensidad para la cuenca del noroeste de Guayaquil Ecuador es de 251.54 mm/h, en (^{Carballo, Paredes, y Guevara, 2013}) la intensidad para la región del estado Cojedes Venezuela es de 1.81 mm/h y para la ciudad de Tunja Colombia es 34.71mm/h. Por esta razón, se evidencian diferencias significativas de la precipitación, debido las múltiples variables que inciden en esta como: la posición geográfica, la altitud, la influencia del cambio climático, periodos de lluvia, crecimiento expansivo de las ciudades, entre otras.

Finalmente, las implicaciones de este estudio permitirán establecer modelos regionales de precipitación de la cuenca del rio Chicamocha que incide en la región central del Boyacá, ubicar zonas propensas a inundaciones, determinar la probabilidad de la intensidad de la precipitaciones para estimar la suficiencia de la necesidad de agua de cultivos de papa y pasto (^{Sepúlveda et al., 2015}), cebada, trigo, maíz, hortalizas, y de esta manera contribuir a la acertada toma de decisiones para la inversión agrícola y ganadera.

CONCLUSIONES

De los resultados mostrados, de su análisis y de su discusión, se pueden obtener las siguientes conclusiones, del modelo IDF de lluvias para la ciudad de Tunja: 1) la característica principal de la metodología empleada, radica en la posibilidad de predecir la ocurrencia de eventos de lluvia que son complejos e imposibles de estimarlos de manera confiable por métodos basados en las leyes de la mecánica o la física; 2) los procedimientos utilizados son de fácil aplicación, ya que basta con disponer de los datos de precipitación, aplicar funciones de probabilidad, regresión potencial y tratamientos matemáticos para obtener los parámetros del modelo; 3) los resultados obtenidos por este método son comparables dentro de intervalos estadísticamente aceptables con otros modelos de lluvia y no deben aceptarse dogmáticamente, asumiendo la disponibilidad de registros porque entre mayor número de registros los resultados serán más confiables; y 4) la metodología empleada para encontrar el modelo puede ser aplicada para analizar otras variables que inciden en el comportamiento del clima.