1.  INTRODUCCIÓN

Los métodos simplificados de enrutamiento de avenidas disfrutan de una posición de preeminencia dentro de la gama de métodos de enrutamiento disponibles. Estos métodos son sencillos de formular y utilizar, y bien pueden explicar la mayoría de los fenómenos de las ondas de inundación en lo referente a aplicaciones prácticas. Las versiones más refinadas simulan el modelo de onda difusiva y en esto radica su ventaja.

Un método simplificado que está recibiendo mayor atención es el método de Muskingum en el cual los parámetros de enrutamiento se calculan en función de las características hidráulicas del canal y las de la malla (1), en este documento denominado modelo difusivo de Muskingum-Cunge. Su comportamiento numérico, sin embargo, no se entiende aún completamente. Por ejemplo, no existe un criterio de resolución espacial generalmente aceptado para preservar la precisión. Una prueba de la contínua controversia es el número de contribuciones recientes, pero el asunto queda aún por esclarecer (2,7,10).

Una consecuencia directa de la resolución espacial insuficiente son los notorios caudales negativos que se han documentado en algunos casos (3,10). Mientras no se identifique la causa de esta anomalía, los métodos simplificados adolecerán de falta de credibilidad. Esto sólo puede obstaculizar su aceptación para aplicaciones prácticas.

El objetivo de este artículo es, por lo tanto: (1) revisar y examinar las contribuciones recientes sobre resolución espacial y temporal en el modelo difusivo de Muskingum-Cunge; y (2) desarrollar un criterio de resolución espacial para preservar la exactitud del método. Aquí se entiende por exactitud lo que se puede lograr dentro del marco del enrutamiento de ondas difusivas.

2.  CONTRIBUCIONES RECIENTES

Koussis (2) fue uno de los primeros en reconocer el efecto de una resolución espacial insuficiente en el hidrograma calculado por el modelo difusivo de Muskingum-Cunge y en vincular este efecto al tamaño de la malla. Propuso, "para obtener resultados razonables", la siguiente restricción:

2 X Δx Δt ≥ ________              c

(1)

en el cual Δt = intervalo de tiempo; Δx = intervalo de espacio; X = parámetro; y c = celeridad de la onda de inundación. Koussis no dio más detalles sobre cómo la Ec. 1 fue obtenida, y se atribuye aquí a su experiencia con el método.

Weinman y Laurenson (10), en su reciente revisión de métodos aproximados de enrutamiento de ondas de inundación, se enfocan en el modelo difusivo de Muskingum-Cunge y sus variantes. Señalan la posibilidad de salidas negativas poco realistas, pero sugieren que, para efectos prácticos, éstas son lo suficientemente pequeñas y de corta duración como para ser ignoradas si se satisface la siguiente desigualdad:

Δx          Tr _____ ≤ ______     c          2X

(2)

en el cual T_r = tiempo de crecida del hidrograma.

Es evidente la sorprendente similitud entre los criterios de Koussis y Weinman y Laurenson. Se puede implementar una base de comparación si T_r se expresa en términos de Δt. Asumiendo que en un caso dado T_r / Δt = n es un número entero mayor de 5, la Ec. 2 se reduce a lo siguiente:

Δx        n Δt _____ ≤ ______     c          2X

(3)

lo que indica que el criterio de Koussis es más conservador en la estimación de Δx que el de Weinmann y Laurenson.

Los autores (7) han informaron los resultados de los hallazgos preliminares que parecían indicar que si se debe preservar la precisión (es decir, evitar las salidas negativas), la siguiente desigualdad debe mantenerse:

CD ≥ ξ

(4)

en el cual C y D son números de Courant y de Reynolds de la malla, respectivamente, definidos de la siguiente manera:

Δt C = c  _____              Δx

(5)

qo D =  _________          So c Δx

(6)

en el cual q_o = caudal de referencia por unidad de ancho; y S_o = pendiente de equilibrio de la superficie de agua. Además Ponce y Theurer afirmaron que un valor de ξ = 0.25 proporcionaba una precisión satisfactoria en la mayoría de los casos.

Kundzewicz (3) aplicó el método de análisis de convolución al modelo de onda difusiva de Muskingum y concluyó que la Ec. 2 no es válida para el caso general. Además, afirmó que la búsqueda de un criterio general para suprimir las salidas negativas podría resultar infructuosa.

3.  INTERVALOS DE ESPACIO Y DE TIEMPO

Los criterios de precisión como las Ecs. 1, 2 o 4 proporcionan pautas para la selección de intervalos de espacio y tiempo si el objetivo es eliminar los flujos de salida negativos. De acuerdo a Ponce y Yevjevich (5), los intervalos de espacio y tiempo se utilizan junto con los propiedades hidráulicas y de la onda de inundación para calcular los parámetros de enrutamiento C y D, es decir, las Ecs. 5 y 6. A su vez, estos parámetros se utilizan en la ecuación de enrutamiento (Fig. 1):

n + 1                     n                   n + 1                   n Qj + 1 = C1 Qj + C2 Qj    + C3 Qj + 1

(7)

Fig. 1  Discretización espacial y temporal en el modelo de difusión de Muskingum.

en la cual C₁, C₂ and C₃ se definen como sigue:

1 + C - D C1 = ____________           1 + C + D

(8)

- 1 + C + D C2 = ____________           1 + C + D

(9)

1 - C + D C3 = ____________           1 + C + D

(10)

Intervalo de espacio. Dadas las Ecs. 1, 2 y 4, la elección del intervalo de espacio no está del todo clara. Un límite superior para Δx es ciertamente una necesidad en cualquier procedimiento numérico y el presente método no es una excepción. La Ecuación 1 parecería indicar que el límite superior debería ser el siguiente:

c Δt Δx ≤ ______            2X

(11)

Sin embargo, esta condición adolece de un inconveniente; no está definida para X = 0. En el enrutamiento de onda difusiva de Muskingum, X = 0 implica que Δx es fijo, y la Ec. 11 ya no podría ser aplicable.

El mismo tipo de crítica se aplica a la Ec. 2. Por otro lado, la Ec. 4 se expresa en términos de D en lugar de X, y la expresión para el límite superior de Δx basada en esta ecuación es:

qo Δt Δx ≤ (_______)1/2              So ξ

(12)

La cuestión de un límite inferior para Δx también ha provocado cierta confusión, principalmente debido a la existencia de una relación entre X y Δx:

qo Δx = 0.5 (1 - _________)                        So c Δx

(13)

Históricamente, el parámetro X se ha interpretado como un factor de ponderación, y esto llevó a Miller y Cunge (4), entre otros, a sugerir un límite inferior para X: X ≥ 0. Este concepto se repite en la revisión de Weinmann y Laurenson (10), que establece lo siguiente (p. 1529): Dado que X no puede disminuir de manera realista por debajo de cero, existe un límite para el grado permisible de subdivisión, lo que implica un límite inferior para Δx. Dada la Ec. 13, este límite inferior se puede calcular estableciendo X ≥ 0, a partir del cual:

qo Δx ≥ ______           S0 c

(14)

El lado derecho de la Ec. 14 se reconoce inmediatamente como la "longitud característica de un tramo" identificada por primera vez por Kalinin y Milyukov (4) en relación con el método de enrutamiento que lleva su nombre. Sin embargo, la experiencia con el método ha demostrado que, de hecho, no existe un límite inferior teórico para Δx. Ponce y Theurer (7), así como otros investigadores (3,9), han informado de la experiencia con valores de Δx que infringen la Ec. 14, y sin embargo, sin irrealidad física alguna adjunta a los resultados del cálculo.

La razón de este comportamiento es evidente cuando se reconoce la característica de igualación del momento en la Ec. 13. Siempre que se satisfaga esta ecuación, Δx puede elegirse arbitrariamente para infringir la Ec. 14, lo que implica la posibilidad de valores negativos de X. Los cálculos prácticos que utilizan tal condición ayudan a disipar el mito de que X es un factor de ponderación que debe restringirse en el rango 0.0 ≤ X ≤ 0.5.

Intervalo de tiempo. El criterio para seleccionar el intervalo de tiempo también ha generado alguna confusión. Un límite superior en Δt es claramente una necesidad si la solución debe permanecer dentro de un marco numérico razonable. A diferencia de la resolución espacial, la resolución temporal es un concepto más fácil de comprender en aplicaciones hidrológicas, ya que los hidrogramas se representan en el dominio temporal. La práctica establecida es establecer un límite superior para Δt tal que se cumpla la siguiente desigualdad:

Tr Δt ≤ ______            n

(15)

en la cual n = número entero que tiene un valor mínimo de 5. En la práctica, los hidrogramas están sesgados hacia la derecha; por lo tanto, la Ec. 15 garantiza un nivel satisfactorio de resolución temporal.

Además del requisito de resolución temporal, existe la duda de si el modelo difusivo de Muskingum está sujeto al criterio de estabilidad numérica de Courant que normalmente se aplica a los problemas hiperbólicos. Si éste es el caso, entonces se debe imponer una restricción adicional sobre Δt:

Δx Δt ≤ ______            c

(16)

que proviene de la siguiente relación:

C ≤ 1

(17)

Los cálculos prácticos, sin embargo, han demostrado que las infracción de las Ecs. 16 y 17 no dan lugar a inestabilidad numérica. Este comportamiento se explica recordando que el método es esencialmente una formulación en diferencias finitas de la ecuación de la onda cinemática. En teoría (6), esta característica proporciona estabilidad incondicional con respecto a la condición de Courant, invalidando así la Ec. 17.

Un límite inferior en Δt es principalmente una cuestión de recursos computacionales. La decisión sobre este asunto se deja al modelador.

Criterios de precisión. Del análisis anterior, se deduce que no existe una justificación teórica para un límite inferior para Δt o Δx. Además, aunque existe un límite superior para Δt dado por el requisito de resolución temporal, un límite superior similar para Δx no es claramente discernible. Dadas las Ecs. 5, 6 y 13, el criterio de Koussis se modifica aquí y se expresa en términos de C y D:

C + D ≥ 1

(18)

De esta forma, la sobredeterminación de la Ec. 11 se elimina. En efecto, para X = 0, D = 1 y la Ec. 18 se reduce a C ≥ 0. Dado que X = 0 implica que Δx es finito (Ec. 6), la condición C ≥ 0 se reduce a la forma trivial Δt ≥ 0.

El criterio de Weinmann y Laurenson, Ec. 3, conduce a:

nC + D ≥ 1

(19)

y se aplican comentarios similares a los de Koussis. El criterio de Ponce y Theurer continua siendo:

C D ≥ 0.25

(20)

La siguiente sección describe un programa de experimentos numéricos diseñado para arrojar luz adicional sobre la cuestión de los criterios de precisión.

Consistencia. Estrechamente relacionado con la precisión está el tema de la consistencia de un método de enrutamiento simplificado como el presente. En el contexto de este estudio, la precisión se refiere a evitar salidas negativas físicamente irreales; de otro modo, la consistencia se refiere a la capacidad del método para producir los mismos resultados independientemente del tamaño de la malla. Por ejemplo, si se obtiene un determinado hidrograma al final de un tramo dado usando diez (10) intervalos de espacio, un método consistente produciría el mismo hidrograma si en cambio sólo se hubieran usado dos (2) intervalos.

La consistencia se utilizará en el marco de un método de enrutamiento simplificado en el que el objetivo es optimizar la difusión numérica haciéndola coincidir con la difusión física. Este concepto no debe confundirse con la propiedad de consistencia de un esquema convencional de diferencias finitas, en el que el objetivo es minimizar la cantidad de error numérico para resolver adecuadamente la ecuación diferencial en consideración.

4.  EXPERIMENTOS NUMÉRICOS

Se llevó a cabo un programa de experimentos numéricos con el fin de evaluar la precisión (y por ende, la consistencia) del enrutamiento de difusión de Muskingum descrito por las Ecs. 7-10. Se utilizaron las siguientes consideraciones en el diseño del programa de ensayos:

1. Proporcionar una amplia variación en los números de Courant y Reynolds de la malla.
   
   
2. Proporcionar una amplia variación en el tamaño de la malla: intervalo de espacio Δx y el intervalo de tiempo Δt.

3. Proporcionar una muestra representativa de las propiedades del canal, específicamente el caudal por ancho unitario q_o, la pendiente de fondo S_o, y la celeridad de la onda de inundación, c.

Se probaron seis canales hipotéticos. Para cada canal, se especificaron dos condiciones de flujo: (1) Un flujo base bajo; y (2) un flujo base alto. Cada canal y condición de flujo se probó calculando el mismo hidrograma a través de doce (12) especificaciones de la malla, lo que hace un total de 144 corridas. Las doce (12) especificaciones de la malla representaron combinaciones de tres intervalos de tiempo y cuatro intervalos de espacio. Para simular hidrogramas naturales, se utilizó una distribución gamma como hidrograma de entrada. Esta función tiene la siguiente forma:

t                     tp - t Qt = Qb + (Qp - Qb) ( ___ )m exp ( ________ )                                     tp                   tg - tp

(21)

en la cual Q_t = flujo de entrada en el tiempo t; Q_b = flujo base; Q_p = flujo pico; t_p = tiempo desde el inicio de la escorrentía directa hasta el pico del hidrograma; t_g = tiempo desde el inicio de la escorrentía directa hasta el centro de gravedad del hidrograma; y m está dado por la siguiente ecuación:

tp m = _______          tg - tp

(22)

La Tabla 1 muestra las características del canal y del hidrograma para las doce (12) series de ensayos. Un resumen detallado de los resultados de los experimentos numéricos se da en la Ref. 8.

5.  ANÁLISIS DE RESULTADOS

Los resultados de los experimentos numéricos pueden ser interpretados a la luz de los dos parámetros de enrutamiento: número Courant, C, y el número de Reynolds de la malla, D. Para empezar, hay que señalar que Δx es proporcional tanto a C como a D. Por lo tanto, cuanto mayor sea el Δx, menores serán los valores de C y D.

La Figura 2 es un gráfico logarítmico de los valores de C y D asociados con los resultados del programa de ensayos. En esta figura también se muestran las líneas para C₁ = 0, C₂ = 0, y C₃ = 0. Un análisis de los resultados numéricos indica claramente que la precisión depende del valor de C₂. Los valores negativos de C₂ están asociados invariablemente con caídas en la parte ascendente del hidrograma de flujo de salida, lo que indica la tendencia de los flujos de salida a volverse "negativos" (es decir, valores menores que el flujo base). Esta tendencia aumenta a medida que aumenta el valor de C₂ en el rango negativo.

Fig. 2   Series de pruebas de números de Courant y Reynolds de la malla para determinar los criterios de precisión.

Los valores de C₂ en el rango 0.0-0.5 no muestran caídas en el hidrograma de salida. Sin embargo, dependiendo del problema, aún pueden permanecer pequeñas inexactitudes en los hidrogramas calculados.

Los valores de C₂ superiores a 0.5 están asociados con hidrogramas que no muestran ninguna de las anomalías anteriores. Además, una vez que se alcanza la precisión, digamos para C₂ ≥ 0.5, se garantiza las consistencia. En otras palabras, los hidrogramas calculados en el rango C₂ ≥ 0.5 muestran una excelente reproducibilidad del flujo máximo, el tiempo al pico y la forma del hidrograma, con la variación de Δx y Δt.

6.  CRITERIOS DE PRECISIÓN EN LA RESOLUCIÓN ESPACIAL

La evidencia experimental sugiere que C₂ y sólo C₂ es determinante de la precisión. De hecho, los valores de C₁ y C₃ fueron monitoreados a lo largo de la serie de ensayos, y se demostró que los valores negativos para cualquiera de ellos no afecta la precisión del hidrograma calculado. Este hallazgo apunta a la falacia de limitar los tres coeficientes de enrutamiento a valores no negativos para lograr una "ponderación adecuada", como sostienen Miller y Cunge (4), entre otros.

La condición C₂ ≥ 0 se examina más a fondo aquí. Esto conduce, con la Ec. 9, a C + D ≥ 1, que es el criterio de Koussis modificado desarrollado en la sección anterior.

Una condición más fuerte viene dada por C₂ ≥ 0.5, lo que conduce a

C + D ≥ 3

(23)

En general, se aplica una condición tal que C₂ ≥ ζ , en la cual ζ  es un número real positivo. Esto lleva a lo siguiente:

1 + ζ C + D ≥  ______  = κ                1 - ζ

(24)

Es instructivo examinar el significado físico de este criterio. Dadas las Ecs. 5 y 6, la Ec. 24 se puede expresar de la siguiente manera:

1 Δx ≤ ___ ( ΔxC + ΔxD )           κ

(25)

en la cual:

ΔxC = c Δt

(26)

es la longitud de Courant, definida como la longitud del canal recorrida por la onda de crecida en un intervalo de tiempo Δt; y

qo ΔxD = _____            So c

(27)

en la "longitud característica del tramo", definida en términos de la fricción y sección transversal del canal.

Además, se observa que:

qo = uo do

(28)

en la cual u_o = velocidad media del flujo; y d_o = profundidad del flujo; y

c = β uo

(29)

en la cual β = exponente de la relación caudal-área de flujo Q = α A^β, de la cual:

do ΔxD = _____            β So

(30)

Dadas las Ecs. 29 y 30, la Ec. 24 alternativamente se puede expresar de la siguiente manera:

1                       do Δx ≤ ___ ( β uo Δt + _____ )           κ                     β So

(31)

en la cual β, u_o, d_o, y S_o = propiedades del canal y del flujo; y κ = un parámetro de precisión determinado experimentalmente. Un valor de κ = 2 es suficiente para preservar la precisión para aplicaciones prácticas.

7.  RESUMEN Y CONCLUSIONES

Se lleva a cabo un análisis del enrutamiento de difusión de Muskingum. Se realizan experimentos numéricos para arrojar luz adicional sobre el criterio que establece un límite superior al intervalo espacial Δx con el fin de preservar la precisión. La experiencia computacional indica que para valores muy grandes de Δx, existe una tendencia a que ocurran flujos negativos, los cuales son físicamente irreales.

Los resultados de los experimentos numéricos indican que el valor de C₂ es determinante de la precisión. La condición C₂ ≥ ζ  es postulada para preservar la precisión (y consistencia) del método de enrutamiento, en el que ζ es un número real positivo. Esto conduce a un límite superior para Δx dado por la siguiente relación:

1 Δx ≤ ___ ( ΔxC + ΔxD )           κ

(25)

Se recomienda un valor de ζ = 0.33, el cual corresponde a κ = 2. Las cantidades Δx_C y Δx_D se definen en términos de características de canal y tamaño de la malla.

AGRADECIMIENTOS

La investigación en la cual se basa este documento fue financiada en parte por el Servicio de Conservación de Suelos (SCS), Departamento de Agricultura de los EE.UU., Lanham, Maryland. Se agradece la ayuda de George H. Comer, conservacionista de suelos.

APÉNDICE I.  BIBLIOGRAFÍA

1. Cunge, J. A. 1969. "On the Subject of a Flood Propagation Computation Method (Muskingum Method)," Journal of Hydraulic Research, Vol. 7, No. 2, 205-230.

2. Koussis, A. 1976. "An Approximate Dynamic Flood Routing Method," Proceedings of the International Symposium on Unsteady Flow in Open Channels, Paper L1, Newcastle-Upon-Tyne, England.

3. Kundzewicz, Z. W. 1980. Discussion of "Approximate Flood Routing Methods: A Review," by P. Erwin Weinmann and Eric M. Laurenson, Journal of the Hydraulics Division, ASCE, Vol. 106, No. HY12, Proc. Paper 15866, Dec., 2072-2075.

4. Miller, W. A., y J. A. Cunge. 1975. "Simplified Equations of unsteady Flow," Unsteady Flow in Open channels, K. Mahmood and V. Yevjevich, eds., Water Resources Publications, Fort Collins, Colo., 183-758.

5. Ponce, V. M., y V. Yevjevich. 1978. "The Muskingum-Cunge Method with Variable Parameters," Journal of the Hydraulics Division, ASCE, Vol. 104, No. HY12, Proc. Paper 14199, Dec., 1663-1667.

6. Ponce, V. M., Y. H. Chen, y D. B. Simons. 1979. "Unconditional Stability in Convection Computations," Journal of the Hydraulics Division, ASCE, Vol. 105, No. HY9, Proc. Paper 14807, Sept., 1079-1086.

7. Ponce, V. M., y F. D. Theurer. 1980. Discussion of "Approximate Flood Routing Methods: A Review," by P. Erwin Weinmann and Eric M. Laurenson, Journal of Hydraulic Division, ASCE, Vol. 106, No. HY11, Proc. Paper 15784, Nov., 1945-1947.

8. Ponce, V. M. 1981. "Development of an Algorithm for the Linearized Diffusion Method of Flood Routing," Research Reposrt, San Diego State University Civil Engineering Series, No. 81144, San Diego, Calif., Sept.

9. Weinmann, P. E. 1977. "Comparison of Flood Routing Methods for National Rivers," Report No. 2/1977, Department of Civil Engineering, Monash university, Clayton, Victoria, Australia.

10. Weinmann, P.E., y E. M. Laurenson.1979. "Approximate Flood Routing Methods: A Review," Journal of the Hydraulics Division , ASCE, Vol. 105, No. HY12. Proc. Paper 15057, Dec., 1521-1536.

11. Weinmann, P. E., y E. M. Laurenson. 1981. Closure to "Approximate Flood Routing Methods: A Review," Journal of the Hydraulics Division, ASCE, Vol. 107, No. HY9, Proc. Paper 16479, Sept., 1112-1114.