Artículo de revisión

Recibido: 12/01/2026

Aceptado: 20/01/2026

Aprobado: 10/02/2026

Matemática tensorial de los procesos hidrodinámicos en el karst: Teoría y aplicaciones en Cuba

Inversiones GAMMA, S.A, Cuba.

*E-mail: especialistaprincipal@gmail.com

Resumen

Contexto/Objetivo: La hidrodinámica de acuíferos kársticos es extremadamente compleja debido a su heterogeneidad y anisotropía. Los modelos tradicionales basados en escalares o vectores simples no capturan adecuadamente las direcciones preferenciales de flujo, lo que lleva a inexactitudes en la gestión de recursos hídricos. Este artículo revisa las aplicaciones de la matemática tensorial para modelar el transporte de masa, momento y energía en medios kársticos, con el objetivo de demostrar su superioridad sobre los enfoques tradicionales. Metodología: Se realizó una revisión integral de la literatura, incluyendo trabajos fundamentales y estudios de caso, con un enfoque en la formulación tensorial de propiedades como la conductividad hidráulica y la dispersividad. Se analizaron cuatro casos de estudio: intrusión marina, corrección de direcciones de flujo, transporte de contaminantes orgánicos y no acuosos, y simulación de espeleogénesis. Resultados clave: La matemática tensorial permite una representación realista de la anisotropía del karst, corrigiendo errores en la determinación de direcciones de flujo y en la predicción del transporte de contaminantes. Los casos de estudio muestran que los modelos tensoriales mejoran la precisión en la gestión de la intrusión marina y en la simulación del desarrollo de cavernas. Sin embargo, se identifican limitaciones debido al desconocimiento de ciertos coeficientes fenomenológicos en los procesos de termodinámica de no equilibrio. Conclusiones/Relevancia: El enfoque tensorial es esencial para una modelación hidrodinámica precisa en el karst. Su adopción puede mejorar significativamente la gestión de recursos hídricos, la protección contra la contaminación y la evaluación de riesgos geotécnicos. Se recomienda investigación futura para determinar los coeficientes fenomenológicos faltantes y para integrar completamente los procesos de no equilibrio en los modelos.

Palabras clave: Matemática tensorial, Hidrogeología kárstica, Anisotropía hidráulica, Transporte de contaminantes, Intrusión marina, Espeleogénesis, Modelado numérico.

Tensor mathematics of the hydrodynamic processes in karst: Theory and applications in Cuba

Abstract

Background/Purpose: Karst aquifers exhibit extreme hydrodynamic complexity due to pronounced heterogeneity and anisotropy, which cannot be accurately described by traditional scalar or vector-based models. This review synthesizes the critical application of tensor mathematics to model mass, momentum, and energy transport in these anisotropic media. Its purpose is to demonstrate the transformative role of tensor calculus in advancing the physical realism of karst hydrogeological models and to highlight the inaccuracies inherent in conventional isotropic approaches for key management issues. Methods: A comprehensive analysis is conducted through a review of foundational tensor theory and its specific application to karst hydrology. The methodology examines four detailed case studies: (1) seawater intrusion dynamics, (2) subsurface flow direction correction, (3) transport of organic pollutants and non-aqueous phase liquids (NAPLs), and (4) speleogenesis simulation. The formulation and solution of tensor-based governing equations for anisotropic hydraulic conductivity and dispersion are central to this analysis. Key Findings: Tensor mathematics is confirmed as an essential framework for capturing the directional dependence of hydraulic properties in karst. The case studies reveal that tensor-based models correct significant errors in flow path identification, provide a more accurate prediction of contaminant plume migration (including NAPLs), and offer a thermodynamic foundation for simulating cavern development. A major constraint identified is the limited knowledge of key phenomenological coefficients within the non-equilibrium thermodynamic processes governing karstification. Conclusions/Implications: The adoption of tensor formalism is imperative for progressing from qualitative descriptions to robust quantitative predictions in karst hydrogeology. This approach directly enhances the accuracy of water resources management, contamination remediation strategies, and geohazard assessment. Future research must prioritize the empirical determination of unknown phenomenological coefficients to fully integrate non-equilibrium thermodynamics into predictive models, thereby closing a critical knowledge gap in the field.

Keywords: Contaminant transport, Hydraulic anisotropy, Karst hydrogeology, Numeric modeling, Seawater intrusión, Speleogenesis, Tensor mathematics.

1. Introducción

La matemática tensorial es una herramienta fundamental para modelar la compleja hidrodinámica de los acuíferos cársicos. La aplicabilidad reside en su capacidad para describir propiedades y procesos que varían en dirección y magnitud dentro de un medio heterogéneo y anisotrópico por excelencia, como lo es el karst. En hidrogeología cársica, la matemática tensorial se utiliza principalmente para modelar el comportamiento anisotrópico y heterogéneo del flujo de agua ya que, a diferencia de los medios porosos homogéneos, el karst presenta direcciones preferenciales de flujo (fracturas, conductos, estratificación) que no pueden describirse con escalares o vectores simples.

La formulación tensorial proporciona un marco matemático necesario y suficiente para modelar la complejidad hidrodinámica del karst, permitiendo pasar de descripciones cualitativas a predicciones cuantitativas robustas para la gestión y protección de sus valiosos recursos hídricos. El enfoque tensorial permite una representación física realista de la anisotropía estructural y la conservación de propiedades de simetría en ecuaciones. Especialmente importante es su formulación invariante bajo rotación del sistema coordenado y, sobre todo, su capacidad para acoplar múltiples procesos (flujo, transporte, calor, mecánica).

Así, en los últimos años se han obtenido excelentes resultados en la modelación del flujo preferencial en los acuíferos cársicos, desde los modelos más simples de doble porosidad/permeabilidad hasta aquellos de múltiple porosidad que incorporan espacios múltiples jerarquizados y la integración de heterogeneidades aplicando la promediación estadística de longitudes, volúmenes o áreas elementales representativas –incluso para las zonas de flujo y no flujo separadamente-, entre otras. Estos recursos analíticos han sustentado importantes soluciones ingenieras para la gestión de este tipo de acuíferos, como son las relacionadas con la interpretación correcta de la dirección y magnitud del flujo, ensayos de caudal o pruebas con trazadores, la identificación de zonas de alta recarga preferencial alineadas con sistemas de dolinas/fracturas y la aplicación de modelos de transporte reactivo para simular e interpretar la dispersión anisotrópica de contaminantes, incluyendo su propia capacidad de autodepuración.

Las principales aplicaciones de este recurso matemático, entonces, son las siguientes:

- Caracterización de la permeabilidad y anisotropía

- Generalización de la ley de Darcy

- Modelación del flujo en medios dobles o múltiples

- Simulación de transporte de contaminantes

- Análisis de esfuerzos y deformación en la roca

- Interpretación de datos geofísicos y de trazadores

- Modelos numéricos avanzados (e.g., EF, VDF)

Las bases teóricas de este tipo de análisis han sido desarrolladas, entre otros, por Schoeller (1962); Bear (1972); Freeze y Cherry (1979); Bear y Bachmat (1990); Bear y Buchlin, 1991; Bear et al., (1968) y Renard y de Marsily (1997). Pero en su aplicación al karst destacan los trabajos de Kiraly (1975); Long et al. (1982); Sauter y Liedl, (2000); Kovacs (2003); Reinmann y Hill (2009); de Röij et al. (2013): Hartmann et al. (2014); Chen et al., (2017), Zhang et al. (2021). Estos últimos cinco autores han desarrollado contribuciones específicas a la modelación numérica de los problemas hidrodinámicos del karst bajo principios de análisis tensorial de algunas propiedades físicas de sistema cársico.

En Cuba los métodos se han aplicado en todas las líneas de actuación recién mencionadas. Así, en temas como la caracterización de la conductividad hidráulica y la anisotropía, la generalización de la ley de Darcy y la modelación del flujo en medios dobles o múltiples (Molerio, 1982, 1985a, 1986a) han constituido la base del Modelo Termodinámico del Desarrollo de Karst (MTDC) desarrollado por el autor (Molerio, 1985b, 1985c, 1988, 1989a, 1989b, 1989c, 1989d, 1990a, 1992a, 1992b, 1992c, 1992d, 2003a, 2003b, 2013a, 2024a; March y Molerio, 1987); modelos de transporte reactivo y transporte de contaminantes y de masas en general, incluyendo trazadores (Molerio et al., 1996; Molerio, 2005, 2021, 2022a, 2025a, 2025b) y en la prospección hidrogeológica en medios anisótropos (Molerio et al., 1998; Molerio, 2005, 2023a; Núñez y Molerio, 2022). Especialmente importantes han sido sus aplicaciones para reajustar las direcciones de flujo subterráneo (Molerio, Alonso y Sardiñas, en prensa) y la extensión de estabilidad de dolinas y cavernas con aplicaciones generalizadas en Ecuador y Filipinas (Molerio et al., 1990; Molerio, 1990b, 1993, 2018a, 2018b, 2019a, 2019b, 2019c, 2020, 2022b, 2024b; Otero y Molerio, 2022). En algunos ejercicios docentes de posgrado en la Universidad de Pinar del Río, los métodos fueron introducidos hace muchos años, fundamentalmente para los estudios de transporte de masas (Antigüedad, et al., 1997; Molerio, 1997).

El tratamiento tensorial de la conductividad hidráulica en el karst es especialmente trascendente y es, de hecho, la base más rigurosa y robusta de prácticamente toda la hidrodinámica del karst. En medios porosos clásicos, la permeabilidad suele ser un escalar; pero en el karst, la permeabilidad es altamente direccional, dominada por fracturas, diaclasas y conductos y se representa como un tensor de segundo orden (3x3 en 3D). De este modo, captura cómo la conductividad hidráulica preferencial cambia según la dirección del flujo y, en el espacio euclidiano, K se describe como una elipse cuyo eje mayor, en función de la velocidad, se orienta a lo largo de una fractura principal (Fig. 1). Perpendicular a ella, la K menor se encuentra cuando se orienta en función de gradiente. Modelar K como un escalar subestimaría o sobreestimaría extraordinariamente los flujos reales; y más aún, es necesario considerar siempre que en karsts de desarrollo extremo -aquellos dominados por el flujo en conductos-, el concepto de medio continuo puede colapsar, requiriendo enfoques híbridos o discretos. Y, como se ha insistido en múltiples ocasiones, el Efecto de Escala es crítico, en tanto puede resultar en extremo distorsionador (Kiraly, 1975, 1978; Molerio, 1984, 1986b).

Fig. 1. Relación entre el gradiente (J) y la velocidad de filtración en medio anisotrópico (Kiraly, 1978).

La solución adecuada de la Ley de Darcy generalizada bajo K tensorial permite calcular el flujo de agua en cualquier dirección dentro de una red de fracturas interconectadas, lo que es esencial para entender la trayectoria y velocidad del agua subterránea en el karst. En la modelación del flujo en medios dobles o múltiples las ecuaciones de intercambio de masa entre estos dominios se acoplan usando parámetros de naturaleza tensorial cuando dependen de la dirección.

Asimismo, la dispersividad hidrodinámica (D), crucial para modelar cómo se propaga un contaminante, es inherentemente tensorial en medios anisotrópicos y, particularmente en un conducto cársico, la dispersión longitudinal (a lo largo del flujo) es muy diferente a la transversal. El tensor de dispersión permite modelar esta difusión direccional, vital para evaluar vulnerabilidad y elaborar planes de remediación.

La estabilidad de los terrenos cársicos es una amenaza de consideración, sobre todo en espacios urbanos y su tratamiento ha sido muy limitado en tanto es tomado, con frecuencia, como un caso fortuito. Solo en aquellos casos en que constituye un peligro inminente se han aplicado análisis más robustos como en los hundimientos sistemáticos de las carreteras de Alabama (Waltham et al, 2005). El estado de esfuerzo (σ) en la roca cársica es un tensor de segundo orden, en tanto afecta la apertura de fracturas y por tanto la conductividad hidráulica. Por ello, modelar cómo cambios en el campo de esfuerzos (por sismos, variaciones de carga hidráulica, diferencias locales en el campo tensional por la diferente competencia y fracturamiento de las rocas) alteran el tensor K, en tanto que, como acoplamiento hidromecánico afecta la variabilidad temporal de la conductividad hidráulica en el karst y, por ende, su papel en la redistribución del campo de esfuerzo/deformación.

Por otro lado, los datos de tomografía de resistividad eléctrica o de georradar se interpretan a menudo asumiendo un tensor de conductividad eléctrica análogo al de permeabilidad, lo que contribuye a inferir la estructura anisotrópica del subsuelo cársico (dirección de conductos preferenciales) a partir de mediciones indirectas, calibrando así los modelos hidrogeológicos. Experiencias recientes en el campo de dolinas de El Cayuelo han combinado técnicas de georradar, tomografía eléctrica, microgravimetría, migración de contaminantes y simulación del desarrollo de cavernas bajo un esquema conceptual que integra el sistema de flujo a una variación en el campo de esfuerzos/deformación y alivio de presiones en el sistema. Las no linealidades pueden tratarse oportunamente introduciendo un tensor de coeficiente de no-linealidad.

Las tendencias actuales se mueven hacia la aplicación de tensores de orden superior para la heterogeneidad no-Gaussiana, el desarrollo de métodos estocásticos tensoriales para incertidumbre en K(x), el aprendizaje automático para estimar parámetros tensoriales de datos indirectos y a computación en GPU para operaciones tensoriales masivas.

En este artículo se revisan con cierto detalle las prestaciones de la matemática tensorial en algunas de las aplicaciones relacionadas con la transformación y generalización de la estructura de la conductividad hidráulica anisotrópica en las ecuaciones de Darcy, Darcy-Weisbach, Richards y Forcheimer en el medio cársico en cuatro casos de estudio y se comentan las inexactitudes que tributa a la gestión de los recursos hídricos el enfoque tradicional y dominante de a) la intrusión marina; b) el transporte de contaminantes orgánicos, c) de líquidos de fase no acuosa y d) la simulación del desarrollo del cavernamiento. Se discuten varias restricciones derivadas del desconocimiento de algunos coeficientes fenomenológicos en los procesos de la termodinámica de no equilibrio que controla todo el proceso de carsificación y cavernamiento.

2. Análisis y discusión

2..1. Formulación tensorial de las ecuaciones hidrodinámicas básicas en el karst

El modelo conceptual del karst descrito por Molerio (1985a), lo define como un ente físico, un medio acuífero peculiar cuyo espacio está caracterizado por las siguientes propiedades:

- Es un sistema termodinámico abierto, es decir, en interacción con el medio exterior;

- Las variables del campo de propiedades físicas exhiben anisotropía tridimensional progresiva;

- El espacio que constituye el medio acuífero se presenta rigurosamente jerarquizado (Tabla 1);

- La existencia de dominios de flujo inherentes a cada espacio;

- Un campo de propiedades físicas que se define y estructura para cada espacio;

- La fuerte influencia del efecto del factor de escala sobre el campo de propiedades físicas (Fig. 2)

- La elevada dependencia respecto al tiempo de las propiedades físicas;

- La formación y desarrollo de estructuras auto reguladas de disipación de energía;

- Un cierto valor de inercia del sistema (memoria), esto es, la modulación de las respuestas a los estímulos inducidos natural o artificialmente, que dependen, sobre todo, de su estado inicial, y, finalmente,

- La irreversibilidad del proceso de carsificación, es decir, su evolución unidireccional.

Tabla 1. Espacios constitutivos del sistema cársico (según Molerio, 1985a)

Categorías			Denominación del espacio	Volumen con relación al sistema	Longitud característica	Régimen de flujo	Diámetro
Espacio cársico		Macro discontinuidades	1. Cavernas	V ≤ 1%	km.	No lineal de alta velocidad	m
	Espacio agrietado no cársico		2. Grietas y discontinuidades en general	1% ≤ V ≤ 5%	km.	No lineal de alta velocidad	mm - m
			- Fallas		km.	Lineal	m
			- Diaclasas		km.	Lineal	mm-cm
			- Planos de estratificación.		km.	Lineal	cm
			- Planos de esquistosidad.		cm.	Lineal	mm
	Espacio poroso	Micro discontinuidades	3. Poros de la matriz rocosa microfracturas	5% < V ≤ 40%	mm.	No lineal de baja velocidad (microflujo)	≤ 0,1 mm
			- Exfoliación		mm.		0.1 mm
			- Esquistosidad		mm.		10 μ
			- Clivaje		mm.		0,1 μ
			4. Matriz sólida	V > 40%	-	No lineal de baja velocidad (no flujo)	< 1 μ

Fig. 2. Curva de efecto del factor de escala sobre la estructura del campo de propiedades físicas

del karst (modificado de Kiraly, 1975; 1978).

La más importante ecuación de base es la Ecuación de Darcy generalizada para el medio poroso anisotrópico, que se expresa del modo siguiente:

Donde, q = vector de velocidad de Darcy [L/T], K = tensor de conductividad hidráulica [L/T] y ∇h = gradiente de carga hidráulica [L/L]. Habida cuenta que el tensor Kij es simétrico y puede representar anisotropía debida a fracturas, estratificación o conductos preferenciales, entonces en notación tensorial, queda:

$q_i\text{=}\text{-}{K}_{\mathit{ij}}\frac{\partial h}{\partial x_j}$

La Ecuación de Conservación de Masa (Continuidad) que en un medio poroso anisotrópico se expresa como

$S_a\frac{\partial h}{\partial t}\text{=}\nabla \bullet \left(K\bullet \nabla h\right)\text{+}q$

Donde h es la carga hidráulica [L]; Sa, el almacenamiento especifico [L-1]; K, el tensor de conductividad hidráulica (simétrico de segundo orden); q, término fuente/sumidero [T-1], $\nabla h$  y $\nabla$  son los operadores divergencia y gradiente, respectivamente.

La Ecuación de flujo en conductos (Ley de Darcy-Weisbach o Manning), unidimensional, régimen no lineal -turbulento- o laminar –lineal- es,

$h_f\text{=}f\frac{L}{D}\frac{v^2}{2g}$

Donde hf, es la pérdida de carga [L]; f, el factor de fricción adimensional de Darcy-Weisbach; L, la longitud del conducto [L]; D, diámetro hidráulico [L]; v, la velocidad media del flujo [LT-1] y g, la aceleración gravitacional [LT-2]

La Ecuación de Manning para flujo en conductos o canales es,

$v\text{=}\frac{1}{n}{R}_h^{\frac{2}{3}}{S}^{\frac{1}{2}}$

En la que v, velocidad media [LT-1]; n, coeficiente de rugosidad de Manning [TL-1/3]; Rh, radio hidráulico (área mojada/perímetro mojado) [L] y S, la pendiente de la línea de energía adimensional; representada por la siguiente ecuación:

$S\text{=}\text{-}\frac{\mathit{dh}}{\mathit{ds}}$

La pérdida de carga se expresa como:

$h_f\text{=}\frac{n^2{v}^{2}L}{R_h^{\frac{4}{3}}}$

Las ecuaciones de Darcy-Weisbach y Manning son, en su esencia, relaciones escalares unidimensionales para calcular pérdidas de carga en la dirección del flujo. No son leyes constitutivas tensoriales como la de Darcy para medios porosos. Sin embargo, para su implementación en modelos 2D o 3D (por ejemplo, en modelos de elementos finitos o volúmenes finitos para simular flujo en ríos o grandes conductos), se pueden escribir en una forma vectorial o tensorial proyectada que generaliza el concepto. La pérdida de carga hf​ es una magnitud escalar, pero el gradiente de pérdida de carga (la fuerza resistiva por unidad de peso) es un vector que se opone al vector de velocidad v.

La forma vectorial generalizada para un flujo 3D se escribe como un término fuente de momento en la ecuación de movimiento:

$S_f\text{=}\text{-}\frac{f}{D}\frac{\left|v\right|}{2g}v$

Para el gradiente de carga,

$\nabla h_f\text{=}\text{-}\frac{f}{D}\frac{\left|v\right|}{2g}\frac{v}{\left|v\right|}\text{=}\text{-}\frac{f}{D}\frac{v\left|v\right|}{2g}$

Que en notación de componentes cuasi tensoriales se escribe:

${\left(\nabla h_f\right)}_i\text{=}\text{-}\frac{f}{D}\frac{v_i\sqrt{v_j{v}_j}}{2g}$

Donde vjvj es la suma sobre j; esto es, la magnitud al cuadrado. De este modo, la forma vectorial generalizada para el gradiente de pérdida de carga es:

$\nabla h_f\text{=}\text{-}\frac{n^2}{R^{4/3}v\left|v\right|}$

Y en notación de componentes:

${\left(\nabla h_f\right)}_i\text{=}\text{-}\frac{n^2}{R^{4/3}}{v}_i\sqrt{v\left|v\right|}$

Rh sigue siendo un escalar que debe calcularse en función de la geometría local. En modelos 2D de aguas someras (profundidad promediada) Rh se aproxima a la profundidad del flujo h.El acoplamiento (coupling) en un modelo de doble porosidad/doble permeabilidad para el flujo en un conducto de gran diámetro, como una galería subterránea (Fig. 3) representada como un elemento 1D expresaría Manning como:

$\frac{\partial Q}{\partial t}\text{+}\dots \text{-}\mathit{gA}\frac{n^2}{R^{4/3}}Q\left|Q\right|\text{-}\mathit{gA}\frac{\partial h}{\partial s}$

Fig. 3. Sección de galería subterránea. Los caudales que circulan por el espacio de las cavernas tanto en flujo hipodérmico como en la zona saturada se expresan siempre en régimen no lineal (turbulento) (Foto del autor).

El medio poroso circundante puede simulare con las ecuaciones de Darcy/Continuidad acopladas por un término de intercambio Γ. Pese a su complejidad física y matemática, el modelo de Doble Porosidad/Doble Permeabilidad (MDP) es la aproximación más simple para un modelo natural que, en el mundo real, integra cuatro espacios (ver Tabla 1). En este caso, la ecuación para la matriz MDP es la siguiente:

$S_m\frac{\partial h_m}{\partial t}\text{=}\nabla \bullet \left(K_m\bullet \nabla h_m\right)\text{-}\Gamma$

Y para los conductos/grietas,

$S_f\text{=}\frac{\partial h_{\mathit{mf}}}{\partial t}\text{=}\nabla \bullet \left(K_f\bullet \nabla h_f\right)\text{-}\Gamma$

La Ecuación General de Transporte con tensor de dispersión se escribe del modo siguiente:

$\frac{\partial C}{\partial t}\text{=}\nabla \bullet \left(D\bullet \nabla C\right)\text{-}v\nabla C\text{+}R$

Donde D es el tensor de dispersión hidrodinámica, que en el karst es altamente anisótropo.

2.2. Presentación de casos de estudio

Algunas aplicaciones de este recurso matemático son especialmente importantes para la gestión adecuada de los recursos hídricos en el karst que, con mucho, es la fuente de aguas subterráneas más importante de Cuba.

El karst constituye el medio acuífero más importante de Cuba y como espacio geográfico, domina el territorio nacional (Fig. 4). Es un geosistema particular que se extiende sobre una superficie aproximada de 66 763 kilómetros cuadrados, casi el 66% del área total del país, de los cuales 62 780 km2 corresponden a tierra firme y el resto al karst de los cayos es isletas del archipiélago (Molerio, 1974).

Fig. 4. Esquema general de distribución de las áreas cársicas de Cuba (basado en la distribución geográfica de las rocas carsificables, básicamente carbonatadas, según el Atlas Nacional de Cuba, 1989).

En estos territorios, de los que 66 495 km2 corresponden a llanuras, se concentran numerosos núcleos poblacionales urbanos y rurales, industrias, centros educacionales, y de turismo y recreación, instalaciones agrícolas y ganaderas, y otros que, en mayor o menor grado, requieren del empleo racional y suministro de agua potable cualitativa y cuantitativamente correspondientes con sus necesidades actuales y perspectivas. Estas llanuras cársicas acuíferas, forman karsts litorales, permanentemente amenazadas por el avance, tierra adentro, de la intrusión marina.

De los 6,3 km3 de recursos de agua subterráneas de que dispone, 80 % se encuentra en estos territorios donde se encuentra el 95% de las cuencas hidrogeológicas cubanas. Por su relación con los ecosistemas que crea, sostiene o interactúa, es transversal a la mayor parte de la sociedad, la economía y el medio ambiente. También en él se encuentran yacimientos de petróleo, importantes minerales útiles, la mayor parte de los suelos agrícolas del país y prácticamente las cuatro quintas partes de la población del país se sirve o aprovecha sus recursos naturales.

Se han seleccionado, para su revisión, varios casos de estudio que caracterizan los siguientes problemas:

- La identificación de la intrusión marina y el pronóstico de avance

- La corrección de la dirección de flujo subterráneo

- Los procesos de transporte de contaminantes (particularmente de los procesos de nitrificación y de líquidos de fase no acuosa, incluidas la migración vertical ascendente y descendente)

- La simulación de los procesos de espeleogénesis y desarrollo de cavernamiento epigenético

2.2.1. Identificación de la intrusión marina y el pronóstico de avance

El problema ha sido examinado con mucho detalle recientemente por el autor (Molerio 2025). El caso de estudio (Molerio, 2006) trata del efecto de las mareas en la modificación del campo de transmisividad en un acuífero litoral y el modo en que ello enmascara tanto las direcciones de flujo como el movimiento de la intrusión marina (Fig. 5).

Fig. 5. Cambios en la dirección y gradiente de flujo como resultado de la penetración cíclica tierra adentro de la marea en un acuífero cársico litoral (Molerio, 2006)

Muchos de los problemas de la ineficacia de la gestión de la intrusión marina se deben a que la conceptualización del problema es imperfecta. Ciertas premisas de la aproximación hidrostática dominante pudieran resultar aplicables bajo determinadas condiciones, pero en realidad es un modelo sumamente limitado para la gestión de los recursos hídricos en acuíferos costeros que conduce a conclusiones falsas y a la proyección de obras de defensa sumamente costosas y de dudosa efectividad aun en medios isotrópicos. En medios anisotrópicos el problema se complica aún más y este es el caso de la absoluta mayoría de los acuíferos costeros de Cuba y muchas regiones como Yucatán (México), Jamaica y Puerto Rico, así como de todas las islas carbonatadas del Caribe, donde constituyen la única fuente de agua segura para todos los usos.

El movimiento natural e inducido de fluidos dependientes de la densidad y los procesos de intercambio que tienen lugar en el karst costero aun requieren de soluciones matemáticas que resuelvan satisfactoriamente el problema y sostengan teóricamente la proyección y adopción de medidas –estructurales y no estructurales- de gestión eficiente. Desde la selección de sistemas de explotación, obras de recarga artificial, diseño y operación de sistemas de monitoreo y sustenten medidas más generales respecto a los efectos negativos del cambio climático.

En acuíferos cársicos, por razón de la heterogeneidad y la anisotropía del campo de propiedades físicas (Palmer, 1999), el problema es aún mayor cuando se pretenden soluciones exactas al pronóstico de la posición de la interfase. La heterogeneidad produce un retroceso de la penetración de la cuña, mientras aumenta en el ancho y la pendiente de la zona de mezcla. La forma de la interfase y el flujo de agua salada son función de la distribución de las permeabilidades en cada realización. Sin embargo, los resultados de la penetración de la cuña y el ancho de la zona de mezcla no presentan grandes fluctuaciones. En caso de heterogeneidad moderada, ambas variables pueden ser reproducidas de forma satisfactoria por un medio homogéneo con dispersión local o efectiva (Abarca, 2006).

El problema está lejos de ser sencillo y ha sido abordado desde dos grandes aproximaciones: a) la hidrostática de Ghyben y Herzberg con modificaciones como la de Lusczinsky que considera ambos fluidos no miscibles separados por una interfase nítida (sharp interface); b) la dinámica, de Hubbert y Bear que, acompañadas de un esquema de intercambio de fluidos de densidad variable es la que representa el mundo real. Soluciones complementarias se encuentran en las interfases variadas de Glover, las soluciones exactas de Van der Veer y las complementarias de Vacher. En ambos casos, para ciertas condiciones de borde se aplica Dupuit-Forcheimer.

Fig. 6. Corte típico de la relación entre las aguas dulces terrestres y salinas marinas en un acuífero litoral (Molerio, 1994)

Ese equilibrio puede variar y de hecho varía en varias escalas de tiempo. El límite entre las aguas dulces y marinas NO es nítido. Ambos líquidos son miscibles y, por lo tanto, el propio concepto carece de validez. En lo absoluto es claro, bien definido y abrupto y la curva de 1 g/L NO ES TAMPOCO indicadora de la posición de la intrusión marina (Figs. 6-8). En el corte vertical existe una gradación de mineralizaciones que dependen de la composición litológica, las barreras hidráulicas y geológicas (acuíferos libres, confinados, semiconfinados, retardados, percolante, acuicludas…) y una distribución vertical diferenciada del campo de velocidad de las aguas asociada a la heterogeneidad y anisotropía de las rocas carsificadas y al relieve y estructura geológica de los acuíferos costeros.

El paso de la zona ocupada por las aguas de mar a la zona ocupada por las aguas dulces ocurre en una zona de transición que depende absolutamente de la dispersión hidrodinámica de la materia disuelta. El coeficiente de dispersión D es un tensor simétrico de segundo orden, resultado del producto escalar del tensor de cuarto orden, dispersividad del medio con el acoplamiento de velocidad. Como mezcla macroscópica de las partículas de fluido, el proceso de dispersión está influenciado por la no uniformidad y la anisotropía del medio (Molerio, 1997). Puede ser muy ancha, como en los acuíferos cársicos donde también se presenta de manera irregular en el corte vertical dependiendo de los flujos hacia y desde el mar asociados a la diferente carga hidráulica de los niveles de cavernamiento en profundidad como es el caso de los acuíferos carbonatados Neógeno-Cuaternarios cubanos y comportarse como una línea muy delgada con en acuíferos de muy baja permeabilidad como aquellos en materiales terrígenos o vulcanógeno sedimentarios como los del Cretácico – Paleógeno cubano.

La hipótesis de que el agua de mar permanece estática de Ghyben-Herzberg no es válida, pues el flujo del agua de mar produce una pérdida de energía que se traduce en una menor penetración de la cuña de agua salina. Por lo tanto, la profundidad de la interfaz se subestima (Pool, Abarca y Carrera, 2007). De igual modo, considerar que el problema se desarrolla en el contacto entre las aguas dulces y las del mar en el borde del litoral y que son válidas las presunciones de Dupuit – Forcheimer para flujo horizontal puede conducir, y de hecho conduce, a resultados inadecuados.

El efecto del desarrollo de las estructuras de conducción del agua en medios anisotrópicos; planos de estratificación, grietas y cavernas produce un efecto de deformación de la cuña de intrusión marina que depende de su distribución, conexión física, geometría y distribución del campo de velocidad del fluido, afectada esta última también por el material que eventualmente rellene o colmate estos espacios.

La penetración de la cuña de intrusión marina depende también de la profundidad del fondo del acuífero (Abarca, 2006) y de la estructura tectonofacial del acuífero. En efecto, los cambios en la forma de la cuenca sedimentarias y las variaciones faciales horizontales y verticales condicionan tanto desarrollos diferenciados de carsificación y el cavernamiento en acuíferos carbonatados, como diferentes distribuciones del espacio de las grietas o de la morfología del dominio de flujo en sistemas deltaicos.

En el caso de estudio del campo de dolinas de El Cayuelo (Molerio, 2006), en virtud de las propiedades tensoriales, se identificó una dirección de la conductividad hidráulica o de las transmisividades en dirección de la velocidad y otra, en dirección del gradiente. Las componentes principales de las máximas transmisividades se encuentran entre los 20 y 45º, en tanto las componentes mínimas se encuentran entre los 70 y 80º. Es evidente, entonces que los vectores de flujo q, y de gradiente hidráulico J, no son paralelos, por lo que puede existir flujo en una dirección diferente a la del gradiente hidráulico. En medios anisotrópicos, el campo de velocidades en el sistema acuífero puede escribirse de la forma siguiente para las tres coordenadas (Fig. 7):

$v_x\text{=}\text{-}{k}_x\frac{\delta h}{\mathit{\delta x}}$

$v_y\text{=}\text{-}{k}_y\frac{\delta h}{\mathit{\delta y}}$

$v_z\text{=}\text{-}{k}_z\frac{\delta h}{\mathit{\delta z}}$

Fig. 7. Proyección estereográfica del campo de transmisividad en la dolina del repositorio.

Para el caso general de anisotropía, la Transmisividad es un tensor simétrico de dos componentes, uno en la dirección de la velocidad y otro en la dirección del gradiente. Para flujo bidimensional, entonces ocurre que

$T\text{=}\left|\begin{array}{c}{T}_{\mathit{xx}}{T}_{\mathit{xy}}\\ T_{\mathit{xy}}{T}_{\mathit{yy}}\end{array}\right|$

y, para los ejes principales $\xi y\eta$ , entonces

$T\text{=}\left|\begin{array}{c}\text{\&}{T}_{\mathit{\xi \xi }}0\\ \text{\&}0T_{\mathit{\eta \eta }}\end{array}\right|$

Donde $T_{\mathit{\xi \xi }}{T}_{\mathit{\eta \eta }}$ son las transmisividades máximas y mínimas, respectivamente. La Fig. 8 muestra el elipsoide de direcciones construido con los valores derivados de los ensayos puntuales de caudal.

Fig. 8. Representación gráfica de la elipse de transmisividad.

En el campo ensayado, $T_{\mathit{\xi \xi }}$  y $T_{\mathit{\eta \eta }}$  no son irrestrictamente ortogonales entre sí, de manera que se está en presencia de una anisotropía axisimétrica o transversal. Esta anisotropía exhibe una fuerte componente en dirección 45º, donde se obtuvo el mejor ajuste de interpolación mediante kriging lineal con pendiente unitaria.

En el elipsoide de direcciones de T se construyeron las componentes de T en dirección del gradiente para el caso de anisotropía transversal, en la que se derivó una componente en dirección 135º-315º. Pueden pronosticarse dos direcciones fundamentales de flujo subterráneo en el sistema ensayado, una en dirección 20-45º 200-225º) y otra en dirección 135-315º. Cambios locales de flujo deben esperarse en las direcciones 70-80º (250-260º), pero asociadas a campos muy locales.

La transmisividad, entonces, parece estar muy bien estructurada, ya que la dirección de la componente del gradiente ( $T_{\mathit{\eta \eta }}$ ) satisface estrictamente la propiedad de ortogonalidad respecto a la dirección del flujo ( $T_{\mathit{\xi \xi }}$ ). Esto es de la mayor importancia, en tanto descarta la presencia de campos ficticios de anisotropía asociados con la transversalidad de la anisotropía (Bear, Zaslavsky e Irmay, 1968) o con el efecto de escala sobre la estructura del campo de transmisividad (Kiraly, 1975, Molerio 1984, 1986).

La solución de la ecuación

$T\text{=}\left|\begin{array}{c}{T}_{\mathit{\xi \xi }}0\\ 0T_{\mathit{\eta \eta }}\end{array}\right|$

puede plantearse en términos de la dirección de referencia que, en este caso, equivale a la diferencia modular entre las direcciones de $T_{\mathit{\xi \xi }}$  y $T_{\mathit{\eta \eta }}$ , tal que:

$T_{11}\text{=}{T}_{\mathit{\xi \xi }}{\cos }^2\alpha \text{+}{T}_{\mathit{\eta \eta }}\mathit{se}{n}^2\alpha$

$T_{22}\text{=}{T}_{\mathit{\xi \xi }}\mathit{se}{n}^2\alpha \text{+}{T}_{\mathit{\eta \eta }}{\cos }^2\alpha$

$T_{21}\text{=}{T}_{12}\text{=}\text{-}\frac{1}{2}(K_{\mathit{\xi \xi }}\text{-}{K}_{\mathit{\eta \eta }})\mathit{sen}2\alpha$

Y, resolviendo adecuadamente, se obtienen los siguientes resultados:

$T_{11}\text{=}{T}_{\mathit{xx}}\text{=}\frac{17312m^2}{d}$

$T_{22}\text{=}{T}_{\mathit{yy}}\text{=}\frac{5379m^2}{d}$

$T_{21}\text{=}{T}_{12}\text{=}{T}_{\mathit{xy}}\text{=}{T}_{\mathit{yx}}\text{=}\text{-}4570m^2\text{/}d$

Por tal motivo, la recarga o descarga (s) instantánea del acuífero se expresan, para las direcciones principales de anisotropía (Papadopoulos, 1967), por la siguiente expresión, que toma en cuenta el almacenamiento S:

$s\text{=}\frac{2.303Q}{4\pi \sqrt{T_{\mathit{\xi \xi }}\mathit{T\eta \eta }}}{\log }_{10}\left[\frac{2.25t}{S}\left(\frac{T_{\mathit{\xi \xi }}{T}_{\mathit{\eta \eta }}}{T_{\mathit{\xi \xi }}{\eta }^2\text{+}{T}_{\mathit{\eta \eta }}{\xi }^2}\right)\right]$

2.2.2. Corrección de la dirección del flujo subterráneo

El trazado de hidroisohipsas en medios anisotrópicos y, con ello en el karst, tiene un valor de uso limitado por razón de la heterogeneidad del medio y, fundamentalmente de su anisotropía. El problema es conocido y ha sido descrito por los autores clásicos hará pronto 100 años. Sin embargo, sin corrección alguna continúa utilizándose y peor aún, se evalúan los recursos de agua subterránea en base a ello, se definen las zonas de protección sanitaria de los pozos de abasto y se diseñan obras de recarga artificial o se explotan las aguas subterráneas.

Fig. 9. Representación vectorial de la Ley de Darcy a) medio isotrópico; b) medio anisotrópico (Hubbert, 1940)

Un libro de texto clásico, como el de Todd y Mays (2005) presenta una opción de transformación a un medio continuo que tampoco es utilizada en la práctica común de gestión de los acuíferos cársicos. Sin discutir dónde y cuándo es válida la reducción al continuo (véanse los trabajos de Bear y Maasland, 1957) señala que “en medios anisotrópicos las líneas de flujo y las equipotenciales no son ortogonales excepto cuando el flujo es paralelo a una de las direcciones principales. Para calcular los flujos para esta situación, los límites de una sección de flujo deben ser transformados de modo que se pueda obtener un medio isotrópico”. Lo cual tampoco es común que se haga.

De acuerdo con Kiraly (1975, 1978) en un medio sencillo de doble porosidad/permeabilidad condicionada el flujo por el patrón de agrietamiento, la conductividad hidráulica tensorial K se expresa como:

La fig. 10 muestra las direcciones reales de flujo derivadas del procesamiento anterior, en la dirección de la velocidad (Kν) y en la dirección del gradiente hidráulico (KI), diferentes a las que pueden deducirse de las hidroisohipsas. En este caso particular, nótese que el vertimiento sistemático de residuales azucareros en la Gran Dolina de Aguacate (A) que ha afectado históricamente la calidad de las aguas en los Manantiales Bello, en Matanzas (B) se acomoda perfectamente a un máximo de Kν asociado al sistema dominante de grietas (líneas discontinuas) y, muy probablemente, del cavernamiento regional. Esta conexión ha sido, del mismo modo, comprobada por ensayos con trazadores fluorescentes e isotópicos; queda pendiente la exploración directa del karst subterráneo.

Fig. 10. Líneas de flujo definidas por la corrección tensorial de la conductividad hidráulica K en la dirección de flujo (Kv, dominante) y del gradiente hidráulico (KI, complementaria). Comentarios en el texto.

2.2.3. Procesos de transporte de hidrocarburos en aguas subterráneas

El recorrido de la eventual migración de aguas petrolizadas desde un repositorio y la proyección de las medidas de protección de las aguas superficiales, subterráneas y, eventualmente, marinas, solamente puede pronosticarse si se conocen las direcciones de flujo subterráneo. En el karst la dispersividad no es una propiedad escalar, sino una cantidad tensorial, escala-dependiente y dominio-específica, lo que exige modelos híbridos (continuo + discreto) para una predicción realista.

La descripción de la dispersión de hidrocarburos (petróleo y aguas petrolizadas) en un acuífero cársico, considerando la naturaleza tensorial de la dispersividad, es un problema complejo debido a la doble (o múltiple) porosidad, la fuerte anisotropía y la no linealidad de estos medios. La dispersión hidrodinámica en medios porosos se describe comúnmente mediante la ley de Fick generalizada, donde el coeficiente de dispersión (D) es un tensor de segundo orden que depende de la dispersividad intrínseca del medio (α), que también es un tensor y la velocidad promedio del flujo (v). Pero el petróleo, como Líquido de Fase No acuosa (LFNA-NAPL) no sigue una ley de dispersión Fickiana clásica; su avance es por gradientes de presión capilar y gravedad y la dispersión es más bien una saturación irregular a lo largo de las líneas/conductos de flujo.

En el karst, estas direcciones no están alineadas con los ejes cartesianos, sino con la dirección de los conductos principales, fracturas y estratificación. Por tanto, se requiere una transformación tensorial mediante una matriz de rotación. Además de las propiedades mencionadas en este artículo para el karst, en el caso particular del petróleo (NAPL - LNAPL o DNAPL) forma fases separadas, con flujo multifásico. Las aguas petrolizadas (emulsiones o fase disuelta) se transportan en la fase acuosa, pero su dispersión está modulada por la retención interfacial y la disolución cinética. Siguen leyes del transporte advectivo-dispersivo, pero con coeficientes de dispersión dependientes de la concentración (viscosidad y densidad variables), de ahí que la dispersividad aparente puede aumentar debido a la heterogeneidad de la distribución de NAPL residual.

En un medio anisotrópico, la dispersividad (α) es un tensor de cuarto orden que relaciona el tensor de dispersión (D) con el tensor de velocidad. Para simplificar, en la práctica se suele reducir a un tensor de segundo orden simétrico (3×3 en 3D) con componentes longitudinal, transversal y vertical.

Un enfoque es aplicar un modelo de doble continuo (o múltiple continuo). Para la matriz/fracturas (medio continuo equivalente) se define un tensor de dispersividad efectiva αeff mediante upscaling (promediación de REV) y se considera anisotropía derivada de la estructura de fracturas (ej.: tensor elíptico). Para los conductos kársticos (elementos discretos o medio continuo preferencial), se toma en cuenta que la dispersión es esencialmente unidimensional a lo largo del conducto, con αL >> αT. La dispersión transversal puede ser despreciable si el conducto está bien desarrollado, pero crucial en zonas de mezcla, como puede ser la confluencia de galerías subterráneas.

El modelo debe considerar la naturaleza tensorial de la dispersividad en medios altamente anisotrópicos y la posibilidad de flujo multifásico. Para un contaminante disuelto en agua en un medio cársico con conductividad hidráulica anisotrópica, la ecuación de transporte para un soluto en fase acuosa con dispersividad tensorial (medio anisotrópico):

Donde Dij es el tensor de dispersividad hidrodinámica efectiva:

En medios cársticos fracturados o con conductos, αL​ y αT​ son tensores de segundo orden debido a la anisotropía estructural, no solo escalares.

Para fase libre de hidrocarburos (no acuosa), se usa el enfoque de flujo multifásico en medios porosos que, adaptado a la doble porosidad del karst, se expresa como, la conservación de masa para la fase de hidrocarburo libre (Fase β)

Con la Ley de Darcy multifásica

Donde, Sβ, la saturación de la fase hidrocarburo, krβ​, permeabilidad relativa (función de Sβ​), K, tensor de permeabilidad intrínseca del medio cársico (anisótropo) y Pβ​, la presión en la fase β. En el karst, K puede ser un tensor diagonal no isótropo o completo si hay anisotropía cruzada

Para un componente disuelto desde la fase hidrocarburo, con intercambio de fases y dispersión tensorial:

Con kdis, Coeficiente de transferencia de masa entre fases [T⁻¹], Cs, la concentración de saturación del soluto en agua [M·L⁻³] y D, tensor de dispersividad. La forma más general para un medio con anisotropía estructural (fracturas, conductos preferenciales) en donde βijkℓ​ es un tensor de cuarto orden como ya ha sido definido, es la siguiente:

que relaciona la velocidad con la dispersión. En la práctica, se simplifica a:

2.2.4. Simulación del desarrollo del cavernamiento

La modelación de la espeleogénesis es una tarea extremadamente compleja debido, entre otras, a las siguientes causas: el elevado número de variables involucradas, la naturaleza estocástica de muchas de estas variables y de los procesos que en ella intervienen, la presencia de procesos no lineales, la ausencia de ecuaciones que describan algunos de los mecanismos que intervienen en el proceso general y la ausencia o formulación incompleta de un modelo conceptual susceptible de ser tratado física y matemáticamente.

El tema ha recibido sistemática atención y, de hecho, comenzó con las primeras teorías sobre el origen de las cavernas enunciadas hace más de tres siglos. Modernamente (véanse generalidades en Gunn, 2004 y White, 2009) y bajo diferentes aproximaciones, se han ocupado del tema, (Curl, 1971; Eraso, 1975, 1982, 1985/86; Gunn, 1986; Eraso y Fernández Rubio, 1990; Paredes et al., 1993; Eraso y Domínguez, 2004; Ewers, 1976, 1978, 1982; White, 1976; Gale, 1984; Lauritzen et al., 1985; Smart, 1988; Oreskes et al, 1994; Fields y Pinsky, 2000; Palmer 1981, 1988, 1991, 1995, 1998; Dreybrodt 1987, 1988, 1989, 1990, 1992, 1993, 1995, 1996, 1998; Dreybrodt y Buhmann, 1991; Dreybrodt y Siemers, 1997; Klimchouk et al., 2000; Jeannine, 2001; Dreybrodt, Gabrovsek y Romanov, 2005; Dreybrodt y Svensson, 1989; Ford, 1971, 1988; Ford y Ewers, 1978; Hanna y Ramajaran, 1998; Gabrovsek y Dreybrodt, 2001; Howard y Groves, 1995; Klimchouk, 2003, 2005; Klimchouk et al., 2000; Kovacs, 2003; Molerio 1982a, 1985a,1985c, 1986a, 1986b, 1988, 1989b, 1989c, 1989d, 1990, 1992c, 1992d, 1995b, 1996, 1997, 2000, 2003a, 2003b; Molerio Guerra y Flores 1984; Mylroie y Carew, 1990, 1995; Palmer, 1975, 1988, 1989, 1991, 2000, 2007; Palmer, Palmer y Sasowsky, 1999; Reinmann, 2012; Reinmann y Hill, 2009; Stephanssson, Hudson y Lanaru, 2004; Szymczak y Ladd, 2006, 2009; Werner, 2004; White, 2002).

El modelo desarrollado por el autor (MTDK: Modelo Termodinámico del Desarrollo del Karst) se basa en los siguientes principios, enunciados inicialmente en 1986 (Molerio, 1986a, 1986b):

- Las galerías subterráneas son espacios lineales y no planares o areales;

- Las topologías lineales por lo común se desarrollan a lo largo de las intersecciones entre superficies;

- En cuanto al desarrollo de las cavernas estas superficies son de dos tipos: la zona de máxima concentración de solvente y la zona de máxima concentración de flujo;

- La superficie máxima de concentración de solvente (MCS) es generalmente horizontal o subhorizontal y depende de la evolución geoquímica del medio, la fuente de aporte y el tiempo de residencia de las aguas en el macizo;

- Las superficies de máxima concentración de flujo (MCF) suelen estar fuertemente inclinadas y, con menor frecuencia, pueden ser completamente horizontales;

- En la intersección de las superficies MCF y MCS se encuentra no solamente el mayor volumen de fluido sino la más elevada concentración de solvente, de ahí que en ella ocurre la mayor probabilidad de disolución y por ello, de desarrollo de cavernas;

- El trabajo que se realiza en el punto de intersección no da lugar al equilibrio térmico. El intercambio de masa y energía conduce al desarrollo progresivo de un sistema abierto en el cual, la entropía crece a partir de un instante inicial to en que el sistema deja de ser cerrado;

- Las superficies MCF se encuentran en la dirección de la componente de conductividad hidráulica en el sentido de la velocidad. Flujo lateral se encuentra en la dirección de la componente de gradiente hidráulico, de manera que la máxima probabilidad de desarrollo puede determinarse conociendo éstas, lo que significa que el desarrollo de la red de cavernas no es un fenómeno aleatorio y por tanto, puede predecirse;

- Conociendo la orientación de las superficies más favorables para el desarrollo de las redes de conductos es imprescindible entonces, determinar la dirección en que ocurre el proceso de excavación. De acuerdo con el segundo principio de la termodinámica, ella ocurre en la dirección del máximo incremento de la entropía;

- Cada espacio del universo cársico, incluidas las cavernas, se desarrolla según los principios de la termodinámica de los procesos de no equilibrio y a ella corresponden valores de la función de disipación de energía crecientes que se expresan como una sumatoria a partir del centro de gravedad del mismo.

- El acoplamiento termodinámico, para la definición de la función de disipación de energía, se plantea a partir de las relaciones clásicas de Onsager que relacionan fuerzas y flujos entre las componentes de calor y fluido, difusión y afinidad química. Estas últimas son acoplamientos en sí mismas y han sido excluidas momentáneamente. Actualmente se estudia una vía adecuada para reducir el número de componentes de la difusividad de 81 a 36 componentes, aceptando que represente un tensor de cuarto orden.

Estos principios básicos controlan las reacciones químicas y los procesos físicos que tienen lugar y que conducen al crecimiento tridimensional de las cavernas. El sistema, en su conjunto, se considera como una estructura autorregulada no lineal sometida a acciones aleatorias (Netushil, 1987; Yevjevich, 1971). Bajo tales principios han sido estructurados los conceptos originalmente expuestos por Egemeier, 1969; Groenevelt y Belt, 1969; Glansdorff y Prigogine, 1971; Kiraly, 1975; Carnahan, 1976 y Eraso, 1982.

El acoplamiento termodinámico, para la definición de la función de disipación de energía, se plantea a partir de las relaciones clásicas de Onsager que relacionan fuerzas y flujos entre las componentes de calor y fluido, difusión y afinidad química. Estas últimas son acoplamientos en sí mismas y han sido excluidas momentáneamente. Actualmente se estudia una vía adecuada para reducir el número de componentes de la difusividad de 81 a 36 componentes, aceptando que represente un tensor de cuarto orden.

El origen de las cavernas está gobernado por un balance de masas tal que la tasa de crecimiento de los conductos, como consecuencia de la remoción de masa de las paredes de la cavidad o el conducto, es igual a la tasa de transporte de masa en solución. El acoplamiento termodinámico, para la definición de la función de disipación de energía, se plantea a partir de las relaciones clásicas de Onsager que relacionan fuerzas y flujos entre las componentes de calor y fluido, difusión y afinidad química.

El algoritmo general ha sido desarrollado a partir de los principios de la termodinámica de no equilibrio, que se considera esencial para determinar la dirección de los procesos de cavernamiento. El algoritmo se basa en un conjunto de ecuaciones de control que describen: a) la continuidad macroscópica del campo de propiedades físicas; b) los mecanismos de triggering en el sistema físico-químico; c) la competencia entre diferentes líneas de flujo y d) la dirección de la evolución en tiempo y espacio, de los procesos de desarrollo del cavernamiento.

Los procesos de transporte de masa, momento y energía entre medios continuos equivalentes fueron resueltos para cada uno de los espacios involucrados, en términos de la derivación de un conjunto de ecuaciones que describieran las correspondientes funciones de transporte. Del mismo modo, el trabajo que tiene lugar en el sistema, y que conduce a la formación de estructuras autorreguladas de disipación de energía, fue examinado a partir de la derivación de un conjunto de ecuaciones que describen la función disipativa de los espacios.

Sin embargo, en primera aproximación, cada forma del relieve posee un valor de la función de disipación de energía controlada, inicialmente, por la tasa volumétrica de producción de entropía. Del mismo modo, cada espacio que integra el sistema debe estar caracterizado por un entorno característico de valores para tal función. Debe suceder, entonces, que la función sea continua en tiempo y espacio, por cuanto cada nueva forma contiene, en el sistema, mayor entropía que la precedente. De aquí se avanza a la caracterización del sistema combinando las ecuaciones de balance de masas, momento y energía con la ecuación de Gibbs para lograr una formulación adecuada que exprese el balance de entropía, En todo caso, se supone la simetría del campo tensorial de los coeficientes fenomenológicos involucrados.

Así, se obtuvo una ecuación general para definir el cambio de entropía en el sistema en función de la diferencia de potenciales y la afinidad química de las reacciones fundamentales. Las fuerzas y flujos considerados en el balance termo hidrodinámico básico son: a/ transporte de calor (ley de Fourier); b/ flujo volumétrico (ley de Darcy); c/ transporte dispersivo-difusivo (ley de Fick); d/ conductividad eléctrica (ley de Ohm) y e/ la afinidad química de las reacciones, enlazadas por las funciones de transferencia y disipación de energía mediante acoplamiento termodinámico múltiple.

Los resultados obtenidos en la simulación de sistemas reales fueron altamente promisorios. Sin embargo, las desviaciones respecto al modelo natural se derivaban de numerosas fuentes de incertidumbre cuya discriminación resulto una tarea ardua. La más importante de estas se reducía a distinguir las condiciones que provocasen que un sistema inicialmente cerrado o aislado, reversible, sin coacciones exteriores, se transformase en un sistema termodinámico abierto, que estuviese caracterizado por las propiedades definidas en el modelo conceptual del karst. En este sentido, se orientó la investigación hacia la caracterización de acciones aleatorias exteriores, o producidas por el sistema, en términos de la adaptabilidad de este para filtrar tales estímulos, definir su efecto en el caso de provocar fluctuaciones termodinámicas, diferenciar la estacionalidad de las señales aleatorias y tratar de resolver la respuesta del sistema en la dirección de los niveles crecientes de entropía. Uno de los aspectos básicos involucrados en el cambio de tipo termodinámico del sistema lo constituyen los mecanismos de triggering cinético.

Cualquier estado macroscópico es función del estado microscópico en virtud de la reversibilidad de las ecuaciones del micromovimiento. Un sistema puede pasar espontáneamente de uno a otro estado microscópico y viceversa, sin requerir de estímulos externos. Lo mismo puede decirse en cuanto a los estados macroscópicos; sin embargo, la irreversibilidad del micromovimiento define que, sin coacciones exteriores, el paso espontáneo es solamente unidireccional. Esta paradoja se resuelve, por lo común, considerando que las ecuaciones macroscópicas del movimiento irreversible para cualquier sistema aislado deben ser validas solo a partir de un instante inicial antes del cual el sistema se encontraba aislado. De lo contrario, tratándose de un sistema permanentemente aislado, el proceso irreversible puede ocurrir solamente como resultado de una fluctuación espontánea.

La cuestión más importante en este sentido, es que el crecimiento de la entropía del sistema ocurre solamente a partir de ese instante inicial. Así, para un sistema aislado, las ecuaciones macroscópicas son tales, que para un intervalo infinito de tiempo, todo es reversible, ya que la entropía primero decrece y después crece. Para un sistema que no está aislado siempre, el instante inicial se destaca físicamente y, a partir de él, las ecuaciones macroscópicas solo pueden dar lugar al crecimiento de la entropía, lo que no contradice la irreversibilidad microscópica.

Se requiere definir un conjunto de mecanismos que produzcan un efecto de alteración del estado cuasi estacionario de equilibrio del sistema. Resulta lógico suponer que el desarrollo privilegiado de algunos conductos en detrimento de otros se deba a una combinación entre la cinética del proceso de disolución y el régimen de flujo en el sistema.

Desde el punto de vista de la cinética del proceso de disolución, toda vez que parece claro que la aparición de régimen no lineal de alta velocidad no es, necesariamente, el único mecanismo de triggering, aun cuando la turbulencia contribuya, significativamente, al incremento en la disolución. En tanto indican dos regímenes de disolución: uno fuertemente insaturado y otro próximo a la saturación de calcita. Las tasas de disolución transformadas en tiempos de tránsito bajo ciertas condiciones iniciales de porosidad, gradiente hidráulico y suministro de dióxido de carbono permiten validar la aproximación de la ley cúbica de distribución de velocidad en capilares. La distancia critica de recorrido varia con la tercera potencia del diámetro del conducto.

Asumiendo válidos tales mecanismos, y tomando en cuenta la naturaleza de las reacciones y el control por difusión-dispersión, el problema se reduce a determinar los coeficientes de difusión y de dispersión. Estos últimos, especialmente importantes por cuanto la variación en el régimen de flujo, que de hecho aparece con el incremento del diámetro efectivo, da lugar a que en el nivel fenomenológico macroscópico las propiedades de ambos coeficientes son diferentes.

El tratamiento de la fluctuación termodinámica provocado por la coincidencia entre la aparición del coeficiente fenomenológico de dispersión, el flujo no lineal de alta velocidad, y el incremento de la tasa de disolución, parecen susceptibles de ser tratados como impulsos únicos de duración to en el instante aleatorio, de manera que satisfaga las condiciones en que el instante inicial sea mucho menor que el tiempo total (to << t) y que la probabilidad de los extremales de la funcional sea de magnitud despreciable (to/2t) cuando el tiempo total tiende a infinito.

La producción de entropía en el sistema es el elemento más importante para pronosticar la dirección en que ocurrirán los procesos de desarrollo del cavernamiento., toda vez que, de acuerdo con el segundo principio de la termodinámica, ella ocurre en la dirección del máximo incremento de la entropía.

La función de disipación de energía

Las ecuaciones fenomenológicas son las siguientes para todos los procesos de termofiltración, termoosmosis, dispersión térmica, termodispersión, efecto termoeléctrico directo e inverso, ósmosis química, ultrafiltración, electroósmosis, flujo de corriente, acoplamiento de dispersión, electroforesis y corriente de sedimentación.

- Para el transporte de calor:

- Para el transporte de la fase líquida:

- Para la difusión:

- Para las reacciones químicas:

- Para las Relaciones Recíprocas de Onsager, de manera genérica:

La ecuación básica de disolución de Dreybrodt, tiene la expresión siguiente (Fig. 11):

Fig. 11. Tasa de disolución para películas de agua en régimen turbulento desarrollado considerando la solución en equilibrio con el CO2 atmosférico. La profundidad  del estrato se expresa en cm (según Dreybrodt, 1996)

Donde α es una constante em cm s-1, c es la concentración de calcio en mmol cm3 en la solución y ceq la concetración de equilibrio respecto a la calcia. La función de transferencia de fluido equivale a:

- Para el par poro – grieta:

- Para el par grieta – caverna:

La aplicación del modelo MTDK al Sistema Cavernario Majaguas Cantera mostró los resultados que se presentan en la Fig. 12. La correspondencia con los sectores conocidos / simulados resultó del 89%. En orientación y del 65% en longitud. Esto último se presume que se deba a que quedan sectores en las galerías conocidas aun no descubiertos e integrados al sistema. La longitud de galerías potenciales cuadruplica la sumatoria de las galerías conocidas hasta el momento (33 km; Pedroso et al., 2001; otero et al., 2016).

Fig. 12. Sistema Cavernario Majaguas Cantera (izquierda: real conocido; derecha: simulado)

La conclusión más importante de esta investigación es que la carsificación y el cavernamiento no son procesos aleatorios. Están regidos por leyes termodinámicas y, por tanto, su origen, evolución en el tiempo y posición en el espacio pueden pronosticarse.

3. Conclusiones

La matemática tensorial es fundamental para modelar con precisión la hidrodinámica del karst, un medio inherentemente heterogéneo y anisotrópico. Los enfoques escalares tradicionales son insuficientes y generan errores significativos. La adopción de la matemática tensorial es indispensable para una hidrogeología kárstica predictiva y cuantitativa, superando las limitaciones de los modelos isotrópicos tradicionales y sentando las bases para una gestión sostenible y precisa de estos acuíferos complejos.

La anisotropía hidráulica (tensor K) controla la dirección del flujo. En el karst, el vector de velocidad de Darcy y el gradiente hidráulico no son colineales, haciendo inválida la ortogonalidad clásica entre líneas de flujo y equipotenciales. Esto requiere correcciones tensoriales para mapear correctamente las direcciones de flujo subterráneo. Por su parte, la dispersión de contaminantes es un proceso tensorial (tensor D). La migración de solutos y Líquidos de Fase No Acuosa (NAPL) está gobernada por un tensor de dispersión altamente anisotrópico, donde la dispersión longitudinal en conductos domina sobre la transversal. Modelarla como un escalar subestima gravemente la propagación de la pluma.

Los modelos tradicionales de intrusión marina son inexactos en karst. La aproximación hidrostática (Ghyben-Herzberg) y la suposición de interfaz nítida fallan al no considerar la anisotropía del medio, la miscibilidad de los fluidos y la dinámica real de la zona de mezcla, llevando a pronósticos erróneos en la gestión costera.

La espeleogénesis puede modelarse como un proceso termodinámico direccional. El Modelo Termodinámico del Desarrollo del Karst (MTDK) propone que el cavernamiento no es aleatorio, sino que ocurre en la dirección de máxima producción de entropía, en la intersección de superficies de máxima concentración de flujo y de solvente.

Existe una brecha crítica de conocimiento en coeficientes fenomenológicos. El avance de los modelos está limitado por la falta de datos empíricos para los coeficientes de acoplamiento en las relaciones de no equilibrio (Onsager) que controlan los procesos disolutivos y de transporte a escala de poro y fractura.

Se requiere un cambio de paradigma en la gestión. La aplicación sistemática del formalismo tensorial permite pasar de descripciones cualitativas a predicciones cuantitativas robustas, esenciales para la protección de recursos hídricos, remediación de contaminación y evaluación de riesgos geotécnicos en terrenos kársticos.

Agradecimientos

A Julio J. Valdés Ramos, quien fuera la primera persona de nuestro grupo de trabajo en aplicar estas herramientas al estudio del karst cubano; Pedro J. Astraín, Jefatura Central de Construcciones Militares; Ismael Herrera, Susana Gómez y Martín Díaz, del Instituto de Investigaciones en Matemáticas Aplicadas y en Sistemas de la Universidad Nacional Autónoma de México (IIMAS-UNAM); muy especialmente a Lászlo Kiraly (Centro de Hidrogeología y Geotermia de la Universidad de Neuchatel, Suiza) y Adofo Eraso (Asociación GLACKMA); Denis Louvat, antes en el Organismo Internacional de Energía Atómica; Odon Radai, del Instituto VITUKI (Hungría). Por su extraordinario apoyo y consultas en los trabajos de campo y gabinete a mis colegas, Manuel Acevedo, Raydel Alonso, Yamilé Bustamante, Andrés Díaz, Alberto García, Lázaro Fiallo, Ernesto Flores, Erik García, Néstor Gómez, Mario Guerra, Lorena Guillén, Claudia Hernández, Cecilia March, Manuel Núñez, Vladimir Otero, Yoemí Portuondo, Manuel Rivero, Ernesto Rocamora y Ana, mi compañera.

Bibliografía

Abarca Cameo, E. (2006). Seawater intrusion in complex geological environments [Tesis doctoral, Universitat Politècnica de Catalunya]. https://www.tdx.cat/handle/10803/6762#page=1

Antigüedad, I., Arellano, M., Calvache, M. L., Esteller, M. V., Fagundo, J. R., Gómez, M. A., González, A., González, A., Gutiérrez, J., Llanusa, H., López, M., Molerio, L. F., Morales, T., Morell, I., Mugerza, I., & Pulido, A. (1997). Curso avanzado sobre contaminación de aguas subterráneas: Monitoreo, evaluación, recuperación (Vols. 1-2).

Bear, J. (1972). Dynamics of fluids in porous media. Dover Publications.

Bear, J., & Bachmat, Y. (1990). Introduction to modeling of transport phenomena in porous media. Kluwer Academic Publishers.

Bear, J., & Buchlin, J.-M. (Eds.). (1991). Modelling and applications of transport phenomena in porous media. Springer.

Bear, J., Zaslavsky, D., & Irmay, S. (1968). Physical principles of water percolation and seepage. UNESCO.

Carnahan, C. L. (1976). Non-equilibrium thermodynamics of groundwater flow systems: Symmetry properties of phenomenological coefficients and considerations of hydrodynamic dispersion. Journal of Hydrology, 31(1–2), 125–150. https://doi.org/10.1016/0022-1694(76)90071-3

Chen, Z., Auler, A. S., Bakalowicz, M., Drew, D., Griger, F., Hartmann, J., Jiang, G., Moosdorf, N., Richts, A., Stevanovic, Z., Veni, G., & Goldscheider, N. (2017). The World Karst Aquifer Mapping project: Concept, mapping procedure and map of Europe. Hydrogeology Journal, 25(3), 771–785. https://doi.org/10.1007/s10040-016-1519-3

Cuba, Academia de Ciencias. (1970). Atlas Nacional de Cuba.

de Rooij, R., Graham, W., & Maxwell, R. M. (2013). From a local fracture model to a continuum description of conductivity in fractured rocks. Water Resources Research, 49(11), 7281–7291. https://doi.org/10.1002/2013WR013838

Dreybrodt, W. (1987). The kinetics of calcite dissolution and its consequences to karst evolution from the initial to the mature state. National Speleological Society Bulletin, 49, 31–49.

Dreybrodt, W. (1988). Processes in karst systems: Physics, chemistry, and geology. Springer Series in Physical Environments, Vol. 4. Springer-Verlag.

Dreybrodt, W. (1989). Karst development in its initial state: A model of speleogenesis. En A. Kósa (Ed.), Proceedings of the 10th International Congress of Speleology (Vol. 1, pp. 174–176).

Dreybrodt, W. (1990). The role of dissolution kinetics in the development of karst aquifers in limestone: A model simulation of karst evolution. The Journal of Geology, 98(5), 639–655. https://doi.org/10.1086/629431

Dreybrodt, W. (1992). Dynamics of karstification: A model applied to hydraulic structures in karst terranes. Applied Hydrogeology, 1, 20–32.

Dreybrodt, W. (1993). A model of karstification in the vicinity of hydraulic structures. En G. Gunay, A. I. Johnson, & W. Back (Eds.), Hydrogeological processes in karst terranes (IAHS Publication No. 207, pp. 33–45). International Association of Hydrological Sciences.

Dreybrodt, W. (1995). Principles of early development of karst conduits under natural and man-made conditions revealed by mathematical analysis of numerical models. Water Resources Research, 31(11), 2923–2935. https://doi.org/10.1029/95WR02252

Dreybrodt, W. (1996). Principles of early development of karst conduits under natural and man-made conditions revealed by mathematical analysis of numerical models. Water Resources Research, 32(9), 2923–2935. https://doi.org/10.1029/96WR01332

Dreybrodt, W. (1998). Limestone dissolution rates in karst environments. Bulletin d'Hydrogéologie, 16, 167–183.

Dreybrodt, W., & Buhmann, D. (1991). A mass transfer model for dissolution and precipitation of calcite from solutions in turbulent motion. Chemical Geology, 90(1–2), 107–122. https://doi.org/10.1016/0009-2541(91)90042-G

Dreybrodt, W., Gabrovšek, F., & Romanov, D. (2005). Processes of speleogenesis: A modeling approach. ZRC Publishing, Karst Research Institute.

Dreybrodt, W., & Siemers, J. (1997). Early evolution of karst aquifers in limestone: Models on two-dimensional percolation clusters. En Proceedings of the 12th International Congress of Speleology (Vol. 2, pp. 75–80).

Dreybrodt, W., & Svensson, U. (1989). The effect of MgCO3 on the solubility of calcite: Enhanced aggressiveness. En A. Kósa (Ed.), Proceedings of the 10th International Congress of Speleology (Vol. 1, pp. 177–178).

Egemeier, S. J. (1969). Origin of caves in eastern New York as related to unconfined groundwater flow. National Speleological Society Bulletin, 31(4), 97–111.

Eraso, A. (1982). Consideraciones sobre el problema de la génesis y evolución del Karst.

Eraso, A. (1985/1986). Método de predicción de las direcciones principales de drenaje en el karst. Kobie (Serie Ciencias Naturales), 15, 15–165.

Eraso, A., & Fernández-Rubio, R. (1990). Predicción de las direcciones de drenaje subterráneo en macizos rocosos: Estado actual de la cuestión. En Livro homenagem ao Prof. Carlos Romariz (Vol. 1, pp. 1–20).

Eraso Romero, A., & Domínguez Álvarez, M. d. C. (2000). Método de predicción de las direcciones principales de drenaje subterráneo en macizos anisotrópicos: “Método Eraso.” Boletín Sedeck, 1, 4–10.

Field, M. S., & Pinsky, P. F. (2000). A two-region non-equilibrium model for solute transport in solution conduits. Journal of Hydrology, 229(3–4), 105–116. https://doi.org/10.1016/S0022-1694(00)00159-1

Ford, D. C. (1971). Geologic structure and a new explanation of limestone cavern genesis. Transactions of the Cave Research Group of Great Britain, 13(2), 81–94.

Ford, D. C., & Ewers, R. O. (1978). The development of limestone cave systems in the dimensions of length and depth. Canadian Journal of Earth Sciences, 15(11), 1783–1798. https://doi.org/10.1139/e78-186

Freeze, R. A., & Cherry, J. A. (1979). Groundwater. Prentice-Hall.

Gale, S. J. (1984). The hydraulics of conduit flow in carbonate aquifers. Journal of Hydrology, 70(1–4), 309–327. https://doi.org/10.1016/0022-1694(84)90130-9

Glansdorff, P., & Prigogine, I. (1971). Thermodynamic theory of structure, stability and fluctuations. Wiley-Interscience.

Groenevelt, P. H., & Bolt, G. H. (1969). Non-equilibrium thermodynamics of the soil-water system. Journal of Hydrology, 7(3), 358–388. https://doi.org/10.1016/0022-1694(69)90053-5

Gunn, J. (1986). A conceptual model for conduit flow dominated karst aquifers. En G. Günay & A. I. Johnson (Eds.), Karst water resources (IAHS Publication No. 161, pp. 587–596). International Association of Hydrological Sciences.

Gunn, J. (Ed.). (2004). Encyclopedia of caves and karst science. Fitzroy Dearborn.

Hanna, R. B., & Rajaram, H. (1998). Influence of aperture variability on dissolutional growth of fissures in karst formations. Water Resources Research, 34(11), 2843–2853. https://doi.org/10.1029/98WR02131

Hartmann, A., Goldscheider, N., Wagener, T., Lange, J., & Weiler, M. (2014). Karst water resources in a changing world: Review of hydrological modeling approaches. Reviews of Geophysics, 52(3), 218–242. https://doi.org/10.1002/2013RG000443

Howard, A. D., & Groves, C. G. (1995). Early development of karst systems: 2. Turbulent flow. Water Resources Research, 31(1), 19–26. https://doi.org/10.1029/94WR02139

Hubbert, M. K. (1940). The theory of ground-water motion. The Journal of Geology, 48(8, Pt. 1), 785–944. https://doi.org/10.1086/624930

Hubbert, M. K. (1953). Entrapment of petroleum under hydrodynamic conditions. Bulletin of the American Association of Petroleum Geologists, 37(8), 1954–2026.

Jeannin, P.-Y. (2001). Modeling flow in phreatic and epiphreatic karst conduits in the Hölloch Cave (Muotatal, Switzerland). Water Resources Research, 37(2), 191–200. https://doi.org/10.1029/2000WR900256

Kiraly, L. (1975). Rapport sur l'état actuel des connaissances dans le domaine des caractères physiques des roches karstiques. Eclogae Geologicae Helvetiae, 68(2), 275–278.

Klimchouk, A. B. (2003). Unconfined versus confined speleogenetic settings: Variations of solution porosity. Speleogenesis and Evolution of Karst Aquifers, 1(2), 1–7.

Klimchouk, A. (2005). Conceptualisation of speleogenesis in multi-storey artesian systems: A model of transverse speleogenesis. International Journal of Speleology, 34(1–2), 45–64. https://doi.org/10.5038/1827-806X.34.1.4

Klimchouk, A. B., Ford, D. C., Palmer, A. N., & Dreybrodt, W. (Eds.). (2000). Speleogenesis: Evolution of karst aquifers. National Speleological Society.

Kovacs, A. (2003). Geometry and hydraulic parameters of karst aquifers: A hydrodynamic modeling approach [Tesis doctoral, Université de Neuchâtel]. https://doc.rero.ch/record/4960

Lauritzen, S.-E., Abbott, J., Arnesen, R., Crossley, G., Grepperud, D., Ive, I., & Johnson, S. (1985). Morphology and hydraulics of an active phreatic conduit. Cave Science, 12(3), 139–146.

Liakopoulos, A. C. (1961). Darcy’s law and the field equations of the flow of underground fluids [Tesis doctoral, University of California, Berkeley].

Liakopoulos, A. C. (1965). Darcy’s coefficient of permeability as symmetric tensor of second rank. Hydrological Sciences Journal, 10(3), 41–48. https://doi.org/10.1080/02626666509493401

Long, J. C. S., Remer, J. S., Wilson, C. R., & Witherspoon, P. A. (1982). Porous media equivalents for networks of discontinuous fractures. Water Resources Research, 18(3), 645–658. https://doi.org/10.1029/WR018i003p00645

Maasland, M. (1957). Soil anisotropy and land drainage. En J. N. Luthin (Ed.), Drainage of agricultural lands (pp. 216–285). American Society of Agronomy.

March Delgado, C., & Molerio León, L. F. (1987). A general approach to an algorithm for groundwater flow in karstic aquifers: Hydro- and thermodynamical considerations. En International Symposium on Groundwater Monitoring and Management.

Molerio León, L. F. (1982). Análisis de un modelo teórico de la conductividad hidráulica en el karst. Boletín Grupo Espeleológico Martel de Cuba, 4, 6.

Molerio León, L. F. (1984). El efecto del factor de escala en la interpretación del campo físico de las rocas carsificadas [Resumen]. En XXVII International Geological Congress (Vol. 7, Secc. 16, pp. 468–469).

Molerio León, L. F. (1985a). Dominios de flujo y jerarquización del espacio en acuíferos cársicos. En Simposio XLV Aniversario Sociedad Espeleológica de Cuba (p. 54).

Molerio León, L. F. (1985b). Pronóstico de vías preferenciales de circulación en el carso. En Simposio XLV Aniversario Sociedad Espeleológica de Cuba (pp. 52–53).

Molerio León, L. F. (1985c). El área elemental representativa (AER) para la evaluación de las propiedades físicas del karst: Modelo teórico. En Simposio XLV Aniversario Sociedad Espeleológica de Cuba (p. 45).

Molerio León, L. F. (1986a). Determinación de la conductividad hidráulica direccional en acuíferos cársicos mediante fotointerpretación y cálculo tensorial. En Encuentro Técnico Ramal de Hidroeconomía (p. 66).

Molerio León, L. F. (1986b). Fundamentos hidrodinámicos y termodinámicos para la predicción empírica de la posición, distribución y geometría de las redes de cavernas. En Encuentro Técnico Ramal de Hidroeconomía (pp. 67–68).

Molerio León, L. F. (1986c). El efecto del factor de escala en la interpretación del campo de propiedades físicas de los acuíferos cársicos. Voluntad Hidráulica, 75, 19–26.

Molerio León, L. F. (1988). Particularidades de la transformación matemática de un modelo conceptual del karst [Conferencia invitada]. Taller Internacional de Hidrología Cársica de la Región del Caribe, UNESCO.

Molerio León, L. F. (1989a). Aproximación multidisciplinaria a un modelo matemático del desarrollo del karst. Ingeniería Hidráulica, 10(2), 133–144.

Molerio León, L. F. (1989b). ¿Podemos simular la evolución del karst? En Jornada Científica del Comité Espeleológico Provincial de Villa Clara (p. 23).

Molerio León, L. F. (1989c). El origen de las cavernas. En Jornada Científica del Comité Espeleológico Provincial de Villa Clara (p. 26).

Molerio León, L. F. (1990a). Simulación matemática del desarrollo de las cavernas. En Congreso 50 Aniversario Sociedad Espeleológica de Cuba (p. 55).

Molerio León, L. F. (1990b). CLEO: Presentación de un algoritmo para el cálculo de la estabilidad de las cavernas. En Congreso 50 Aniversario Sociedad Espeleológica de Cuba (p. 89).

Molerio León, L. F. (1992a). Modelo de transporte de masa en la zona no saturada de los acuíferos cársicos: 1/Algoritmo ADRIANA (versión 2.91). En 1er Taller Iberoamericano de Informática y Geociencias (p. 34). Academia de Ciencias de Cuba.

Molerio León, L. F. (1992b). Procesos de transporte de masa en la zona no saturada de los acuíferos cársicos tropicales. En GTICEK. Taller Internacional sobre Cuencas Experimentales en el Karst (pp. 1–15). Universitat Jaume I.

Molerio León, L. F. (1992c). Complementos de un modelo de simulación matemática del desarrollo del karst. En GTICEK. Taller Internacional sobre Cuencas Experimentales en el Karst (pp. 83–92). Universitat Jaume I.

Molerio León, L. F. (1992d). Distribución del cavernamiento en las Sierras del Pesquero, San Carlos, Resolladero y Mesa, Pinar del Río Cuba. En II Congreso Espeleológico de Latinoamérica y el Caribe (pp. 19–20).

Molerio León, L. F. (1993). Problemas ingenieriles en áreas cársicas: La estabilidad de las cavernas. En II Jornadas Venezolanas de Geología Ambiental (p. 15).

Molerio León, L. F. (1997). Transporte de masa, momento y energía. En I. Antigüedad et al., Curso avanzado sobre contaminación de aguas subterráneas (pp.78-??). https://www.researchgate.net/publication/341490688

Molerio León, L. F. (1998). Mathematical simulation of karst development. En International Symposium on Hydrology in the Humid Tropic Environment (pp. 315–325). AIHS.

Molerio León, L. F. (2000). Cien años de teorías sobre el origen de las cavernas [Resumen]. En Congreso Internacional 60 Aniversario Sociedad Espeleológica de Cuba (p. 43). Fundación Antonio Núñez Jiménez.

Molerio León, L. F. (2003a). Modelo del desarrollo de cavernas y conductos cársicos. En V Congreso Cubano de Geología y Minería. Memorias Geomin 2003 (pp. 84–91).

Molerio León, L. F. (2003b). Simulación matemática de los procesos de carsificación, cavernamiento (espeleogénesis) y migración de hidrocarburos en sistemas hipogenéticos. En VI Taller Internacional Informática y Geociencias (p. 6).

Molerio León, L. F. (2005). Prospección hidrogeológica en rocas de baja permeabilidad: Modelo RBP. Mapping: Revista Internacional de Ciencias de la Tierra, 106, 30–38.

Molerio León, L. F. (2006). Efecto de la anisotropía de la transmisividad sobre el campo de flujo en un acuífero cársico litoral. Mapping: Revista Internacional de Ciencias de la Tierra, 108, 77–80.

Molerio León, L. F. (2007). Migración vertical de Líquidos de Fase No Acuosa (LFNA) en acuíferos cársicos. En 2ª Convención Cubana de Ciencias de la Tierra. VII Congreso de Geología (p. 98). ISBN: 978-959-7117-16-2.

Molerio León, L. F. (2013a). Descripción general de la aproximación termodinámica a la simulación del desarrollo de cavernas y conductos cársicos (MTDC). Mundos Subterráneos UMAE, 25, 9–19. http://www.umae.org/paginas/mundos.php

Molerio León, L. F. (2013b). Resumen de la tipología hidrogeológica del karst cubano. Aragonito, 1, 1–25.

Molerio León, L. F. (2015). Disposición final de aguas producidas tratadas de yacimientos gasopetrolíferos carbonatados en acuíferos cársicos litorales someros salinizados. Ciencias de la Tierra y el Espacio, 16(1), 75–87.

Molerio-León, L. F. (2018a). Neotectónica y patrones de cavernamiento en Punta Guanos, Matanzas, Cuba (IV). Estabilidad estructural del macizo. Gota a Gota, 15, 58–63. https://www.resear­chgate.net/publica­tion/325389367_NEOTECTONI­CA_Y_PATRO­NES_DE_CAVERNAMIEN­TO_EN_PUNTA_GUA­NOS_MATAN­ZAS_CUBA_V_TIPIFICACION_INGENIERA_Y_ANALISIS_ESTRUCTURAL_15

Molerio-León, L. F. (2018b). Subsidencia asociada al campo de deformación de galerías subterráneas: El Hundimiento de la Escuela La Inmaculada, Zaruma, Ecuador. Argentina Subterránea, 18(43), 29–51. www.fade.org.ar

Molerio León, L. F. (2019a). Inestabilidad crítica de cavernas (I): Inventario de evidencias. Gota a Gota, 18, 17–31. https://drive.google.com/file/d/12LavS6JAIoQmsLkc7Y9iIA1emwgOl18h/view

Molerio León, L. F. (2019b). Inestabilidad crítica de cavernas (II): Procesos de dolinización. Gota a Gota, 18, 73–79. https://drive.google.com/file/d/1T6mEAxu-W5o5jjEwG16HWN2-3VZbsq4O/view

Molerio León, L. F. (2019c). Modelación matemática analítica de los procesos de colapso de cavernas. Argentina Subterránea, 19(46), 27–44.

Molerio León, L. F. (2020). Inestabilidad crítica de cavernas (III): Procesos de clastificación. Gota a Gota, 19, 10–16.

Molerio-León, L. F. (2021). Petróleo y karst 2: el papel del cavernamiento en el control de la migración de la inyección somera de aguas producidas en medios cársicos. Argentina Subterránea, 21(50), 21–46. https://www.researchgate.net/publication/356893720

Molerio-León, L. F. (2022a). Avances en los algoritmos de base de la tecnología de inyección somera de aguas producidas de yacimientos gasopetrolíferos (TISPM) en medios cársicos. Revista Maya de Geociencias, Edición Especial, 21–36. https://www.researchgate.net/publication/367088104

Molerio-León, L. F. (2022b). El “Efecto de Vecindad” en la superposición de la deformación y el colapso de galerías subterráneas en el karst. Argentina Subterránea, 22(51), 45–58. http://www.fade.org.ar/Bibliografía.html

Molerio-León, L. F. (2023a). Modelo RBP para la prospección de aguas subterráneas en rocas de baja permeabilidad: algoritmo para flujo y transporte en fractura única (single fracture). Geoinformativa, 16(1), 1–14. https://rgi.igp.minem.cu/index.php/rgi

Molerio-León, L. F. (2023b). Tipos de cuevas cubanas y procesos diferenciados de carsificación y cavernamiento. Gota a Gota, 29, 74–88. https://www.researchgate.net/publication/373213056

Molerio León, L. F. (2024a). Resumen conceptual del modelo termodinámico de desarrollo del karst (MTDK). Argentina Subterránea, 24(56), 41–64. https://www.researchgate.net/publication/387498415

Molerio León, L. F. (2024b). Protección estructural de la cueva de Pinagrealan (Norzagaray, Filipinas): Identificación del campo natural de deformaciones elasto-plásticas. Gota a Gota, 32, 31–40. https://www.researchgate.net/publication/384446356

Molerio León, L. F. (2025). Algoritmos acoplados de flujo-transporte con densidad variable para intrusión marina en medios anisotrópicos. Revista Maya de Geociencias, Edición Especial, 66–92. https://www.researchgate.net/publication/395597764

Molerio León, L.F., Raydel Alonso Montero, Ana Margarita Sardiñas (en prensa). Los sistemas cársicos del Polje Jaruco y el problema hidrogeológico de su terminación orienta

Molerio León, L. F., Guerra Oliva, M. G., Flores Valdés, E., & Rocamora Alvarez, E. (1998). Hydrogeological prospecting with geomathematical tools in karstic and fissured non-karstic aquifers. En International Symposium on Hydrology in the Humid Tropic Environment (pp. 337–345). AIHS.

Molerio León, L. F., Hernández Moret, M., Velázquez Sánchez, O., Guerra Oliva, M., & Labrada Cortés, M. (1990). Estudio experimental de resistencia a cargas dinámicas en la bóveda de una caverna. En Congreso 50 Aniversario Sociedad Espeleológica de Cuba (p. 88).

Molerio León, L. F., Portuondo López, Y., & Bustamante Allen, Y. (1996). Migración de hidrocarburos en acuíferos cársicos: Factores de control de transporte y pronóstico del movimiento. Taller Nacional Manejo de Desastres por Derrame de Hidrocarburos.

Montoriol Pous, J., & Assens Caparrós, J. (1957). Sobre el papel desempeñado por el efecto salino en la génesis de ciertas cavidades kársticas desarrolladas en las líneas de costa. Revista de Ciencias, Universidad de Oviedo, VII(1), 81–88.

Mylroie, J. E., & Carew, J. L. (1990). The flank margin model for dissolution cave development in carbonate platforms. Earth Surface Processes and Landforms, 15(5), 413–424. https://doi.org/10.1002/esp.3290150505

Mylroie, J. E., & Carew, J. L. (1995). Karst development on carbonate islands. En D. A. Budd, P. M. Harris, & A. Saller (Eds.), Unconformities and porosity in carbonate strata (AAPG Memoir 63, pp. 55–76). American Association of Petroleum Geologists.

Netushil, A. (1987). Teoría del mando automático: Sistemas no lineales y mando en el caso de acciones aleatorias. Mir.

Núñez Lafitte, M., & Molerio-León, L. F. (2022). Aproximación geomatemática a la prospección de aguas subterráneas en rocas de baja permeabilidad: Modelo RBP. Geoinformativa, 15(1), 17–35. https://rgi.igp.minem.cu/index.php/rgi

Otero Collazo, V., & Molerio-León, L. F. (2022). Condiciones de estabilidad de la dolina y el salón de Los Dibujos de Cinco Cuevas, Boca de Jaruco. Mundo Subterráneo, 7, 42–85. https://sites.google.com/view/mundosubterraneo/geolog%C3%ADa

Oreskes, N., Shrader-Frechette, K., & Belitz, K. (1994). Verification, validation, and confirmation of numerical models in the Earth sciences. Science, 263(5147), 641–646. https://doi.org/10.1126/science.263.5147.641

Otero-Collazo, V. et al. Aspectos geomorfológicos del sistema cavernario Majaguas-Cantera. Campañas de exploración y topografía 2012-2016 (Pinar del Río, Cuba occidental). En: B. Andreo y J. J. Durán (Eds.), El karst y el hombre: las cuevas como Patrimonio Mundial. Nerja (Málaga). Asociación de Cuevas Turísticas Españolas, pp. 407-422

Palmer, A. N. (1991). Origin and morphology of limestone caves. Geological Society of America Bulletin, 103(1), 1–21. https://doi.org/10.1130/0016-7606(1991)103<0001:OAMOLC>2.3.CO;2

Palmer, A. N. (1999). Anisotropy in carbonate aquifers. En A. N. Palmer, M. V. Palmer, & I. D. Sasowski (Eds.), Karst modeling (Karst Waters Institute Special Publication 5, pp. 223–227). Karst Waters Institute.

Palmer, A. N. (2000). Digital modeling of individual solution conduits. En A. B. Klimchouk, D. C. Ford, A. N. Palmer, & W. Dreybrodt (Eds.), Speleogenesis: Evolution of karst aquifers (pp. 194–200). National Speleological Society.

Palmer, A. N., Palmer, M. V., & Sasowsky, I. D. (Eds.). (1999). Karst modeling (Karst Waters Institute Special Publication 5). Karst Waters Institute.

Pankow, J. F., & Cherry, J. A. (1996). Dense chlorinated solvents and other DNAPLs in groundwater: History, behavior, and remediation. Waterloo Press.

Papadopoulos, I. S. (1967). Nonsteady flow to a well in an infinite anisotropic aquifer. En Symposium of Dubrovnik: Hydrology of fractured rocks (Vol. 1, pp. 21–31). UNESCO/AIHS.

Paredes, C., Eraso, A., & Chacón, E. (1993). Aplicación de una metodología de análisis, predicción y modelización de las direcciones principales de drenaje en el karst de Sorbas (Almería). Espeleotemas, 3, 3–16.

Pedroso Herrera, I.I., Pajón Morejón, J.M., Fundora Granda, M.J. y Sagnotti L. 2001. Magnetismo ambiental de sedimentos cuaternarios en el sistema cavernario Majaguas-Cantera. Pinar del Río. IV Congreso de Geología y Minería Geología del Cuaternario, Geomorfología y Carso, 63-73

Polubarinova-Kochina, P. Y. (1962). Theory of ground water movement (J. M. R. de Wiest, Trad.). Princeton University Press. (Trabajo original publicado 1952).

Pool, M., Abarca, E., & Carrera, J. (2007). Simplificaciones en la modelación de la intrusión marina: Validez y alcance. Boletín Geológico y Minero, 118(Núm. Especial), 593–608.

Reimann, T. (2012). Adaptation of numerical modeling approaches for karst aquifer characterization [Tesis doctoral, Technische Universität Dresden].

Reimann, T., & Hill, M. E. (2009). Modelling karst hydrodynamics. En W. B. White (Ed.), Encyclopedia of caves (2ª ed., pp. 506–512). Elsevier/Academic Press.

Reimann, T., Geyer, T., & Birk, S. (2011). Transport and dispersion in karst aquifers: Turbulent effects and tensor representation. Geophysical Research Abstracts, 13, EGU2011-2939.

Renard, P., & de Marsily, G. (1997). Calculating equivalent permeability: A review. Advances in Water Resources, 20(5–6), 253–278. https://doi.org/10.1016/S0309-1708(96)00050-4

Sauter, M., & Liedl, R. (2000). Modelling karst hydrology. En A. B. Klimchouk, D. C. Ford, A. N. Palmer, & W. Dreybrodt (Eds.), Speleogenesis: Evolution of karst aquifers (pp. 291–298). National Speleological Society.

Schoeller, H. (1962). Les eaux souterraines. Masson.

Smart, C. C. (1988). A deductive model of karst landscape evolution based on hydrological probability. Earth Surface Processes and Landforms, 13(3), 271–288. https://doi.org/10.1002/esp.3290130307

Stephansson, O., Hudson, J. A., & Jing, L. (Eds.). (2004). Coupled thermo-hydro-mechanical-chemical processes in geo-systems (Vol. 2). Elsevier.

Szymczak, P., & Ladd, A. J. C. (2006). A network model of channel competition in fracture dissolution. Geophysical Research Letters, 33(5), L05401. https://doi.org/10.1029/2005GL025334

Teutsch, G., & Sauter, M. (1998). Distributed parameter modelling approaches in karst hydrological investigations. Bulletin d'Hydrogéologie, 16, 99–109.

Waltham, T., Bell, F., & Culshaw, M. (2005). Sinkholes and subsidence: Karst and cavernous rocks in engineering and construction. Springer. https://doi.org/10.1007/b138363

Werner Spiller, M. (2004). Physical and numerical experiments of flow and transport in heterogeneous fractured media: Single fracture flow at high Reynolds numbers, and reactive particle transport [Tesis doctoral, Institut National Polytechnique de Toulouse].

White, W. B. (2002). Karst hydrology: Recent developments and open questions. Engineering Geology, 65(2–3), 85–105. https://doi.org/10.1016/S0013-7952(01)00116-8

White, W. B. (Ed.). (2012). Encyclopedia of caves (2ª ed.). Elsevier/Academic Press.

Yevjevich, V. (1971). Sources of stochasticity in geophysical processes. En Proceedings of the Moscow Symposium on Mathematical Models in Geophysics (IAHS Publication No. 116, pp. 83–100). International Association of Hydrological Sciences.

Zhang, Z., Chen, X., Cheng, Q., & Soulsby, C. (2021). Storage dynamics, hydrological connectivity and flux ages in a karst catchment: Conceptual modelling using stable isotopes. Hydrology and Earth System Sciences, 25(2), 901–923. https://doi.org/10.5194/hess-25-901-2021

Sintesis biográfica

Leslie f. Molerio león (ORCID 0000-0001-9667-2258). Cubano. Geólogo de Yacimientos Minerales e Hidrogeólogo-Hidrotécnico, MSc. en Hidrología Isotópica. Ex-Director del Servicio Hidrogeológico de la República de Cuba y del Centro de Hidrología y Calidad de las Aguas (Instituto Nacional de Recursos Hidráulicos), Especialista Principal en Ingeniería (INVERSIONES GAMMA, S.A.), Especialista en Obras Subterráneas y Auditor Ambiental titulado. Experto para agencias de Naciones Unidas, del sistema interamericano y consultoras internacionales. Con 59 años de vida laboral activa, ha participado en unos 490 proyectos de su especialidad en 50 países. Autor o coautor de unos 370 artículos científicos y de 19 libros de su especialidad. Es Miembro de la Academia de Ciencias de Nueva York, Embajador de su Bicentenario y de 34 sociedades científicas nacionales, extranjeras e internacionales.

Conflicto de intereses: El autor de este trabajo declara no presentar conflicto de intereses.

Contribución de autores: Leslie F. Molerio león: Conceptualización, Metodología, Investigación, Análisis formal, Recursos,Curación de datos, Redacción – borrador original, Redacción – revisión y edición, Visualización, Administración del proyecto.

1 Una mala práctica derivada de la traducción “sea water intrusion” asociada a “saline water” confunde intrusión “marina” asociada al intercambio agua dulce-salada del mar en los acuíferos costeros con la intrusión salina que puede estar derivada de otras fuentes de mineralización de las aguas, como la fuga de aguas de capa de yacimientos gasopetrolíferos o la migración desde estructuras diapíricas. La intrusión “marina” es un caso particular de salinización de las aguas terrestres derivados de un intercambio debido exclusivamente entre las aguas de mar y las continentales.