Fig. 1.  Esquema de una cuenca hidrográfica como un libro abierto (20).

Una fuente de incertidumbre en este tipo de modelado es la exactitud del valor estimado de fricción. El flujo sobre los planos suele ser de profundidades tan pequeñas que es posible que prevalezca el flujo laminar. Por el contrario, es probable que el flujo en el canal sea turbulento. Sin embargo, en ciertos casos, especialmente cerca del extremo aguas abajo de los planos, el flujo puede estar en el régimen mixto laminar-turbulento. Para aplicaciones de rutina, los coeficientes de fricción generalmente se estiman de tal manera que tengan en cuenta el flujo laminar en los planos (6).

Las ecuaciones de onda cinemática son una opción común en la simulación de flujo en cuencas. Estas ecuaciones están bien establecidas y generalmente conducen a soluciones físicamente realistas, sin las complejidades numéricas y de otro tipo que suelen ser parte de las formulaciones de ondas dinámicas completas. Para derivar la ecuación de la onda cinemática, el enunciado usual de conservación de la masa en un volumen de control se combina con una forma simplificada de conservación de la cantidad de movimiento que sólo tiene en cuenta las fuerzas de fricción y del volumen (masa). El enunciado de conservación de la masa es el siguiente:

∂Q          ∂A            _____   + ______ = qL   ∂x           ∂t

(1)

en el cual Q = caudal, en m³/s; A = área, en m²; y q_L = flujo lateral, en m³/s/m.

El enunciado simplificado de la cantidad de movimiento es: S_o = S_f , en el cual S_o = pendiente de fondo, y S_f = pendiente de fricción. De hecho, éste es un enunciado del flujo uniforme, que conduce a una forma alternativa de expresar la cantidad de movimiento:

Q = α A β

(2)

en el cual α y β son un coeficiente y un exponente, respectivamente. Los valores de α y β contienen información sobre las características de fricción y sección transversal [por ejemplo, β = 3 para un canal ancho con fricción laminar; β = 5/3 para un canal ancho con fricción de Manning; β = 3/2 para un canal ancho con fricción de Chézy; y β = 4/3 para un canal triangular con fricción de Manning (11)].

La velocidad de propagación de las ondas cinemáticas se obtiene de la Ec. 2:

∂Q            Q            c = _____   = β ____ =  βv         ∂A             A

(3)

en la cual c = celeridad de la onda cinemática; y v = velocidad media del flujo.

La ecuación de onda cinemática se obtiene multiplicando las Ecs. 1 y 3 y haciendo uso de la regla de la cadena para obtener:

∂Q             ∂Q            _____  + c  ______ = cqL   ∂t               ∂x

(4)

La Ecuación 4 contiene términos de convección en el lado izquierdo y un término fuente en el lado derecho. Significativamente, la Ec. 4 no contiene un término de difusión. Esta ecuación puede resolverse por métodos analíticos o numéricos. Una clase típica de soluciones numéricas es aquélla en la cual los valores de ∂Q/∂t, ∂Q/∂x, c y q_L se expresan en términos de todos o algunos de los cuatro valores discretos adyacentes de Q, c y q_L en el espacio y tiempo (Fig. 2). Alternativamente, el valor de c podría mantenerse constante, lo cual da lugar a la solución lineal de la onda cinemática.

Fig. 2.  Discretización espacial y temporal de la ecuación de la onda cinemática.

Una solución numérica de la Ec. 4 siguiendo el procedimiento antes mencionado no está exento de inconvenientes. A menos que las derivadas espaciales y temporales están perfectamente centradas en la malla, se generará una cierta cantidad de difusión en la solución numérica (Fig. 2). Esta difusión numérica puede eliminarse centrando las derivadas espaciales y temporales, pero aún quedarí una pequeña cantidad de dispersión numérica (16). Esto último suele ser simplemente una molestia, pero en ciertos casos puede ser de tal magnitud que invalide la solución. Además, dado el hecho de que normalmente hay una pequeña cantidad de difusión en el problema físico, la eliminación completa de la difusión numérica es generalmente no recomendable.

Los esquemas totalmente descentrados siguen siendo populares entre los profesionales (6, 7, 10). Ponce et. al. (16) han demostrado que de los tres posibles esquemas totalmente descentrados (ver Fig. 3), los esquemas I y II son condicionalmente estables y convergentes a la solución analática en cuanto el número de Courant (definido como la relación entre la celeridad física c y la celeridad de la malla Δs/Δt) se aproxime a 1. El esquema I, estable para números de Courant menores o iguales a 1, equivale a "diferenciación a barlovento". El esquema II, estable para números de Courant mayores o iguales a 1, es exactamente la imagen opuesta del esquema I. El esquema III es incondicionalmente estable, pero no convergente para cualquier valor del número de Courant (16).

Fig. 3.  Tres esquemas numéricos totalmente descentrados para resolver esa ecuación de la onda cinemática (16).

Los esquemas I y II son complementarios, con versiones en amplio uso, por ejemplo, en la opción de enrutamiento de la onda cinemática utilizada en los modelos HEC-1 y su secuela, HEC-HMS (6). El esquema III (y varias versiones relacionadas), aunque no convergente, también se ha visto favorecido por su característica de estabilidad incondicional (7, 10). En cualquier aplicación práctica, es probable que los números de Courant locales varíen y sean diferentes de 1. En este caso, ninguno de estos tres esquemas resolvería la ecuación de la onda cinemática, y ciertamente no el esquema III, dada su demostrada falta de convergencia. En teoría, para números de Courant cercanos a 1, en cuanto la discretización se hace extremadamente pequeña, los esquemas I y II pueden aproximarse asintóticamente a la solución analítica de la ecuación de la onda cinemática. Con las potentes computadoras de hoy, esto es ciertamente una posibilidad. Sin embargo, la preocupación subsiste si la eliminación completa de la difusión numérica es la estrategia apropiada para este tipo de problema.

En resumen, el modelador hidrológico se enfrenta al siguiente dilema: si busca eliminar la difusión numérica a cualquier costo, la solución puede ser puramente cinemática y, por lo tanto, carecer de la propiedad de difusión natural inherente a los flujos físicamente realistas. Por otro lado, si el modelador usa un esquema totalmente descentrado como I, II o III, el resultado final serí una cantidad incontrolada de difusión numérica que sólo se asemejarí vagamente al proceso físico de difusión.

La solución de esta dificultad es "igualar" razonablemente las difusiones física y numérica, un procedimiento que se ha utilizado con éxito en el enrutamiento del flujo en canales y otras aplicaciones (9, 14, 18). La clara ventaja de este enfoque es que hace que la solución numérica sea esencialmente independiente del tamaño de la malla. En otras palabras, la convergencia se satisface con una malla mucho más gruesa que para los esquemas I o II. Dada la variación natural en los números de Courant que probablemente están presentes en aplicaciones prácticas, un método que converje para una amplia gama de tamaños de malla es ciertamente uno a tener muy en cuenta.

3.  FORMULACIÓN DE LA ONDA DIFUSIVA

La solidez del método de "igualación de difusividades" se ilustra aquí comparando los resultados de dos soluciones para la dinámica de la cuenca. La primera es una aplicación del esquema III a la ecuación de onda cinemática, Ec. 4, con el efecto de difusión numérica estando "descontrolado". El segundo es una aplicación del concepto de "igualación de difusividades". Por simplicidad, ambas soluciones se formulan en modo lineal, manteniendo constantes los parámetros de celeridad y difusividad. El enrutamiento de ondas lineales, comparado al de ondas no lineales, ha indicado que para cuencas pequeñas las características no lineales son generalmente pequeñas y pueden ignorarse al comparar la precisión numérica de esquemas alternativos (14).

La aplicación del esquema III a la Ec. 4 en un modo lineal conduce a la siguiente ecuación de enrutamiento:

n + 1                 n + 1            n               Q j + 1     =   C0Qj       +  C2Qj + 1  + C0QL

(5)

en la cual C₀ = C/(1 + C); C₂ = 1/(1 + C); y Q_L = flujo lateral promedio, C = número de Courant de la celda computacional, C = c (Δt /Δs), en la cual Δs es Δx o Δy.

Para los cálculos de flujo superficial, Q_L es el producto de la precipitación neta por el área de la celda. Para el enrutamiento en el canal, Q_L es la contribución lateral promedio en el tiempo, integrada a lo largo de la longitud de una celda. Con una estimación de la celeridad media de onda para el plano y el canal, junto con la especificación adecuada de las condiciones iniciales y de contorno, la Ec. 5 permite el cálculo de la respuesta de la cuenca para una precipitación neta dada.

La derivación del método de igualación de difusividades es paralela a la del método Muskingum-Cunge, con la adición de flujo lateral (3, 14, 15). Las ecuaciones de enrutamiento son las siguientes:

n + 1                n + 1          n            n               Q j + 1     =   C0Qj     +  C1Qj  +  C2Qj + 1   + C3QL

(6)

en la cual C₀ = (-1 + C + D)/(1 + C + D); C₁ = (1 + C - D)/(1 + C + D); C₂ = (1 - C + D)/(1 + C + D); C₃ = (2C)/(1 + C + D); y D = número de Reynolds de la malla, definido como la relación de difusividad física q/(2S_o) a difusividad numérica (cΔs)/2, como sigue (14):

q                  D  =  __________             So c Δs

(7)

en el cual q = flujo unitario (q = vd; d = profundidad de flujo); S_o = pendente de fondo; c = celeridad de la onda cinemática, y Δs = tamaño de la malla.

Debe tomarse en cuenta que hay grandes diferencias entre las Ecs. 5 y 6 con respecto a sus requisitos de datos. La ecuación 5 requiere sólo valores de c tanto para el plano como para el canal (o alternativamente: la velocidad media v y un valor adecuado de β, ver la Ec. 3). De otro modo, la Ec. 6 requiere valores de c, q y S_o tanto para el plano como para el canal (o alternativamente: β, v, d, y S_o ).

Dado que la difusividad física es v = q /(2S_o), se deduce que, a diferencia de la Ec. 6, la Ec. 5 no contiene ninguna información física al respecto de difusión. Sin embargo, no es sorprendente que la aplicación de la Ec. 5 muestra un comportamiento similar a la difusión. Esta última no está relacionada con el problema físico, es directamente atribuible al tamaño de la malla, y es una función de éste.

4.  FORMULACIÓN DE ONDA DIFUSIVA CON COMPONENTE DINÁMICO

El análisis anterior ha demostrado las ventajas de utilizar la onda difusiva en lugar de la onda cinemática en el modelado de la respuesta de una cuenca. Aparte del supuesto de linealidad, la onda de difusión se puede mejorar mediante la inclusión de la dinámica del flujo no permanente. De acuerdo con Dooge (4), la formulación de la onda difusiva con componente dinámico para el caso de fricción de Chézy en canales hidráulicamente anchos es la siguiente:

en la cual F = número de Froude, definido como F = v/(gd)^1/2, y g = aceleración de la gravedad. La ecuación 8 tiene la propiedad importante de que la difusividad física v es nula para F = 2, un umbral de estabilidad neutra en el que las ondas cinemáticas y dinámicas se propagan a la misma velocidad (para la fricción de Chézy en canales hidráulicamente anchos) (11, 12)

La ecuación 8 se puede generalizar para cualquier formulación de fricción y forma de sección transversal considerando que el número de Froude de estabilidad neutra es igual a 1/(β - 1). Por lo tanto, la difusividad física expresada en términos de β es la siguiente:

q                     v  =  _____ [1 - (β - 1)2F2]           2So

(9)

Dada la ecuación 9, se puede formular un modelo de onda difusiva con componente dinámico en base a la Ec. 6, con la Ec. 7 modificada para tener en cuenta el número de Froude y la dependencia de β de la difusividad física, Ec. 9. En efecto, para la onda de difusión con componente dinámico, la Ec. 7 se modifica a lo siguiente:

q                     D  =  ___________  [1 - (β - 1)2F2]            So c Δs

(10)

Con el número de Reynolds de la malla definido con la Ec. 10, la Ec. 6 representa una descripción lo más completa posible de la dinámica de las ondas, dentro del contexto de una formulación difusiva lineal. Sin embargo, la utilidad de la Ec. 10 en comparación con la Ec. 7 dependerá de si los fenómenos no permanentes que se están simulando se pueden interpretar correctamente como una onda dinámica (13). La experiencia indica que para una gran clase de problemas de interés práctico en hidrología, la escala de ondas bien puede estar en el rango de la onda difusiva, o aún cinemática.

5.  EJEMPLO ILUSTRATIVO

Las metodologías descritas en la sección anterior se comparan aquí aplicándolas a un problema típico de dinámica de cuencas. Siguiendo un uso generalizado, el esquema de la Ec. 5 se denomina "onda cinemática", aunque se reconoce su difusión numérica apreciable e incontrolada. El método de "igualdad de difusividades", Ecs. 6 y 7, se denomina en el presente artículo "onda de difusión". La onda de difusión con componente dinámico, Ecs. 6 y 10, aunque propiamente dentro de la familia de modelos de onda difusiva, se denomina aquí "onda dinámica" (véanse las Tablas 2 y 3).

El ejemplo consiste en una cuenca conceptualizada como un libro abierto con dos planos adyacentes a un canal (Fig. 1). Se asume que los planos son impermeables, con un exceso de precipitación igual a la precipitación total. El flujo de entrada a los planos es el exceso de precipitación y la escorrentía es por flujo laminar superficial en una dirección perpendicular al alineamiento del canal. El flujo de entrada al canal es por contribución lateral de los planos y el flujo de salida del canal es la respuesta de la cuenca.

Las dimensiones son las siguientes: longitud del plano (en la dirección del flujo en plano) = 120 pies (36.58 m); ancho del plano = 240 pies (73.15 m); longitud del canal = 240 pies (73.15 m). La precipitación es 3 pulg/hr (76.2 mm/hr), que, cuando se multiplica por el área de la cuenca, conduce a una concentración de flujo de 4 pies³/s (0.113 m³/s) en la salida (asumiendo sólo traslación y despreciando la difusión). Las pendientes de fondo son: pendiente de los planos S_op (en la dirección x) = 0.01; pendiente del canal S_oc = 0.01.

El flujo laminar en los planos y el flujo turbulento en el canal se calculan siguiendo las prácticas establecidas en la hidráulica (2). Las condiciones de flujo promedio en los planos son: descarga por unidad de ancho q_p= 0.004 pies²/s (0.00037 m²/s); velocidad media v_p = 0.5 pies/s (0.15 m/s); profundidad media de flujo d_p = 0.008 pies (0.0024 m); exponente de la curva de gasto (flujo laminar en secciones hidráulicamentes anchas) β_p = 3; y celeridad media de la onda c_p = 1.5 pies (0.46 m/s). Las condiciones de flujo promedio en el canal son la descarga por unidad de ancho q_c = 1 pie²/s (0.092 m²/s); velocidad media v_c = 3 pies (0.91 m/s); profundidad media de flujo d_c = 0.333 pies (0.1 m); exponente de la curva de gasto (flujo turbulento en canales triangulares) β_c = 4/3; y celeridad media de la onda c_c = 4 pies (1.22 m/s).

Estos cálculos conducen a un tiempo de concentración t_c = (120 pies/1.5 pies) + (240 pies/4 pies) = 140 s. Por lo tanto, se elige una duración de tormenta de 3 min para asegurar el logro del equilibrio en la salida de la cuenca. El volumen total de lluvia es de 0.15 pulg. (3.81 mm), o 720 pies³ (20.39 m³).

6.   RESULTADOS

Los tres métodos se probaron con respecto a su comportamiento con varios niveles de resolución de malla. Se eligieron cinco niveles de resolución, que varían de muy grueso a muy fino, como se muestra en la Tabla 1. Los resultados de las simulaciones se muestran en las Tablas 2 y 3, expresados ​​en términos de: (1) flujo máximo; y (2) tiempo pico.

La Figura 4 muestra los hidrogramas calculados para: (1) método cinemático (Esquema I); y (2) método difusivo (igualación de difusividades). Como era de esperar, los resultados del método cinemático variaron sustancialmente con la resolución de la malla. Con una resolución muy gruesa (A), el flujo máximo fue de 2.3507 pies³/s (0.0666 m³/s), mientras que para una resolución muy fina (E), el flujo máximo fue de 3.9641 pies³/s (0.1122 m³/s). Los resultados del método difusivo mostraron propiedades de convergencia mucho mejores que el método cinemático. Con una resolución muy gruesa (A), el flujo máximo fue de 3.9336 pies³/s (0.1113 m³/s), mientras que para una resolución muy fina (E), fue de 3.9971 pies³/s (0.1132 m³/s). Se concluye que el método difusivo proporciona independencia de malla para una amplia gama de niveles de resolución, mientras que esto no es el caso del método cinemático. La demostrada independencia de malla del método difusivo es decididamente una ventaja.

Tabla 2.  Caudal máximo calculado por tres métodos (en pies cúbicos por segundo).
Nivel de resolución (1)	Cinemático (2)	Difusivo (3)	Dinámico (4)
A	2.3507	3.9336	3.9390
B	3.0721	3.9797	3.9834
C	3.5921	3.9932	3.9954
D	3.8612	3.9964	3.9979
E	3.9641	3.9971	3.9985
Nota:  Consultar el texto para obtener una descripción del método. 1 pie3/s = 0.0283 m3/s.

Tabla 3.  Tiempo pico calculado por tres métodos (en segundos).
Nivel de resolución (1)	Cinemático (2)	Difusivo (3)	Dinámico (4)
A	240	180	180
B	210	180	180
C	195	180	180
D	187	180	180
E	184	180	180
Nota:  Consultar el texto para obtener una descripción del método.

Fig. 4.  Hidrogramas de flujo de salida calculados en función de la resolución de la malla (ver la explicación del parámetro de la curva en la Tabla 1): (a) método cinemático; (b) método difusivo.

Los resultados del método dinámico no se representaron en la Fig. 4 debido a su gran parecido con los del método difusivo (ver Tablas 2 y 3). Esta semejanza subraya la falta de una contribución dinámica significativa en el ejemplo mostrado. El criterio adimensional de Ponce et. al. (13) se puede utilizar para demostrar si la onda es realmente una onda difusiva. Según este criterio, la onda es difusiva si se satisface la siguiente desigualdad:

g    1/2                           TSo (_____)    ≥ 30             d

(11)

en la cual T = escala de tiempo representativa de la onda de avenida. En este caso, se puede suponer que la escala de tiempo es igual al doble de la duración de la precipitación T = 360 s. Para el flujo del canal, con S_o = 0.01 y d = 0.333 pies (0.1 m), la Ec. 11 conduce a un lado izquierdo que tiene un valor de 35.4, satisfaciendo así la desigualdad. Esto explica por qué la solución de onda dinámica no es muy diferente de la solución de onda difusiva para este problema en particular.

La conservación de la masa se controló mediante la integración de los hidrogramas de flujo de salida para todos los ensayos descritos en las Tablas 2 y 3. Todos los métodos y ensayos conservaron la masa (volumen) en un 100%. Además, se utilizó el método difusivo con niveles de resolución espacial que satisfacen los criterios de precisión correspondientes (17).

En suma, cabe señalar que si la discretización se hubiera hecho extremadamente pequeña, los caudales máximos calculados por los métodos difusivo y dinámico (ver Tablas 2 y 3) habráan convergido al valor de flujo máximo obtenido al concentrar el flujo en la salida, es decir 4 pies³/s (0.1132 m³/s). Esto se atribuye al hecho de que la translación (es decir, la convección) es el mecanismo físico predominante para el ejemplo dado, y la difusión juega un papel relativamente mucho menor. Por lo tanto, un método difusivo (y difusivo con componente dinámico) en el que la difusión es pequeña y controlada con parámetros físicos es ciertamente preferible a un método cinemático en el que la difusión numérica podría ser excesiva y no tener relación con la física del problema.

7.  RESUMEN Y CONCLUSIONES

Se formula un método de ondas difusiva para calcular la dinámica de las cuencas hidrográficas. El método tiene mejores propiedades de convergencia que los modelos de ondas cinemáticas que utilizan derivadas descentradas en sus formulaciones numéricas. Al descentrar las derivadas se introduce difusión numérica en cantidades incontroladas, y los resultados de la simulación dependen en gran medida de la resolución de malla elegida.

A diferencia de los esquemas descentrados, el esquema de onda difusiva se formula igualando las difusividades física y numérica. Esto da como resultado un control efectivo de la difusión numérica y conduce a simulaciones que son esencialmente independientes del tamaño de la malla.

El método de la onda difusiva se amplía al ámbito de las ondas dinámicas al incluir la dependencia del número de Froude de la difusividad física. La formulación resultante representa una descripción lo más completa posible de la dinámica de las ondas dentro del marco de la teoría de las ondas difusivas.

Los tres métodos: (1) cinemático, con difusión numérica incontrolada; (2) difusión, con difusividades físicas y numéricas iguales; y (3) dinámicos, con difusividad dependiente del número de Froude, se prueban para convergencia y precisión general. Se muestra que el método difusivo tiene propiedades de convergencia mucho mejores que el método cinemático. En particular, la demostrada independencia de la malla del método difusivo es una ventaja a tener en cuenta. Al menos para el problema de ejemplo probado aquí, se muestra que la solución del método dinámico no constituye una mejora sustancial sobre la del método difusivo.

APÉNDICE I.  REFERENCIAS