1. Introducción
A lo largo del tiempo, el Planeta se ha visto atacado por una serie de amenazas, como los contaminantes del suelo, que lo han llevado a su acelerado deterioro. Por ejemplo, la contaminación de suelos agrícolas debido a infiltraciones de agua contaminada, entre otros. Ante ello, surgen interrogantes como: ¿De qué manera se contaminan los suelos?, ¿Cuáles son los niveles de contaminación de los suelos a 10, 15 y 20 años?, ¿Qué medidas deberían adoptarse para evitar consecuencias negativas futuras? ¿Cuáles son los niveles de efectividad al aplicar medidas de prevención?
La búsqueda de respuestas conlleva al análisis de problemas producidos por infiltraciones de agentes contaminantes en el suelo. (Mulligan et al, 2004). La problemática planteada exige tomar medidas, una de ellas es el desarrollo de modelos matemáticos, modelos numéricos y programas computacionales que caractericen y simulen el fenómeno, analicen posibles escenarios y determinen cómo actúa y se desenvuelve el contaminante en el medio poroso. Todo ello en favor de garantizar la preservación de los recursos naturales y la vida del planeta (Gonzales de Vallejo, 2002).
El objetivo de este trabajo es desarrollar un modelo matemático que represente el fenómeno de infiltración de contaminante en el suelo basado en la ecuación de Richards no lineal, su modelo numérico del método de los elementos finitos y su implementación en un programa computacional que permita realizar simulaciones en diferentes escenarios y facilite una mejor comprensión del proceso de infiltración. Para que sirva como una herramienta para la toma de decisiones en beneficio de nuestro mundo.
La Ley de Darcy
Darcy (1856), demostró que el caudal 
que atraviesa el permeámetro era proporcional a la sección A y a la diferencia de alturas 
dividido por el diferencial de longitud 
es decir:
Donde k es la constante de proporcionalidad. Además determinó que al utilizar otro tipo de arena (más gruesa o más fina, o a su vez la mezcla de estas) y combinando todas las variables, se volvía a cumplir la ecuación anterior, pero la constante de proporcionalidad lineal era distinta. Por tanto esta, constante era una característica propia del material y la llamo permeabilidad. Dado que el caudal esta en 
, la sección en 
, y 
y son longitudes, se tiene que las unidades de permeabilidad 
son las de la velocidad 
. Por otra parte tomando incrementos infinitesimales se deduce que, la Ley de Darcy se expresa de esta forma:
En donde:
es el cambio de altura respecto a la horizontal
caudal que circula por de sección
conductividad hidráulica
gradiente hidráulico expresado en incrementos infinitesimales
El signo menos se debe a que el caudal es una magnitud vectorial, cuya dirección es hacia los decrecientes; es decir, que o 
es negativo y, por tanto, el caudal será positivo.
Ley de Darcy-Buckingham
La Ley de Darcy se ha obtenido para un solo fluido, sin embargo los suelos no saturados en sus poros, presentan dos fluidos (agua y aire), lo cual conlleva a la no aplicabilidad de este principio para estos casos. El efecto de taponamiento de los poros en los suelos, son causados por las burbujas de aire, que retienen el líquido impidiendo su natural permeabilidad, esta es la razón principal para que un suelo parcialmente saturado tenga menor permeabilidad que otro totalmente saturado.
Buckingham (1907), propuso que el grado de saturación va en aumento en la misma medida en que las burbujas de aire disminuyen en el tiempo ya que son arrastradas por las corrientes de agua, lo cual permite concluir que la permeabilidad de un suelo parcialmente saturado aumenta significativamente con el paso del tiempo.
El incremento de la constante de proporcionalidad (permeabilidad) en suelos semi - saturados es mayor mientras la presión del líquido aumenta, ya que esto provoca una considerable disminución del volumen que ocupan las burbujas de aire. Luego en este tipo de suelo la permeabilidad K depende de h y la ley de Darcy se corrige de la siguiente forma:
Ecuación de la continuidad
La densidad de un flujo (líquido o gas) homogéneo puede depender de mucho factores tales como la temperatura y la presión a la que está sometido. Para lo líquidos, la densidad varía muy poco dentro de amplios límites de presión y temperatura.
Con fines prácticos en la contaminación ambiental, los líquidos que no están sujetos a una fuente contaminante serán considerados de densidad constante. Por el contrario, si el líquido está sujeto a una fuente contaminante cambiará de densidad. La densidad de los gases es muy sensible a los cambios de temperatura y presión, la densidad se define como: 
o 
siendo sus unidades de denominación 
o 
.
Sea 
un volumen de control que se mueve con una velocidad de 
. La densidad del fluido contenido en es 
, de donde 
e integrando sobre , tenemos:
Derivando con respecto al tiempo, nos da:
y por el teorema de transporte de Reynolds, obtenemos, (Richards,1931)
Sistema conservativo
Un sistema se dice conservativo si la masa del sistema no se crea, ni se destruye; es decir que la masa permanece constante durante todo el proceso, esto es 
para todo 
es decir,
Como la integral precedente es válida para cualquier volumen de control, se demuestra que
Ecuación de Richards
La ecuación de Richards se obtiene de reemplazar la ley de Darcy-Buckingham (3) en la ecuación de conservación de masa o ecuación de continuidad para un sistema conservativo 
sobre 
, (Darcy, 1856),
Así la ecuación quedaría de la forma:
En el presente trabajo a fin de darle mayor generalidad a nuestro estudio se considerará la siguiente ecuación del tipo parabólico no lineal:
Sean 
abiertos, acotado, de frontera 
lipschisiana a trozos, 
. Ponemos , 
, 
.
Consideraremos el problema siguiente:
Formulación débil
Según Sayas (2007), para la obtención de una formulación débil, procedemos del modo siguiente, multipliquemos a la ecuación diferencial del problema (10) por una función 
tal que 
, e integremos sobre 
, dando como resultado:
Por el teorema de derivación bajo el signo de integración, (Richards, 1931) se tiene:
Ahora por el teorema de Green -Gauss, que se puede encontrar en (Brezis, 1983);
Y debido a que , se tiene:
Reemplazando estos resultados se obtiene la siguiente ecuación, para todo 
tal que .
Para cada 
, se define:
Con estas notaciones la formulación débil del problema parabólico se escribe:
Modelo numérico y computacional
Sea 
un subespacio de dimensión finita 
una base de 
La formulación aproximada se establece del modo siguiente: para cada 
, hallar 
solución de:
Sean 
y 
La derivada 
en 
se aproxima por diferencias finitas progresivas del modo siguiente:
Para el resto de la ecuación usamos el método theta con 
. Luego el problema consiste en hallar 
solución de
.
Para 
(el método se llama de Crank Nicholson), se obtiene
, de donde
Definimos:
, 
,
Se tiene el siguiente problema:
Hallar solución de : 
La derivada de Fréchet de 
en 
viene dada por:
Aplicando el método de Newton se sigue:
y 
tal que 
El problema discreto se expresa del modo siguiente:
De forma más general el algoritmo a implementar seria:
Algoritmo numérico
ENTRADAS:
y N (Topología de la malla)
(Solución inicial)
(Solución inicial arbitraria)
(Tiempo de simulación)
(Número de divisiones para la discretización temporal)
(Número máximo de iteraciones)
(Precisión con la que deseamos la solución)
PASOS:
PARA 
REPETIR
Calcular la matriz 
con información de 
y 
Calcular el vector 
con información de y
Resolver el sistema 
HASTA QUE 
o 
FIN DEL PARA
SALIDA: Solución del problema 
Pruebas numéricas
En esta sección se va a tomar un ejemplo particular de la ecuación de Richards. Es decir, se fijan todas las funciones involucradas en el problema. De tal manera que la función 
cumpla con las condiciones de frontera y condiciones iniciales y se reemplazan en la ecuación para obtener la función 
,con lo cual se ha construido un ejemplo particular de ecuación de Richards cuya solución es conocida, lo que nos permite validar el programa computacional y se obtiene algunos resultados numéricos.
Adicionalmente, se deja constante el tiempo y variamos el número de elementos (triángulos) de la discretización espacial, se calcula el error con la norma de 
para lo cual se integran en valor absoluto la diferencia de la solución exacta y aproximada sobre cada triángulo y luego se suman., Con lo cual se obtiene una muy buena aproximación de esta norma. Así mismo se registra el tiempo que se demora la computadora en realizar los cálculos, y finalmente se registran los resultados, como lo muestra la Tabla 1.
Tabla 1
Resultados numéricos
Como se puede apreciar, mientras aumentan el número de elementos en el mallado, el error disminuye de manera considerable, aunque el tiempo crece de manera exponencial, Sin embargo, este inconveniente se puede solucionar según las características del computador, en este caso la simulación se realizó en una laptop core i3.
Como segunda prueba, se analiza la convergencia del método de Newton en la solución de la ecuación de Richards y se evalúa la diferencia entre la solución numérica en el k-ésimo paso y en el k-1 ésimo paso, los resultados se verifican en la Tabla 2.
Tabla 2
Convergencia del método de Newton
2. Resultados
En esta sección se encontrarán algunos resultados importantes sobre todo del simulador, corridas del programa con una comparación entre las soluciones exactas versus las aproximadas así como también graficas de en 2D y 3D del proceso de infiltración.
Se presentan los resultados del programa computacional que resuelve la ecuación de Richards bidimensional, se exponen las soluciones aproximadas y exactas en los nodos de la malla para lo cual se ha usado un problema cuya solución se conoce.
En la Figura 1 se presentan los resultados de la solución exacta y la aproximada; donde se nota que los resultados son muy similares, así validamos tanto el modelo numérico y su implementación.
El simulador tiene dos maneras de visualizar los resultados, la primera, en la Figura 2, muestra el resultado en tres dimensiones es decir, usa la malla para graficar la función resultante. Lógicamente como es un problema evolutivo dicha función va cambiando en el tiempo, lo cual nos una buena idea de cómo evoluciona el fenómeno.
Fig. 1:
Comparación de resultados de la solución exacta y aproximada
Fig. 2:
Proceso de Infiltración en 3D
La segunda forma de visualización es en dos dimensiones es decir a cada posición (x, y) del plano le corresponde un código de colores que varía de amarillo que significa 0 contaminación, hasta rojo intenso, que es donde la función toma los valores más altos. Lo que implica una mayor contaminación. En el caso del fenómeno de infiltración de agua en el suelo se puede observar en un corte transversal como el fluido va poco a poco ingresando en el suelo. En la Figura 3, se aprecia tres gráficas del simulador en tres tiempos diferentes; después de 1 hora la infiltración de agua en el medio poroso es leve, después de 5 horas aumenta, mientras que al cabo de 10 horas los niveles de infiltración son muy considerables.
Fig. 3:
Corte transversal de la superficie: proceso de infiltración en el tiempo
3. Conclusiones y recomendaciones
Es posible demostrar la existencia y unicidad de soluciones de la ecuación de Richards tanto para el problema tipo elíptico coma para el problema tipo parabólico cuando las funciones involucradas cumplen con determinadas características.
Para resolver numéricamente la ecuación de Richards bidimensional usando el método de elementos finitos, se debe usar primeramente el método de Crak Nicolson en la variable temporal para lo cual se realiza una discretización en el tiempo. Luego se linealiza la ecuación resultante usando el método Newton. Finalmente se debe discretizar la variable espacial y de esta manera se consigue el esquema numérico que posteriormente se transforma en un programa computacional.
Para realizar pruebas con el simulador se usó un problema cuya solución se conocía de antemano, aquí se pudo determinar que dejando fijo el tiempo y variando el número de elementos de la malla, el error en norma disminuye considerablemente conforme aumenta el número de triángulos de la discretización espacial. Pero, así mismo el tiempo empleado por el computador aumenta considerablemente al aumentar el número de elementos.
Como era de esperarse la discretización de la variable tiempo también influye en la calidad de la solución obtenida, es decir si el paso temporal es demasiado fino, el número de operaciones elementales crece notablemente, para lo cual se necesitaría un computador con mayor capacidad de procesamiento, caso contrario el tiempo de obtención de resultados es demasiado grande.
La infiltración en un medio poroso de cualquier líquido contaminante depende básicamente de dos variables, el tipo de suelo y el tiempo de exposición de éste al contaminante.
A partir de los resultados del simulador, se recomienda realizar pruebas con datos experimentales en diferentes tipos de suelo y con diferentes tipos de contaminantes con el objetivo de contrastar resultados y de ser necesario realizar una calibración del simulador.
En base a los resultados que se obtiene en 
y se recomienda formular un modelo tridemensional y realizar el estudio de la existencia y unicidad de la solución de la ecuación de Richards.
