Machine learning and deep learning for numerical simulation of contaminant transport in water systems: a narrative review

1 Introducción

El transporte de contaminantes en sistemas hídricos condiciona la calidad del agua superficial y subterránea, afectando el abastecimiento humano y los ecosistemas asociados. Anticipar la distribución espacial, los tiempos de llegada y la concentración de solutos en ríos, suelos y acuíferos permite evaluar riesgos de exposición y diseñar esquemas de remediación ambiental.

En términos físico-matemáticos, estos fenómenos se describen mediante ecuaciones acopladas de flujo y transporte reactivo. En medios porosos saturados, el marco clásico articula la ley de Darcy con balances de masa que incorporan advección, difusión molecular, dispersión hidrodinámica, sorción y reacciones de degradación [1]. En ríos y canales abiertos, el transporte se acopla a las ecuaciones de flujo en lámina libre o formulaciones unidimensionales de dispersión longitudinal, mientras que en medios fracturados intervienen la geometría de apertura y la conectividad hidráulica del medio.

La resolución práctica de estos problemas se apoya en métodos numéricos consolidados. Esquemas de diferencias finitas, volúmenes finitos y elementos finitos, implementados en códigos estándar como MODFLOW, MT3D-USGS y STOMP, permiten resolver dominios heterogéneos y condiciones de frontera complejas [1,2]. Cada discretización numérica agrega compromisos operativos: difusión numérica en mallas gruesas, dispersión artificial ante gradientes abruptos, restricciones de paso temporal por estabilidad y pérdidas o ganancias artificiales de masa global. Pese a su confiabilidad física, la carga computacional de estos simuladores se vuelve restrictiva en tareas multiconsulta que exigen miles de evaluaciones del modelo, tales como la calibración inversa de parámetros, el análisis de sensibilidad, la propagación de incertidumbre y la optimización de sistemas de remediación [2,3].

Frente a estas restricciones, el aprendizaje automático (machine learning, ML) y el aprendizaje profundo (deep learning, DL) ofrecen alternativas para mitigar la carga de cálculo. En hidrología y calidad del agua, las técnicas de ML y DL se han empleado en predicción de series temporales, mapeo de contaminantes e inferencia espacial [4,5,6]. Sin embargo, la simulación del transporte de solutos no se reduce al ajuste estadístico de curvas observadas: un modelo puede registrar un error cuadrático medio reducido y, aun así, vulnerar la conservación de masa, desplazar frentes advectivos o generar concentraciones negativas sin sustento físico [7]. Por otra parte, cuando un modelo de ML se entrena con datos generados por un simulador numérico, asimila la respuesta del esquema discretizado con sus correspondientes aproximaciones y sesgos de malla, en lugar de la solución física continua del problema.

Por consiguiente, la integración de ML y DL en la mecánica de fluidos ambientales busca complementar las capacidades de los códigos numéricos tradicionales en ámbitos computacionalmente exigentes. La literatura reciente se organiza principalmente en dos grandes vertientes metodológicas: los modelos sustitutos o emuladores, diseñados para aproximar la relación entrada–salida de simuladores costosos para tareas multiconsulta [2,8], y las redes neuronales informadas por física (physics-informed neural networks, PINNs), que incorporan los residuos de las ecuaciones diferenciales y las condiciones de frontera dentro de la función de pérdida para guiar el aprendizaje en regímenes con datos escasos o problemas inversos [9,10,7]. De forma complementaria, emergen formulaciones basadas en redes neuronales de grafos (graph neural networks, GNN) y operadores neuronales, concebidas para capturar conectividades hidrogeológicas no estructuradas y mapeos continuos independientes de malla.

La pregunta central de esta revisión es: ¿en qué condiciones técnicas y bajo qué criterios de consistencia física pueden integrarse los modelos sustitutos y las PINNs en la simulación e inversión del transporte de contaminantes en sistemas hídricos? Para responderla, el presente trabajo examina los supuestos matemáticos, la evidencia empírica, las demandas computacionales y los límites de extrapolación de ambas familias. El análisis parte de que la validez de estos modelos exige evaluar de forma conjunta las métricas estadísticas de error, la conservación de masa, la estabilidad en regímenes dominados por advección, el balance de costo computacional y la interpretabilidad de los parámetros inferidos.

2 Alcance y metodología de revisión

Esta revisión narrativa examina el desarrollo reciente de ML y DL en la simulación y caracterización del transporte de solutos en sistemas hídricos, con énfasis en medios porosos, acuíferos y aguas superficiales gobernados por la ecuación de advección–dispersión–reacción. A diferencia de las revisiones centradas exclusivamente en predicción hidrológica basada en datos observacionales, el análisis se sitúa en la intersección entre la formulación física, la consistencia numérica y el aprendizaje computacional. El corpus bibliográfico se seleccionó a partir de artículos indexados, revisiones temáticas y desarrollos metodológicos con código o formulación reproducible, priorizando estudios que reportan comparaciones directas con simuladores numéricos o datos de laboratorio y campo.

Para estructurar el análisis con claridad, la discusión se concentra en dos familias: (1) modelos sustitutos y emuladores, que emplean redes neuronales, procesos gaussianos y reducción de orden para acelerar la evaluación entrada–salida de simuladores numéricos en regímenes multiconsulta; y (2) PINNs, que integran las ecuaciones gobernantes y condiciones de contorno en la función de pérdida para resolver problemas directos e inversos con datos escasos.

La Figura 1 ilustra el flujo metodológico en el que se insertan estas herramientas. Mientras los sustitutos emulan la respuesta del simulador a partir de conjuntos de datos previamente calculados, las PINNs incorporan los operadores diferenciales directamente en la fase de optimización. En ambas familias, la evaluación técnica debe verificar la reproducibilidad del esquema generador, el costo computacional total y el comportamiento ante variaciones en las condiciones hidrodinámicas del medio.

Figure 1: Integración conceptual de física, datos y ML/DL en el modelado del transporte de contaminantes. Los modelos sustitutos emulan la respuesta entrada–salida de simuladores numéricos para agilizar tareas multiconsulta, mientras que las PINNs incorporan residuos de las ecuaciones gobernantes y condiciones de frontera durante el entrenamiento.

3 Marco físico y numérico del transporte de contaminantes

Para ilustrar el modelo matemático que subyace a las metodologías de ML y DL, considérese el problema prototipo del transporte reactivo de un soluto en un acuífero saturado. En este sistema, la velocidad intersticial promedio de poro ${\textstyle \mathbf {v} }$ se deriva del campo de flujo hidrodinámico según la ley de Darcy:

$\mathbf {v} =-{\frac {\mathbf {K} }{\phi }}\nabla h,$

(1)

donde ${\textstyle \mathbf {K} }$ es el tensor de conductividad hidráulica del medio poroso, ${\textstyle \phi }$ es la porosidad efectiva y ${\textstyle h}$ representa la carga hidráulica [1].

A partir del balance de masa sobre un volumen representativo elemental, el transporte transitorio de la concentración disuelta ${\textstyle C(\mathbf {x} ,t)}$ queda gobernado por la ecuación de advección–dispersión–reacción:

$R_{f}{\frac {\partial C}{\partial t}}+\mathbf {v} \cdot \nabla C=\nabla \cdot (\mathbf {D} \nabla C)-\lambda C+S(\mathbf {x} ,t),$

(2)

donde ${\textstyle R_{f}=1+{\frac {\rho _{b}K_{d}}{\phi }}}$ cuantifica la retardación por sorción lineal (${\textstyle \rho _{b}}$ es la densidad aparente del suelo y ${\textstyle K_{d}}$ el coeficiente de distribución isotérmica), ${\textstyle \mathbf {v} \cdot \nabla C}$ representa el transporte advectivo a lo largo de las líneas de flujo, ${\textstyle \lambda }$ es la tasa de decaimiento reactivo de primer orden y ${\textstyle S(\mathbf {x} ,t)}$ describe la tasa neta de inyección o extracción volumétrica de soluto. En medios anisótropos, el tensor de dispersión hidrodinámica ${\textstyle \mathbf {D} }$ combina la difusión molecular y la dispersión mecánica debida a la variabilidad de velocidades a escala de poro según la formulación de Bear–Scheidegger [1]:

$\mathbf {D} =\left(\alpha _{T}\|\mathbf {v} \|+D_{m}\right)\mathbf {I} +{\frac {\alpha _{L}-\alpha _{T}}{\|\mathbf {v} \|}}\mathbf {v} \mathbf {v} ^{T},$

(3)

donde ${\textstyle \alpha _{L}}$ y ${\textstyle \alpha _{T}}$ son las dispersividades longitudinal y transversal del medio, ${\textstyle D_{m}}$ es el coeficiente de difusión molecular efectivo e ${\textstyle \mathbf {I} }$ es el tensor identidad.

La resolución de la ecuación (2) requiere condiciones iniciales y de contorno: una distribución inicial ${\textstyle C(\mathbf {x} ,0)=C_{0}(\mathbf {x} )}$, concentraciones prescritas en zonas de recarga o entrada (frontera de Dirichlet, ${\textstyle C=C_{b}}$), flujos nulos en límites impermeables (frontera de Neumann homogénea, ${\textstyle \nabla C\cdot \mathbf {n} =0}$) y balances de flujo en puntos de descarga (frontera de Cauchy). En aplicaciones inversas de caracterización ambiental, la fuente ${\textstyle S(\mathbf {x} ,t)}$ o el campo de conductividad ${\textstyle \mathbf {K} }$ suelen ser incógnitas que deben reconstruirse a partir de mediciones dispersas de concentración en pozos de monitoreo [11,12].

La dificultad analítica y numérica de este modelo matemático se cuantifica mediante dos números adimensionales: el número de Péclet (${\textstyle Pe}$) y el número de Damköhler (${\textstyle Da}$), definidos como:

$Pe={\frac {UL}{D}},\qquad Da={\frac {\lambda L}{U}},$

(4)

donde ${\textstyle U}$ es una velocidad característica de flujo, ${\textstyle L}$ es la longitud característica del dominio y ${\textstyle D}$ es el coeficiente de dispersión efectivo. Valores elevados de Péclet (${\textstyle Pe\gg 1}$) identifican regímenes dominados por advección, caracterizados por frentes de pluma pronunciados y gradientes locales agudos. En estas condiciones, los métodos de aproximación numérica y las redes neuronales estándar tienden a sufrir difusión numérica espuria u oscilaciones artificiales [13,14]. Este problema prototipo sirve como base física común para evaluar las dos familias analizadas a continuación: los modelos sustitutos orientados a acelerar simulaciones repetidas de este sistema, y las PINNs diseñadas para resolverlo o estimar sus parámetros.

4 Modelos sustitutos y emuladores

4.1 Aproximación funcional del operador en regímenes multiconsulta

Un modelo sustituto aproxima el mapeo entrada–salida evaluado por un simulador numérico de alto costo. Sea ${\textstyle {\mathcal {G}}:{\mathcal {P}}\to {\mathcal {Y}}}$ el operador de referencia que transforma un vector paramétrico ${\textstyle {\boldsymbol {\mu }}\in {\mathcal {P}}\subset \mathbb {R} ^{p}}$ (conductividades hidráulicas, dispersividades, cinéticas reactivas o concentraciones de frontera) en un espacio de respuestas ${\textstyle \mathbf {y} \in {\mathcal {Y}}}$ (distribuciones de concentración ${\textstyle C(\mathbf {x} ,t)}$ o curvas de ruptura). El emulador aproxima dicha relación mediante una estructura funcional parametrizada por ${\textstyle {\boldsymbol {\eta }}\in \mathbb {R} ^{m}}$:

${\widehat {\mathcal {G}}}_{\boldsymbol {\eta }}:{\mathcal {P}}\longrightarrow {\mathcal {Y}},\qquad {\widehat {\mathcal {G}}}_{\boldsymbol {\eta }}({\boldsymbol {\mu }})\approx {\mathcal {G}}({\boldsymbol {\mu }}).$

(5)

Este planteamiento impone una limitación directa: el sustituto asimila el comportamiento del simulador sobre un conjunto finito de corridas ${\textstyle \{{\boldsymbol {\mu }}_{i}\}_{i=1}^{N_{\mathrm {sim} }}}$, reproduciendo los sesgos de discretización y dispersión del código generador en lugar de las leyes continuas de conservación. La literatura clasifica estos esquemas en tres vertientes principales [2]: modelos basados en datos (redes densas, recurrentes, procesos gaussianos y expansiones de polynomial chaos), métodos de reducción de orden (proper orthogonal decomposition y bases reducidas que proyectan la dinámica sobre variedades de baja dimensión) y enfoques multifidelidad que corrigen simulaciones simplificadas con evaluaciones de alta fidelidad. El diseño experimental, mediante muestreo de hipercubo latino o estrategias adaptativas guiadas por incertidumbre, define la capacidad del modelo para explorar configuraciones físicamente plausibles en acuíferos [2,3].

4.2 Evidencia hidrogeológica y límites de extrapolación

La aplicación de sustitutos se concentra en tareas iterativas de calibración, análisis de incertidumbre, optimización de remediación e identificación de fuentes contaminantes [3]. En acuíferos 2D con advección, dispersión y sorción lineal, redes de perceptrón multicapa (MLP) para concentraciones puntuales y redes de memoria a corto y largo plazo (LSTM) para curvas de ruptura alcanzan ${\textstyle R^{2}\approx 0.99}$ frente a simulaciones de MODFLOW/MT3D [8]. En transporte reactivo multiespecie acoplado (contaminante primario y subproducto), redes entrenadas con diferencias finitas logran aceleraciones inferenciales de decenas a cientos de veces con errores cuadráticos medios (RMSE) inferiores a ${\textstyle 0.03}$, aunque con variabilidad local de precisión según el número de Péclet y desajustes cerca de la fuente [15].

En problemas inversos con condiciones de contorno inciertas, esquemas basados en aprendizaje automático automatizado (AutoML) identifican fuentes puntuales y distribuidas; no obstante, el error relativo crece al pasar de fuentes concentradas a fuentes distribuidas, lo que expone las dificultades de identificabilidad física bajo observaciones incompletas [11]. De manera similar, la emulación en fracturas rugosas acelera la predicción del flujo preferencial, pero experimenta pérdidas de precisión ante zonas de contacto y rugosidades fuera del rango entrenado [16]. Estos resultados confirman que los coeficientes ${\textstyle R^{2}}$ globales pueden encubrir desviaciones locales relevantes en frentes de avance o pozos de captación, y que la extrapolación fuera del espacio paramétrico ${\textstyle {\mathcal {P}}}$ permanece sujeta a fallas imprevistas ante modificaciones en el régimen hidrodinámico.

4.3 Balance computacional y umbral analítico de amortización

La conveniencia operativa de un emulador exige distinguir el costo de preparación fuera de línea (offline) del costo de consulta en línea (online). Resolver ${\textstyle N_{\mathrm {eval} }}$ consultas mediante el código numérico demanda ${\textstyle C_{\mathrm {direct} }(N_{\mathrm {eval} })=N_{\mathrm {eval} }C_{\mathrm {sim} }}$, con ${\textstyle C_{\mathrm {sim} }}$ como tiempo medio por corrida. La vía del sustituto requiere generar ${\textstyle N_{\mathrm {sim} }}$ simulaciones numéricas, un esfuerzo de entrenamiento ${\textstyle C_{\mathrm {train} }}$ e inferencias aceleradas ${\textstyle C_{\mathrm {surr} }}$ por consulta:

$C_{\mathrm {total} }(N_{\mathrm {eval} })=N_{\mathrm {sim} }C_{\mathrm {sim} }+C_{\mathrm {train} }+N_{\mathrm {eval} }C_{\mathrm {surr} }.$

(6)

La condición de rentabilidad ${\textstyle C_{\mathrm {total} }(N_{\mathrm {eval} })<C_{\mathrm {direct} }(N_{\mathrm {eval} })}$ define el umbral analítico de amortización:

$N_{\mathrm {eval} }>N_{\mathrm {eval} }^{*}={\frac {N_{\mathrm {sim} }C_{\mathrm {sim} }+C_{\mathrm {train} }}{C_{\mathrm {sim} }-C_{\mathrm {surr} }}}\approx N_{\mathrm {sim} }+{\frac {C_{\mathrm {train} }}{C_{\mathrm {sim} }}},$

(7)

dado que ${\textstyle C_{\mathrm {surr} }\ll C_{\mathrm {sim} }}$ (con razones habituales ${\textstyle C_{\mathrm {surr} }/C_{\mathrm {sim} }\leq 10^{-3}}$). La relación (7) demuestra que el sustituto incurre en un sobrecosto computacional si el número de evaluaciones no supera el lote de entrenamiento más el equivalente en simulaciones del tiempo de ajuste. Esta condición se cumple en propagación de incertidumbre Monte Carlo o calibración estocástica (${\textstyle N_{\mathrm {eval} }\sim 10^{4}}$ a ${\textstyle 10^{6}}$), pero resulta deficitaria en análisis aislados o cuando modificaciones en la geometría del dominio exigen regenerar la base de datos [2].

5 Aprendizaje informado por física y redes neuronales

5.1 Formulación continua y optimización multiobjetivo

Las PINNs integran las ecuaciones diferenciales gobernantes en la función de pérdida de la optimización. En lugar de limitarse al aprendizaje supervisado sobre pares discretos de entrada–salida, se parametriza la solución continua del campo de transporte mediante una red diferenciable ${\textstyle C_{\boldsymbol {\theta }}(\mathbf {x} ,t)}$, con parámetros ${\textstyle {\boldsymbol {\theta }}\in \mathbb {R} ^{m}}$. El entrenamiento ajusta simultáneamente observaciones, condiciones iniciales, restricciones de frontera y residuos de la ecuación diferencial en el dominio espaciotemporal ${\textstyle \Omega \times (0,T]}$. La función de pérdida se define como:

${\mathcal {L}}({\boldsymbol {\theta }})=\lambda _{d}{\mathcal {L}}_{\mathrm {datos} }({\boldsymbol {\theta }})+\lambda _{b}{\mathcal {L}}_{\mathrm {frontera} }({\boldsymbol {\theta }})+\lambda _{i}{\mathcal {L}}_{\mathrm {inicial} }({\boldsymbol {\theta }})+\lambda _{f}{\mathcal {L}}_{\mathrm {PDE} }({\boldsymbol {\theta }}),$

(8)

donde los coeficientes ${\textstyle \lambda _{d},\lambda _{b},\lambda _{i},\lambda _{f}>0}$ ponderan cada término. Cada componente se calcula sobre conjuntos discretos de puntos:

${\mathcal {L}}_{\mathrm {datos} }({\boldsymbol {\theta }})={\frac {1}{N_{d}}}\sum _{k=1}^{N_{d}}\left|C_{\boldsymbol {\theta }}(\mathbf {x} _{k}^{d},t_{k}^{d})-C_{\mathrm {obs} }^{k}\right|^{2},$ (9) ${\mathcal {L}}_{\mathrm {PDE} }({\boldsymbol {\theta }})={\frac {1}{N_{f}}}\sum _{j=1}^{N_{f}}\left|r_{\boldsymbol {\theta }}(\mathbf {x} _{j}^{f},t_{j}^{f})\right|^{2},$ (10) ${\mathcal {L}}_{\mathrm {inicial} }({\boldsymbol {\theta }})={\frac {1}{N_{i}}}\sum _{i=1}^{N_{i}}\left|C_{\boldsymbol {\theta }}(\mathbf {x} _{i}^{i},0)-C_{0}(\mathbf {x} _{i}^{i})\right|^{2},$ (11)

mientras que ${\textstyle {\mathcal {L}}_{\mathrm {frontera} }({\boldsymbol {\theta }})}$ penaliza la discrepancia en operadores de frontera en ${\textstyle N_{b}}$ puntos sobre ${\textstyle \partial \Omega \times (0,T]}$. Para el transporte de solutos establecido en la Sección 5, el residual continuo ${\textstyle r_{\boldsymbol {\theta }}(\mathbf {x} ,t)}$ corresponde a:

$r_{\boldsymbol {\theta }}(\mathbf {x} ,t)=R_{f}{\frac {\partial C_{\boldsymbol {\theta }}}{\partial t}}+\mathbf {v} \cdot \nabla C_{\boldsymbol {\theta }}-\nabla \cdot \left(\mathbf {D} \nabla C_{\boldsymbol {\theta }}\right)+\lambda C_{\boldsymbol {\theta }}-S(\mathbf {x} ,t),$

(12)

con ${\textstyle R_{f}}$ como factor de retardación, ${\textstyle \mathbf {v} }$ el campo de velocidad intersticial, ${\textstyle \mathbf {D} }$ el tensor de dispersión definido en (3), ${\textstyle \lambda }$ la tasa de decaimiento y ${\textstyle S(\mathbf {x} ,t)}$ el término de inyección o extracción volumétrica de soluto. Las derivadas parciales se evalúan sin mallas mediante diferenciación automática (automatic differentiation, AD) sobre el grafo de la red [9,10].

Durante el entrenamiento surge un compromiso entre el ajuste de datos y la consistencia física: si ${\textstyle \lambda _{f}}$ predomina, la red puede satisfacer la ecuación en promedio pero divergir de las mediciones; si ${\textstyle \lambda _{d}}$ domina, actúa como un interpolador empírico. La adimensionalización de las variables espaciotemporales y los coeficientes hidrogeológicos evita disparidades numéricas entre las magnitudes de los gradientes [10,7]. Habitualmente se emplea el optimizador Adam para una fase exploratoria inicial, seguido de refinamiento local mediante L-BFGS [17,13].

5.2 Acoplamiento hidrodinámico, problemas inversos e incertidumbre

En sistemas subterráneos, el transporte está acoplado al campo de flujo. Las formulaciones multifísicas permiten estimar simultáneamente la conductividad hidráulica ${\textstyle K(\mathbf {x} )}$, la carga hidráulica ${\textstyle h(\mathbf {x} )}$ y la concentración ${\textstyle C(\mathbf {x} ,t)}$ a partir de datos dispersos, superando a las redes puramente basadas en datos [17]. Este esquema abarca problemas directos, inversos y retrospectivos en el tiempo (backward-in-time) gobernados por el flujo de Darcy [18]. En medios heterogéneos, el uso de funciones de activación periódicas (${\textstyle \sin }$) reduce el error frente a ${\textstyle \tanh }$ en uno o dos órdenes de magnitud, facilitando la captura de gradientes abruptos [19].

Para abordar la incertidumbre paramétrica y observacional, las PINNs se han integrado con caos polinomial arbitrario, optimización bayesiana y datos multifidelidad [20], así como esquemas de pre-entrenamiento sobre familias de problemas de transporte en suelos [21]. En problemas inversos, los parámetros hidrogeológicos y las fuentes se tratan como variables optimizables diferenciables. No obstante, hacer entrenable un parámetro no asegura su identificabilidad: la convergencia exige ponderaciones adaptativas y reescalamiento de gradientes para evitar combinaciones espurias de parámetros que reproduzcan el campo de concentración [12]. Estas formulaciones se han aplicado a la identificación de fuentes en canales fluviales 1D [22], advección–difusión 2D [23] y acuíferos mediante esquemas bayesianos en dos etapas [24].

5.3 Regímenes advectivos, balance de masa y criterios de reporte

En regímenes dominados por advección (${\textstyle Pe\gg 1}$), el sesgo espectral de las redes neuronales favorece frecuencias espaciales bajas, lo que tiende a suavizar frentes pronunciados. En problemas hiperbólicos con discontinuidades (como el flujo de Buckley–Leverett), las formulaciones estándar divergen a menos que se introduzcan restricciones de entropía y condiciones admisibles de flujo [13]. Aunque el uso de activaciones adaptativas mitiga oscilaciones espurias en problemas de advección–dispersión–reacción [14], la minimización del residual diferencial (10) sobre puntos discretos no garantiza la conservación global de la masa; la inclusión de una penalización integral de balance reduce esta deficiencia [25].

El reporte riguroso de PINNs en transporte hídrico exige consignar: formulación adimensional, arquitectura y funciones de activación, esquema de muestreo de puntos de colocación, evolución de las ponderaciones de pérdida ${\textstyle \lambda _{k}}$, costo computacional de entrenamiento e inferencia, verificación del balance de masa y evaluación fuera de distribución.

6 Criterios de evaluación y comparación metodológica

La evaluación de modelos de ML y DL en transporte de contaminantes no puede limitarse a métricas estadísticas globales. Coeficientes como el RMSE o el coeficiente de determinación (${\textstyle R^{2}}$) promedian las desviaciones espaciotemporales sobre la totalidad del dominio de cálculo ${\textstyle \Omega }$. Dado que las plumas de soluto ocupan habitualmente una fracción reducida del acuífero, un modelo que prediga concentraciones nulas en zonas no perturbadas puede alcanzar con facilidad valores de ${\textstyle R^{2}>0.98}$, encubriendo subestimaciones severas en frentes de avance o retrasos en la llegada a pozos de abastecimiento.

Para garantizar consistencia física y utilidad hidrogeológica, un protocolo riguroso de validación integra cuatro dimensiones analíticas:

1. Conservación de masa y positividad: Comprueba que el modelo no genere ni destruya masa de soluto de manera artificial. El error relativo instantáneo de balance de masa en el tiempo ${\textstyle t}$ se define formalmente como:

   $\varepsilon _{\mathrm {masa} }(t)={\frac {\left|\int _{\Omega }\phi (\mathbf {x} )R_{f}(\mathbf {x} )C_{\boldsymbol {\theta }}(\mathbf {x} ,t)\,d\Omega -M_{\mathrm {teorico} }(t)\right|}{M_{0}+\int _{0}^{t}{\dot {M}}_{\mathrm {in} }(\tau )\,d\tau }},$

   (13)

donde ${\textstyle M_{0}=\int _{\Omega }\phi R_{f}C_{0}\,d\Omega }$ es la masa inicial disuelta y adsorbida en el medio poroso, ${\textstyle {\dot {M}}_{\mathrm {in} }}$ denota el flujo másico inyectado, y la masa teórica acumulada ${\textstyle M_{\mathrm {teorico} }(t)}$ integra los flujos advectivos y dispersivos en la frontera junto con la tasa de decaimiento reactivo:

$M_{\mathrm {teorico} }(t)=M_{0}+\int _{0}^{t}\left[{\dot {M}}_{\mathrm {in} }(\tau )-{\dot {M}}_{\mathrm {out} }(\tau )-\int _{\Omega }\phi (\mathbf {x} )\lambda C_{\boldsymbol {\theta }}(\mathbf {x} ,\tau )\,d\Omega \right]d\tau {.}$

(14)

A la par de (13), debe verificarse la condición de positividad ${\textstyle C_{\boldsymbol {\theta }}(\mathbf {x} ,t)\geq 0}$, cuya vulneración genera artefactos numéricos inadmisibles en ingeniería ambiental.

2. Respuesta ante regímenes advectivos ($Pe$): El comportamiento predictivo debe examinarse en función del número de Péclet (${\textstyle Pe=UL/D}$). En regímenes con ${\textstyle Pe\gg 1}$, se evalúa si la arquitectura amortigua el frente contaminante mediante difusión artificial o si induce oscilaciones espurias aguas arriba del frente de avance.
3. Amortización del costo computacional: Evaluación formal del balance entre el costo preparatorio ${\textstyle C_{\mathrm {offline} }}$ (generación de lotes de simulación y convergencia de gradientes) y el costo inferencial ${\textstyle C_{\mathrm {online} }}$, aplicando el umbral analítico de consultas ${\textstyle N_{\mathrm {eval} }^{*}}$ deducido en (7).
4. Generalización y extrapolación fuera de distribución: Verificación del desempeño ante escenarios no contemplados en la etapa de optimización, tales como variaciones en el gradiente regional de flujo, desplazamientos imprevistos de la fuente o campos de conductividad hidráulica con estructuras de heterogeneidad no muestreadas.

La Tabla 1 sintetiza las diferencias operativas y conceptuales entre las dos familias metodológicas analizadas.

Table. 1 Comparación metodológica entre modelos sustitutos y PINNs en el transporte de contaminantes.

Dimensión	Modelos sustitutos y emuladores	PINNs
Entrada y salida	Parámetros globales o de diseño ${\textstyle \to }$ campos de concentración o series en pozos.	Coordenadas ${\textstyle (\mathbf {x} ,t)}$ y parámetros ${\textstyle \to }$ concentración local ${\textstyle C(\mathbf {x} ,t)}$ y propiedades del medio.
Requerimiento de datos	Cientos a miles de simulaciones numéricas previas mediante diseño de experimentos.	Mediciones dispersas de concentración y puntos de colocación en el dominio; no requieren simulador previo.
Tratamiento de la física	Indirecto; la física proviene del simulador numérico que genera la base de entrenamiento.	Directo; el residuo diferencial de advección–dispersión–reacción se minimiza en la función de pérdida.
Costo computacional	Alto costo inicial (${\textstyle C_{\mathrm {offline} }}$); inferencia en milisegundos (${\textstyle C_{\mathrm {online} }}$).	Optimización prolongada para cada escenario; cada nuevo dominio o condición exige reentrenamiento.
Conservación y extrapolación	Conservación aproximada; degradación rápida fuera del espacio paramétrico entrenado ${\textstyle {\mathcal {P}}}$.	Susceptible a difusión numérica y pérdidas de masa si ${\textstyle Pe\gg 1}$; regularización continua en interpolación.
Uso principal	Análisis de sensibilidad global, optimización de pozos y propagación de incertidumbre Monte Carlo.	Problemas inversos, estimación de parámetros hidrogeológicos y localización de fuentes contaminantes.

7 Desafíos abiertos para una adopción confiable

La adopción de estas técnicas en la toma de decisiones hidrogeológicas enfrenta retos metodológicos específicos que condicionan su fiabilidad práctica en escenarios de ingeniería ambiental.

En primer término, persiste una falta de problemas de referencia (benchmarks) compartidos en hidrogeología computacional. La literatura evalúa las arquitecturas sobre geometrías ad hoc, fuentes simplificadas o condiciones de flujo homogéneo. Esta disparidad en los casos de prueba impide comparar de forma concluyente la precisión, la demanda de memoria y el tiempo de entrenamiento entre distintas alternativas. La comunidad requiere bancos de prueba abiertos basados en configuraciones representativas (como medios porosos con heterogeneidad espacial correlacionada o sistemas fracturados), provistos de mallas y particiones de datos estandarizadas.

En segundo término, existe una brecha notable entre los datos sintéticos de simulaciones numéricas y las condiciones del monitoreo ambiental de campo. Mientras que los simuladores numéricos entregan mallas densas y regulares, los acuíferos reales se registran mediante redes limitadas de pozos con muestreos asincrónicos, perturbados por ruido instrumental e incertidumbre en la estratigrafía subsuperficial. El desarrollo de modelos debe priorizar la asimilación de datos escasos y dispersos, contemplando las incertidumbres asociadas a las condiciones de frontera hidrogeológicas.

En tercer lugar, la conservación estricta de masa en regímenes dominados por advección continúa siendo un problema abierto. En formulaciones estándar de PINNs, el residuo diferencial se penaliza como una pérdida suave en la optimización; en presencia de frentes de soluto pronunciados (${\textstyle Pe\gg 1}$), esto propicia difusión artificial o variaciones no físicas en la masa total de contaminante [25]. Para mitigar este defecto, se investigan esquemas con conservación intrínseca, tales como redes neuronales de volúmenes finitos (FINN) [26] o restricciones que imponen el balance de flujos mediante multiplicadores de Lagrange.

Asimismo, la cuantificación de incertidumbre en problemas inversos exige mayor rigor. La identificación de fuentes de contaminación o la estimación de conductividades hidráulicas presentan una acusada no unicidad: diversas combinaciones de parámetros generan curvas de concentración equivalentes en los pozos de observación [23,12]. Las predicciones puntuales resultan insuficientes para la gestión de riesgos ambientales. La integración de redes neuronales bayesianas, ensambles y esquemas de remuestreo (bootstrapping) acoplados a modelos con estructura física [27,24] proporciona una vía para estimar intervalos de confianza reproducibles.

Finalmente, surgen arquitecturas complementarias dirigidas a resolver limitaciones topológicas y funcionales. Las GNN permiten modelar la conectividad hidráulica en redes de monitoreo con geometría irregular [28]; los operadores neuronales, como DeepONet y las redes de operadores de Fourier (FNO), abordan la aproximación de operadores no lineales entre espacios de funciones continuas, con capacidad para predecir familias completas de soluciones parametrizadas por condiciones de frontera, campos de conductividad o tasas de inyección variables con independencia de malla [29,30,31,32], incluyendo extensiones hacia flujos reactivos multiespecie en medios heterogéneos [33]; y los modelos basados en transformers exploran la captura de dependencias de largo alcance en series hidrogeológicas complejas [34]. Su consolidación operativa dependerá de verificar que preservan los principios de balance físico y garantizan estabilidad numérica con una inversión computacional competitiva frente a los simuladores consolidados.

8 Conclusiones

El análisis bibliográfico desarrollado en este trabajo permite formular tres conclusiones principales sobre la función y los límites de ML y DL en el transporte de contaminantes hídricos:

En primer lugar, los métodos analizados complementan a los simuladores numéricos tradicionales sin reemplazarlos. Los modelos sustitutos justifican su inversión de desarrollo en problemas de optimización de pozos, diseño de remediación y análisis de incertidumbre donde se evalúan miles de escenarios; en contraste, las PINNs aportan ventajas en la estimación de parámetros hidrogeológicos y la localización de fuentes a partir de datos dispersos, al integrar las ecuaciones gobernantes sin requerir una discretización espacial rígida.

En segundo término, la precisión y la estabilidad física de estas herramientas están sujetas al régimen de transporte del sistema. En regímenes advectivos caracterizados por valores elevados del número de Péclet (${\textstyle Pe\gg 1}$), así como en acuíferos con heterogeneidades marcadas, los métodos sufren pérdidas no físicas de masa o difusión espuria del frente contaminante si se limitan a funciones de pérdida estándar. La inclusión de formulaciones de balance estricto y discretizaciones de tipo volumen de control permite evitar predicciones hidrogeológicamente incompatibles con los principios de conservación.

Finalmente, la consolidación metodológica en este ámbito demanda superar las evaluaciones limitadas a métricas estadísticas globales como el coeficiente de determinación o el error cuadrático medio. La utilidad real de estos modelos en la protección de recursos hídricos exige someterlos a bancos de pruebas estandarizados de la comunidad hidrogeológica, reportar balances de masa, evaluar su comportamiento fuera del dominio de entrenamiento y cuantificar sistemáticamente la incertidumbre en los parámetros inferidos ante el ruido y la escasez propios del monitoreo ambiental.

Agradecimientos

Los autores agradecen a la Secretaría de Ciencias, Humanidades, Tecnología e Innovación (SeCiHTI) de México, a la Coordinación de la Investigación Científica de la Universidad Michoacana de San Nicolás de Hidalgo (CIC-UMSNH), a SIIIA MATH: Soluciones en Ingeniería y al Aula CIMNE-Morelia por el financiamiento brindado para el desarrollo, la revisión bibliográfica y la redacción del presente artículo.

Asimismo, expresan su agradecimiento a los editores y revisores del Boletín, cuyos comentarios y sugerencias contribuyeron a consolidar la estructura y precisión técnica de este artículo.

BIBLIOGRAPHY

[1] Bear, Jacob and Cheng, Alexander H.-D. (2010) "Modeling Groundwater Flow and Contaminant Transport". Springer

[2] Asher, M. J. and Croke, B. F. W. and Jakeman, A. J. and Peeters, L. J. M. (2015) "A review of surrogate models and their application to groundwater modeling", Volume 51. Water Resources Research 5957–5973

[3] Luo, Jiannan and Ma, Xi and Ji, Yefei and Li, Xueli and Song, Zhuo and Lu, Wenxi. (2023) "Review of machine learning-based surrogate models of groundwater contaminant modeling", Volume 238. Environmental Research 117268

[4] Shen, Chaopeng. (2018) "A Transdisciplinary Review of Deep Learning Research and Its Relevance for Water Resources Scientists", Volume 54. Water Resources Research 8558–8593

[5] Haggerty, Ryan and Sun, Jianxin and Yu, Hongfeng and Li, Yusong. (2023) "Application of machine learning in groundwater quality modeling - A comprehensive review", Volume 233. Water Research 119745

[6] Zhi, Wei and Appling, Alison P. and Golden, Heather E. and Podgorski, Joel and Li, Li. (2024) "Deep learning for water quality", Volume 2. Nature Water 228–241

[7] Ma, Qingshan and Gong, Qixin and Ge, Weiya and Jing, Miao and Guo, Fei and Butscher, Christoph and Taherdangkoo, Reza. (2026) "Physics-informed neural networks for groundwater: evidence, limits, and a roadmap", Volume 85. Environmental Earth Sciences 123

[8] Yu, Xiayang and Cui, Tao and Sreekanth, J. and Mangeon, Stephane and Doble, Rebecca and Xin, Pei and Rassam, David and Gilfedder, Mat. (2020) "Deep learning emulators for groundwater contaminant transport modelling", Volume 590. Journal of Hydrology 125351

[9] Raissi, M. and Perdikaris, P. and Karniadakis, G. E. (2019) "Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations", Volume 378. Journal of Computational Physics 686–707

[10] Karniadakis, George Em and Kevrekidis, Ioannis G. and Lu, Lu and Perdikaris, Paris and Wang, Sifan and Yang, Liu. (2021) "Physics-informed machine learning", Volume 3. Nature Reviews Physics 422–440

[11] Xu, Yaning and Lu, Wenxi and Pan, Zidong and Luo, Chengming and Bai, Yukun and Qiu, Shuwei. (2024) "Groundwater contaminant source identification considering unknown boundary condition based on an automated machine learning surrogate", Volume 15. Geoscience Frontiers 101732

[12] Berardi, Marco and Difonzo, Fabio V. and Icardi, Matteo. (2025) "Inverse Physics-Informed Neural Networks for transport models in porous materials", Volume 435. Computer Methods in Applied Mechanics and Engineering 117628

[13] Fraces, Cedric G. and Tchelepi, Hamdi. (2021) "Physics Informed Deep Learning for Flow and Transport in Porous Media". SPE Reservoir Simulation Conference

[14] Hou, Qingzhi and Du, Honghan and Sun, Zewei and Wang, Jianping and Wang, Xiaojing and Wei, Jianguo. (2023) "PINN-CDR: A Neural Network-Based Simulation Tool for Convection-Diffusion-Reaction Systems", Volume 2023. International Journal of Intelligent Systems

[15] Nguyen, Thu-Uyen and Suk, Heejun and Liang, Ching-Ping and Ho, Yu-Chieh and Chen, Jui-Sheng. (2025) "Using Machine Learning to Develop a Surrogate Model for Simulating Multispecies Contaminant Transport in Groundwater", Volume 12. Hydrology 185

[16] Xue, Haoyu and Xu, Wenjie and Lou, Yihuai and Hu, Yingtao and Gao, Wu and Li, Jinlong and Zhuang, Duanyang and Liu, Wei and Zhan, Liangtong and Chen, Yunmin. (2026) "Deep-learning-based surrogate model for solute transport in rough-walled fractures", Volume 191. Computers and Geotechnics 107813

[17] He, QiZhi and Barajas-Solano, David and Tartakovsky, Guzel and Tartakovsky, Alexandre M. (2020) "Physics-informed neural networks for multiphysics data assimilation with application to subsurface transport", Volume 141. Advances in Water Resources 103610

[18] He, QiZhi and Tartakovsky, Alexandre M. (2021) "Physics-Informed Neural Network Method for Forward and Backward Advection-Dispersion Equations", Volume 57. Water Resources Research

[19] Faroughi, Salah A. and Soltanmohammadi, Ramin and Datta, Pingki and Mahjour, Seyed Kourosh and Faroughi, Shirko. (2024) "Physics-Informed Neural Networks with Periodic Activation Functions for Solute Transport in Heterogeneous Porous Media", Volume 12. Mathematics 1 63

[20] Chen, Xiaoli and Duan, Jinqiao and Karniadakis, George Em. (2020) "Learning and meta-learning of stochastic advection-diffusion-reaction systems from sparse measurements", Volume 32. European Journal of Applied Mathematics 397–420

[21] Ke, Ze-Wei and Wei, Sheng-Jie and Yao, Shi-Yuan and Chen, Si and Chen, Yun-Min and Li, Yu-Chao. (2025) "Pre-trained Physics-Informed Neural Networks for Analysis of Contaminant Transport in Soils", Volume 180. Computers and Geotechnics 107055

[22] Zhao, Xu and Liu, Haifei and Leng, Fei and Yang, Wei and Yin, Xinan and Yi, Yujun. (2026) "Source Tracing for Pollutants in River Channels Based on a Physics-Informed Neural Network", Volume 62. Water Resources Research

[23] Chuprov, Ivan and Derkach, Denis and Efremenko, Dmitry and Kychkin, Aleksei. (2025) "Application of Physics-Informed Neural Networks for Solving the Inverse Advection-Diffusion Problem to Localize Pollution Sources"

[24] Zhou, Zhengchen and Yan, Xuesong and Hu, Chengyu. (2026) "Two-Stage Bayesian Physics-Informed Neural Networks for Groundwater Contaminant Source Identification", Volume 152. Journal of Water Resources Planning and Management

[25] Sigalingging, Eszra Forenita and Liu, Kin Men and Setianto, Setianto and Faizal, Ferry. (2026) "Mass-Conserving Physics-Informed Neural Networks For The One-Dimensional Advection-Diffusion Equation"

[26] Praditia, Timothy and Karlbauer, Matthias and Otte, Sebastian and Oladyshkin, Sergey and Butz, Martin V. and Nowak, Wolfgang. (2022) "Learning Groundwater Contaminant Diffusion-Sorption Processes With a Finite Volume Neural Network", Volume 58. Water Resources Research

[27] Scheurer, Stefania and Frenner, Riccardo and Bruennette, Tim and Oladyshkin, Sergey and Nowak, Wolfgang. (2026) "Efficient confidence interval computation for physics-aware machine learning of diffusion-sorption models", Volume 8. Frontiers in Water 1813791

[28] Pang, Min and Du, Erhu and Zheng, Chunmiao. (2024) "Contaminant Transport Modeling and Source Attribution With Attention-Based Graph Neural Network", Volume 60. Water Resources Research

[29] Lu, Lu and Jin, Pengzhan and Pang, Guofei and Zhang, Zhongqiang and Karniadakis, George Em. (2021) "Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators", Volume 3. Nature Machine Intelligence 218–229

[30] Li, Zongyi and Kovachki, Nikola and Azizzadenesheli, Kamyar and Liu, Burigede and Bhattacharya, Kaushik and Stuart, Andrew and Anandkumar, Anima. (2020) "Fourier Neural Operator for Parametric Partial Differential Equations"

[31] Kovachki, Nikola and Li, Zongyi and Liu, Burigede and Azizzadenesheli, Kamyar and Bhattacharya, Kaushik and Stuart, Andrew and Anandkumar, Anima. (2023) "Neural Operator: Learning Maps Between Function Spaces With Applications to PDEs", Volume 24. Journal of Machine Learning Research 89 1–97

[32] Azizzadenesheli, Kamyar and Kovachki, Nikola and Li, Zongyi and Liu-Schiaffini, Miguel and Kossaifi, Jean and Anandkumar, Anima. (2024) "Neural operators for accelerating scientific simulations and design", Volume 6. Nature Reviews Physics 320–328

[33] Karimi, Mina and Bhattacharya, Kaushik. (2024) "A Learning-Based Multiscale Model for Reactive Flow in Porous Media", Volume 60. Water Resources Research e2023WR036303

[34] Bai, Tao and Tahmasebi, Pejman. (2022) "Characterization of groundwater contamination: A transformer-based deep learning model", Volume 164. Advances in Water Resources 104217