Download 1. INTRODUCCIÓN El estudio del flujo en rocas fracturadas se

1. INTRODUCCIÓN
El estudio del flujo en rocas fracturadas se constituyó, a partir de la mitad del siglo XX,
como un tópico fundamental para la investigación, debido a que grandes obras civiles se
construyen sobre rocas o constituyen formaciones geológicas almacenadoras de crudo o
de agua, o pueden ser consideradas aptas para la disposición de desechos radioactivos u
otros desechos. Los requerimientos anteriores, más la complejidad del medio, han
conducido a los países desarrollados a promover la investigación en esta área y a realizar
importantes inversiones para la adquisición del conocimiento y el manejo más adecuado
del fenómeno.
El flujo, a través de los espacios vacíos de una roca, ocurre fundamentalmente en un
medio heterogéneo. El estudio del flujo en rocas se ha enfocado hacia la búsqueda de
una mejor representación morfológica del espacio vacío, hacia la estimación de las
propiedades macroscópicas de transporte y hacia la aplicación de ecuaciones de flujo a
dicho medio, para finalmente predecir comportamientos y evaluar problemas en
ingeniería. Obviamente que esto implica la búsqueda de idealización del medio y la
suposición de existencia de una escala característica. Es posible hacer una aproximación
de un medio continuo y definir sus propiedades macroscópicas promedias en estudios
que consideren escalas de longitud mayores a la escala característica. El concepto de la
propiedad macroscópica o parámetro de flujo ha sido utilizado en diferentes teorías de
flujo, ya sea en forma determinística o en forma estocástica, considerando la variación
espacial de la propiedad.
El grado de heterogeneidad y la física del fenómeno cambian con el tipo de roca y con el
tipo de espacios vacíos dentro de ella. En este trabajo se estudia el flujo en rocas
consolidadas que poseen discontinuidades o planos de falla interconectados entre sí, a
través de los cuales fluye y se filtra un fluido. A este medio lo denominaremos medio
fracturado.
El carácter heterogéneo en el medio fracturado es más fuerte que en el medio poroso, lo
que da lugar a que las escalas representativas, si ellas existen, sean mayores que las del
medio poroso, para que la aproximación del continuo sea válida. La existencia de
escalas representativas en el medio fracturado es una aproximación que ha sido
extensamente debatida en la literatura.
Durante las dos últimas décadas, se han
presentado evidencias desde diversas pruebas de laboratorio y de campo, que permiten
vislumbrar la existencia de rasgos heterogéneos que evolucionan con la escala de
observación, lo cual da lugar a problemas en la observación de sus rasgos más
importantes.
Esto conduce a que algunos parámetros geométricos clásicos utilizados
para la caracterización del medio fracturado como son la densidad de fracturas, abertura,
tamaño y orientación, entre otros, pierdan importancia en escalas mayores donde ocurre
el fenómeno de flujo a través de los sistemas de fracturas.
En este trabajo, se utilizan las propiedades especiales que las discontinuidades adquieren
en su génesis para hacer una aproximación de este medio como producto de una
transición de fase o fenómeno crítico.
Esta aproximación permite expresar las
propiedades de flujo del medio fracturado con leyes de escalamiento que se expresan en
función de exponentes críticos, de acuerdo al valor de una escala característica. La
universalidad de estos exponentes constituye una a especial que permite hacer una
interpretación más general, análoga a otros fenómenos críticos.
Los yacimientos o acuíferos, formaciones almacenadoras de fluidos, han sido estudiados
y evaluados a partir de pruebas hidráulicas de diversos tipos y alcances. Éstas pruebas
constituyen un mecanismo confiable y ampliamente utilizado en la estimación de los
parámetros hidráulicos; sin embargo, a pesar de que ellas permiten adquirir información
sobre sus características físicas e hidráulicas en escalas de campo, las pruebas se
encuentran limitadas por los altos costos que implican, tanto su implementación como la
realización de los pozos de observación. Esto obliga a limitar, de algún modo, la
duración y alcance de dichas pruebas.
2
La heterogeneidad es una característica de los materiales de la naturaleza y depende de
la escala y del lente con que se observe.
Existen dos formas de entender la
heterogeneidad, la primera es preguntarse por la presencia de escalas naturales, esta se
estudia con toma de medidas alrededor de un sitio de muestreo y un instrumento
especifico, en este caso es valido hablar de un efecto de escala. La segunda es una
heterogeneidad que no se puede observar a simple vista y debe medirse con instrumentos
que permitan variar el lente de observación bajo una escala especifica. En la practica
para observar este tipo de heterogeneidad se tiene limitación del equipo. Con el fin de
tener un mecanismo que permita observar con diferente escala de observación y
diferente resolución o lente, se diseñó un experimento numérico que permite modelar
pruebas hidráulicas en medios fracturados. El medio fracturado, idealizado como un
sistema de transición de percolación, permite ejecutar diversos tipos de pruebas
numéricas y observar comportamientos transitorios en tantos puntos de observación
como se desee.
Las pruebas numéricas realizadas permiten adquirir suficiente
información sobre un medio con diversos grados de heterogeneidad. Se obtienen, así,
las leyes de escalamiento de los parámetros de flujo, difusividad y conductividad
hidráulica de acuíferos, las cuales reflejan dos formas de observación e interpretación
del medio de acuerdo con su escala característica. Este tema se presenta en el capítulo 4.
En aplicaciones prácticas, como son los bombeos o explotación de pozos en acuíferos y
yacimientos, se presenta un efecto hasta ahora no considerado importante, el efecto no
elástico de las rocas. Este efecto da lugar a que la solución de la ecuación de flujo no
sea la solución clásica de fuente puntual instantánea, solución de Theis o solución Ei,
sino otro tipo de soluciones. Este tema será tratado en la parte final del capítulo 3 con el
interés de presentar los rasgos físicos y matemáticos importantes para la solución de
este problema.
Se utilizan los conceptos de análisis dimensional y grupos de
renormalización para asumir una nueva suposición de autosemejanza, siguiendo las
ideas de Barenblat (1996), la cual llevará a la consideración de una nueva escala
relevante en la solución, la escala del radio del pozo, y su exponente anómalo. Este es
un caso de flujo no lineal en rocas porosas.
3
2. REVISIÓN DE LITERATURA
La importancia que reviste el flujo en rocas fracturadas obliga a realizar una extensa
revisión bibliográfica, que incluye la revisión de trabajos y resultados en las diferentes
escalas espaciales y temporales en las que se presenta el fenómeno. A pesar del esfuerzo
por plasmar en este trabajo el producto de cincuenta años de investigación de la
comunidad internacional, resulta imposible hacerlo en forma más completa y detallada,
sin embargo se espera en lo posible atender con esta revisión, las principales inquietudes
de los lectores, las cuales pueden complementarse en cada una de las referencias dadas.
2.1. EL FLUJO EN FRACTURAS SIMPLES
Inicialmente la investigación del flujo en este medio se enfocó sobre el plano de la
fractura, considerando que el flujo es análogo al que ocurre a través de láminas
paralelas, sobre las cuales se aplican las ecuaciones de Navier-Stokes. Esta ecuación da
un gran peso a la variable geométrica, abertura de la fractura (el flujo es proporciona al
cubo de la abertura), lo cual dio origen a la expresión "ley cúbica", utilizada
extensamente en la literatura. A partir de pruebas de laboratorio en fracturas artificiales
de diferentes materiales en algunos casos y en rocas en otros casos, Lomize (1951),
Romm (1966), Louis (1967), Sharp (1970) investigan el papel que juega la rugosidad de
las paredes de la fractura en la expresión de flujo anterior, expresando el factor de
fricción en forma análoga a como se expresa en el flujo en conductos.
En este caso la escala representativa es la abertura de la fractura y es considerada en los
trabajos anteriores que su dimensión es del orden de micras o algunos milímetros. La
aplicación de los anteriores factores de fricción, considerando flujo laminar o turbulento,
permite hallar el campo de velocidades en la sección de la fractura y permite estimar la
4
permeabilidad de macizos fracturados que presentan una familia de fracturas
horizontales o cuasihorizontales, agregando el efecto de cada una de ellas.
La
estimación se realiza a partir de pruebas de inyección, tipo Lugeon, para las cuales
Heitfeld (1965), Maini (1971), Rissler (1977) y Ewert (1977), entre otros, proponen
métodos de estimación de permeabilidad y evaluaciones de las pruebas considerando
pérdidas de energía en los accesorios y entrada de las fracturas. Para mayor ilustración
se sugiere ver Carrillo (1997) quien realiza una revisión de este tipo de pruebas.
La
pruebas Lugeon se han realizado frecuentemente para estudiar los requerimientos de
impermeabilización de las rocas sobre las cuales se cimentan las presas.
Posterior a estos estudios surgió la necesidad de estudiar la viabilidad de utilizar las
rocas como depósitos de desechos radioactivos, lo cual llevó al incremento de la
investigación en esta área a partir de la década de los años 80. Witherspoon y Wang,
(1980), entre otros autores, proponen la evaluación de la ley cúbica para fracturas
sometidas a esfuerzos, y concluyen que existen áreas de contacto entre las paredes de la
fractura (abertura cero), y que la ley cúbica es válida para cierto rango de abertura y
esfuerzos. Raven and Gale (1985) y Pyrak-Nolte (1987) miden la relación entre los
esfuerzos normales sobre las fracturas y el cerramiento de ellas, y definen un porcentaje
de áreas de contacto. Ellos encuentran que bajo un cierto rango de esfuerzos y aberturas
muy pequeñas, el flujo disminuye más rápidamente que el cubo de la abertura, además
que cuando se excede dicho rango el flujo tiende a un valor que es independiente de los
esfuerzos.
El desarrollo de nuevas tecnologías para la medición de aberturas de fracturas en campo
(pruebas de trazadores, pruebas hidráulicas y medición directa de flujo)
permitió
reportar valores mucho más pequeños que los reportados por Louis (1967) y Romm
(1966). En granitos fracturados en Stripa Mine, Suecia y Fanay-Augere, Francia, sitios
donde se localizaron diaclasas con grupos de rutas de flujo, se estimaron los siguientes
valores: Neretnieks (1982), reportó aberturas de diaclasas del orden de 2 - 10 micras y
Raven (1988) aberturas del orden de 100-200 micras. Neretnieks (1982), reportó áreas
transversales netas de flujo (canales de flujo) del orden de 52 cm2 y Tsang (1987) del
5
orden de 3 cm2. Los caudales de filtración medidos en esos sitios fueron del orden de 115 cm3/min.
Existe suficiente evidencia para concluir que sobre la superficie de la fractura existe una
variación importante de la resistencia al flujo, producto de la variabilidad de la abertura
y de las áreas de contacto, lo cual genera tortuosidad de las líneas de flujo sobre el plano
de la fractura. Tsang (1987) modeló este fenómeno teniendo en cuenta una función de
distribución de aberturas y encontró que el flujo del modelo difiere del flujo calculado
con la ecuación cúbica; Tsang empieza a hablar de canales de flujo en el plano de la
fractura. Desde una técnica de inyección de una aleación especial dentro de fracturas,
Pyrak-Nolte (1987), obtiene moldes de los espacios vacíos de una fractura en laboratorio
y muestra que ellos se encuentran en forma irregular y forman canales tortuosos que
constituyen rutas o canales de flujo. A partir de estudios de laboratorio en núcleos de
rocas con fracturas (Gentier, 1989) se muestra una variación log normal de la abertura
sobre el plano de la fractura.
El modelamiento numérico de este fenómeno lleva a
aplicar diferentes distribuciones que representan variabilidad de los planos de la fractura
ya sean como superficies fractales (Brown, 1987) o superficies con variación espacial
representada con geoestadística (Tsang, 1998).
2.2. SISTEMAS DE FRACTURAS
La necesidad de estudiar el flujo en escalas mayores llevó a considerar la existencia de
familias de fracturas generadas por campos regionales de esfuerzos, que pueden
presentarse con diversos direccionamientos. En este caso se dio importancia a variables
geométricas como distancia entre fracturas, persistencia de la fractura y por supuesto la
orientación de dichas familias. Estas variables dan lugar a la interpretación que se usa
en geología con el diagrama de Schmidt (proyección estereográfica) para caracterizar los
grupos de fracturas y su densidad espacial.
Los primeros modelos geométricos que representan la morfología de los planos de
fracturas, son modelos bastante idealizados. El primero es un modelo ortogonal, que
6
considera los sistemas de fracturas como grupos mutuamente ortogonales, los
planteamientos de Snow (1965) se basan en este modelo.
El modelo de Beacher (1978) considera que las fracturas tienen formas circulares o
elípticas y están representadas por un grupo de discos de radio que puede ser constante o
tener alguna variación representada con una distribución de probabilidad, la localización
del centro de cada disco puede tener un patrón regular o distribuirse como un proceso
Poisson. Este modelo fue ampliamente utilizado, entre otros autores, por Long (1982), y
Cacas, Ledoux y Marsily (1990). Otros modelos menos utilizados son los de Veneziano
(1978) y Dershowitz (1984), los cuales se basan en la similaridad entre la geometría de
las fracturas observadas en campo y la geometría de los planos y líneas que siguen un
proceso Poisson. Estos modelos son utilizados para plantear los modelos de flujo que
hacen una consideración discreta de las fracturas.
Desde el punto de vista del flujo, un macizo fracturado puede ser considerado
homogéneo y anisotrópico cuando se puede construir un tensor con sus valores
principales de conductividad hidráulica (Bear, 1972). Esto se puede lograr de dos
maneras, considerando el medio fracturado equivalente a un medio poroso y aplicando
los métodos convencionales para obtener sus parámetros a partir de pruebas hidráulicas
(Hsieh, 1985) o considerando un medio con fracturas continuas (infinitas) perfectamente
formadas y utilizando metodologías para el estudio de macizos discretamente
fracturados planteando el tensor de conductividad propuesto por Snow (1969).
El
primer enfoque ha sido ampliamente utilizado, pero requiere la suposición de que el
medio pueda ser representado con un volumen elemental de referencia.
Las diaclasas de longitud infinita, fueron tratadas por Snow (1969), quien desarrolló una
expresión matemática para obtener el tensor de conductividad hidráulica de una fractura
simple con orientación arbitraria y abertura constante, bajo un sistema de coordenadas
dado.
El tensor para una red de fracturas, se forma agregando las respectivas
componentes de los tensores de cada fractura individual, o sea que la conductividad
7
hidráulica equivalente puede encontrarse simplemente sumando la contribución de cada
una de las fracturas discretas.
Las diaclasas finitas, son tratadas por Long (1982, 1981), quien genera un modelo
numérico de arreglos de fracturas a partir del modelo geométrico de Beacher. Sobre
estos modelos se imponen ciertas condiciones de frontera para generar un gradiente
lineal sobre una región en la cual se generan sistemas de diaclasas y resuelve la
ecuación de Navier-Stokes para cada uno de los elementos, obteniendo tensores de
permeabilidad en régimen permanente.
El tratamiento del medio fracturado mediante el modelo de un medio continuo
equivalente ha sido ampliamente utilizado. Sin embargo, la discusión sobre la existencia
de un volumen elemental de referencia que permita dar validez al anterior enfoque ha
generado controversia; vale la pena estudiar los planteamientos de Neuman (1987) y
Marsily (1990).
Las observaciones de campo y de laboratorio demostraron que los sistemas de fracturas
en rocas están lejos de comportarse como un medio equivalente al medio poroso
homogéneo o como un medio con fracturas continuas perfectamente formadas y
conectadas entre sí.
Las evidencias son claras en demostrar que las fracturas no
conforman planos de flujo, Pyrak-Nolte (1987) y que ellas no se conectan
completamente entre sí, Hsieh (1998). El flujo ocurre a través de rutas preferenciales,
ya sea a través de canales tortuosos que se interceptan con otros canales, formando
grupos de canales de flujo (clusters) o formando especies de canales que se interceptan
en menor grado, en cuyo caso podría considerarse que el flujo ocurre esencialmente a
través de pocos canales (Tsang y Neretnieks, 1998).
Por esta razón, se piensa que el flujo en una red de fracturas puede aproximarse mejor a
una red de líneas idealizadas como elementos conductores más o menos conectados en
dos o tres dimensiones (Billaux, 1989; Cacas y Marsily, 1990; Odling, 1991). Los
modelos de redes basados en las ideas anteriores son comunes. Long y Whiterspoon
8
(1987) desarrollan una relación entre la permeabilidad de un modelo de fracturas como
conductores de longitud y densidad dada por una distribución Poisson.
La difusividad y la conductividad hidráulica, principales parámetros que permiten
desarrollar los modelos de flujo e interpretar las pruebas hidráulicas, han sido estudiadas
a escalas de laboratorio, a escalas de ensayos puntuales y a escalas de campo,
presentando alta variabilidad. Gelhar (1986), Dagan (1986) y Neuman (1998) modelan
estos parámetros como variables estocásticas que poseen variación espacial
representada mediante variogramas y geoestadística. Existen diversos modelos de este
tipo aplicados a medios fracturados heterogéneos como los desarrollados por Folli
(1990) y Niemi (1994).
2.3 CONECTIVIDAD DE FRACTURAS
La alta heterogeneidad que presenta el medio fracturado y la falta de correlación entre
los parámetros geométricos tradicionalmente medidos (abertura, distanciamiento,
longitud, densidad de fracturas) y los parámetros de flujo, sugieren que se está tratando
un medio básicamente desordenado, el cual, difícilmente puede ser representado en
forma determinística por parámetros geométricos observados a simple vista. Se hizo por
lo tanto indispensable acudir a otro parámetro que permitiera obtener información
adicional sobre los sistemas de fracturas como es su grado de conexión o conectividad.
Este ha sido un término poco aceptado por las escuelas tradicionales de investigación,
pero imposible de ignorar dadas las evidencias encontradas.
Billaux (1989) habla de la conectividad de las redes de fracturas, y propone algunas
metodologías para tratar de evaluar este efecto. Utilizando modelos de sistemas de
fracturas y generando un ensayo numérico de inyección, Billaux, 1990, determina el
número de intersecciones alcanzadas, esto da lugar a una curva llamada de conectividad
y define la conectividad como el numero de rutas a través de las cuales el agua puede
fluir cuando es inyectada en un punto del sistema.
9
Para Billaux existen tres casos representativos de conectividad:
a. El sistema es pobremente conectado:
El número de fracturas que pueden ser
alcanzado desde un punto de inyección es limitado, y la red se encuentra por debajo
de la densidad crítica.
b. La red está muy bien conectada: El número de fracturas que puede ser alcanzado es
grande, y la inyección llega hasta las fronteras de la red. Se considera que la red
está por encima de la densidad critica.
c. Entre estos dos casos extremos, existe el caso en el cual se alcanza cierto número de
fracturas y la respuesta de la inyección no llega hasta todas las fronteras.
Se
considera que la red está cercana a la densidad crítica. En este caso densidad es el
número de fracturas por unidad de área y densidad crítica es el número mínimo
necesario para que ocurra el flujo..
La complejidad del flujo en rocas fracturadas se debe a la geometría compleja que
presentan los sistemas, ellas no tienen longitud infinita sino que pueden presentarse en
muy diversas escalas, por lo tanto la permeabilidad de la roca depende del grado de
interconexión del sistema de fracturas.
La permeabilidad debe tener una dependencia
de la escala, la cual ha sido observada desde escalas de laboratorio hasta la escalas de
campo (Brace, 1980, 1984; Clauser, 1992). Otras aproximaciones de la conectividad se
estudian con la aplicación de la teoría de percolación a los sistemas de fracturas.
2.4 TEORÍA DE PERCOLACIÓN EN FRACTURAS
El estudio de los sistemas de fracturas con la teoría de percolación se inicio con los
trabajos de Robinson (1983, 1984) y Englman (1983), quienes modelaron las fracturas
como segmentos de líneas de posición aleatoria y distribución estadística de la longitud
y la orientación. Ellos relacionan la densidad de líneas en una región con la densidad de
sitios ocupados en una malla de percolación.
10
Berkowitz y Balberg (1993) presentan una revisión de la literatura de percolación
aplicada al flujo en
medio poroso y fracturado. Stauffer (1992) presenta los
fundamentos básicos de la teoría de percolación. Sahimi (1994 y 1995) presenta una
revisión extensiva de los conceptos de la teoría de percolación y la aplicación a medios
porosos y fracturados, incluyendo también la aplicación de la teoría del medio
equivalente desarrollada por Kirkpatrick (1983) la cual se origina igualmente en las
ideas de percolación, pero considera sistemas alejados del umbral de percolación.
Robinson (1984) desarrolla simulaciones numéricas para encontrar algunos exponentes
críticos y la densidad crítica de las fracturas; en el límite de sistemas grandes la densidad
crítica es aproximadamente 5.6 en sistemas de dos dimensiones, valor confirmado
posteriormente por Berkowitz (1995) y Bour and Davy (1997). Hestir and Long (1990)
utiliza la teoría de percolación a partir de un modelo de fracturas mapeado sobre mallas
de percolación para obtener una expresión que relacione la permeabilidad con la
conectividad del sistema la cual se expresa en función de la probabilidad de ocupación.
Los resultados obtenidos se comparan bastante bien con los resultados obtenidos con el
modelo de Robinson. A partir de estos modelos halla el exponente de escalamiento de la
permeabilidad desde simulaciones de percolación.
Los modelos anteriores utilizan
elementos o fracturas de longitud constante. Bour and Davy (1997) utiliza la teoría de
percolación en sistemas de fracturas de longitudes variables que siguen una ley de
potencia, ellos plantean las propiedades de escalamiento de la red y definen longitudes
criticas.
Charlaix (1987) trabaja exponentes críticos en la escala de la fractura y
comportamientos de permeabilidad en escalas de la longitud de correlación
considerando una distribución espacial de resistencia hidráulica y de resistencia crítica, y
expresa la longitud de correlación en función de esta distribución.
La medida de la conectividad, en este contexto, ha sido propuesta como el número
promedio de intersecciones por fracturas, en un proceso aleatorio de generación de
líneas.
Si se escoge una determinada línea se cuenta el número de líneas que la
interceptan, y se define la conectividad
como el valor promedio de líneas que
interceptan dicha línea seleccionada.
11
2.5 PRUEBAS HIDRÁULICAS
La distribución de la heterogeneidad en el espacio es reflejada por pruebas hidráulicas de
campo. Hsieh (1998), Neretnieks (1987), entre otros, concluyeron que no todas las
fracturas conducen fluido, y que ellas se agrupan en especies de clusters conductores
que llevan la mayor cantidad de flujo, o se agrupan en clusters aislados, que no
conducen. Estas pruebas son realizadas en ensayos que incluyen varios pozos de
observación además del pozo de bombeo o inyección, en los cuales se aíslan las zonas
de fracturas conductoras con empaques dentro del pozo. Vale la pena comentar que la
tecnología para los ensayos hidráulicos en campo ha avanzado considerablemente y se
ha aplicado a rocas fracturadas en varias regiones de Estados Unidos y Europa.
Para simular el fenómeno de flujo que ocurre a través de rutas preferenciales formando
caminos intrincados de flujo, se han realizado simulaciones numéricas. Polek (1990) y
Acuña y Yortsos (1995), representaron este medio con sistemas fractales conocidos a
los cuales les aplicaron pruebas de flujo y pudieron observar comportamientos y
respuestas de esas simulaciones.
Barenblat, Zheltov and Kochina (1960) introdujeron el concepto de doble porosidad para
estudiar el flujo en rocas fracturadas en el caso en que la porosidad de la roca no fuera
despreciable.
Rofail (1967) y Barenblat (1990) presentan una metodología para
interpretar pruebas de pozos basados en el concepto anterior de doble porosidad.
Warren and Root (1963) presentaron otros parámetros para caracterizar las formaciones
anteriores. Todos estos autores asumieron flujo cuasi permanente en la matriz de la
roca. Gringarten (1982) introdujo nuevos grupos adimensionales y presentó curvas tipo
para analizar datos de pruebas de pozos. Una revisión completa de análisis de pruebas
de pozos en el área de petróleos se puede encontrar en Earlougher (1977) y Streltsova
(1988). Karasaki (1987) aplica las ecuaciones de flujo radial en acuíferos fracturados
considerando dos dominios diferentes, el primero es un dominio afectado por la
influencia inmediata del pozo donde las suposiciones tradicionales del flujo radial hacia
12
el pozo no son válidas; en el segundo dominio, alejado del pozo, puede asumirse un
medio homogéneo sobre el cual se aplican las ecuaciones tradicionales. Karasaki, Long
y Whiterspoon (1988) realizan una revisión de los modelos analíticos para interpretar las
pruebas Slugs, las cuales se basan en la medida de recuperación de niveles en el pozo.
Barker (1988) considera la dimensión no entera en el dominio de flujo desde un pozo y
la involucra en las ecuaciones de flujo resolviendo con transformadas de Laplace. Polek
(1990) hace la comparación entre la dimensión no entera del flujo y la dimensión no
entera del medio, sugiriendo que la primera es siempre menor a la segunda.
La necesidad de incorporar escalas adicionales, que no están involucradas en las
ecuaciones de flujo sino quizás en la condición inicial, ha sido propuesta por Barenblat
(1987) en el área del flujo de fluidos y por Goldenfeld (1999) en el área de las
transiciones de fase, a partir del planteamiento de grupos adimensionales que siguen una
ley especial de autosemejanza y para su solución se debe acudir a consideraciones
adicionales al análisis dimensional.
Estas teorías consideran las propiedades de
renormalización y homogeneidad dimensional de las ecuaciones de flujo.
Como
ejemplo, se tiene el problema de difusión de una masa de fluido dentro de un medio
poroso considerando retención residual, en cuyo caso, la solución involucra el radio de
la recarga como escala adicional (Barenblat, 1987) generando un componente no lineal
que da lugar a un problema de coeficiente discontinuo y solución de un problema de
autovalor (Gómez y Mesa, 1997).
La incorporación de escalas adicionales se debe a la
no linealidad del flujo en
fenómenos transitorios, lo cual da lugar al planteamiento y solución de problemas con
coeficientes discontinuos (Barenblat, 1987). La no linealidad se puede presentar en el
flujo hacia el pozo debido a la presencia de turbulencia en el dominio cercano al mismo
acuífero (Pérez, 1995).
13
3. FUNDAMENTOS BÁSICOS DEL FLUJO EN ROCAS
Los tipos de rocas existentes en la naturaleza se describen de acuerdo a su origen
geológico como rocas ígneas, sedimentarias y metamórficas. Este trabajo se enfoca
principalmente sobre rocas ígneas y metamórficas que han sido sometidas a campos de
esfuerzos, generando en ellas discontinuidades a las que hemos dado el nombre de
fracturas. Se pueden resaltar los granitos, neises, areniscas, basaltos, esquistos y otras,
las cuales conforman el llamado medio o rocas fracturadas.
El flujo en rocas fracturadas constituye un componente importante en el estudio de las
aguas subterráneas pues pueden llegar a formar acuíferos de alta capacidad y convertir al
flujo en un recurso hídrico relevante. En otros casos, el flujo se presenta en cantidades
limitadas, observándose solamente algunos hilos de agua en túneles y afloramientos de
roca, que no tienen ningún atractivo como recurso hídrico, pero que ocasionan serios
problemas de estabilidad de grandes masas rocosas e igualmente, pueden generar
filtraciones indeseables a través de presas, túneles y depósitos de desechos radioactivos.
Cuando el fluido es petróleo, se habla de yacimientos y se está interesado en su
explotación, la cual debe ser optimizada, dado el valor que tiene este recurso.
El estudio del flujo en medio fracturado, rocas consolidadas, presenta mayor
complejidad que el estudio del flujo en medio poroso pues sus parámetros hidráulicos
están sometidos a mayor incertidumbre, debido a la fuerte heterogeneidad que
caracteriza este medio.
Por lo tanto, se exige una mayor profundidad en el
reconocimiento de los rasgos importantes del fenómeno y en la identificación de las
escalas representativas, ya sea que se estudie el flujo a través de una fractura simple o a
través de un sistema de fracturas.
14
Para describir los fenómenos de transporte, como son
la conducción de calor, la
conducción eléctrica o la transferencia de masa, se utilizan las leyes de Fourier, Ohm y
Fick, respectivamente, leyes de flujo que tienen básicamente la misma forma (los flujos
son proporcionales a los gradientes o potenciales). Éstas son válidas en el régimen
considerado cerca al equilibrio, es decir, donde los flujos son funciones lineales de las
fuerzas. Estas leyes dan origen a la formulación matemática de los problemas y al
planteamiento de ecuaciones de flujo cuyas soluciones pueden ser obtenidas en forma
analítica si se conocen las condiciones iniciales y de frontera. Dichas soluciones se
pueden obtener desde la aplicación del análisis dimensional, cuyo concepto fundamental
es la autosemejanza de los fenómenos físicos.
3.1 CONCEPTOS FUNDAMENTALES
En el medio fracturado, igual que en el medio poroso, el estudio del flujo requiere la
determinación de la conductividad hidráulica o coeficiente de permeabilidad K, el cual
se define como la capacidad de un medio para permitir el paso del agua.
Este
coeficiente refleja ambas propiedades, las del medio y las del fluido y puede ser
expresado como
K =
kρ g
μ
(3.1)
Donde k es la permeabilidad del medio, llamada permeabilidad intrínseca, g es la
aceleración de la gravedad, ρ es la densidad del fluido y μ es la viscosidad dinámica del
fluido.
El flujo en el subsuelo, en la zona saturada, ha sido tradicionalmente estudiado
utilizando la Ley de Darcy (ley de flujo), la cual originalmente se postuló para flujo en
una dirección y cuyo rango de validez se encuentra en el régimen de flujo laminar,
V = KJ,
15
(3.2)
donde V es la velocidad de descarga, K es la conductividad hidráulica y J es el gradiente
hidráulico.
Para describir el flujo en un dominio de medio poroso, se ha utilizado el concepto de
Volumen Elemental de Referencia, VER (volumen para el cual el promedio de la
conductividad hidráulica se estabiliza) el cual se observa en forma esquemática en la
Figura 3.1. De existir, a este volumen de roca se le puede asignar un valor promedio de
conductividad hidráulica (Marsily,1986).
k
Tamaño del Ver (Vol., L)
L
Figura 3.1. Esquema explicativo del VER. Bear (1972).
El concepto del VER permite reemplazar el medio real, por un modelo teórico de
MEDIO CONTINUO sobre el cual se pueden asignar variables promedias, como la
porosidad, conductividad hidráulica, etc., y permite describir el flujo dentro de un
dominio por medio de ecuaciones diferenciales.
En medios que presentan alta
heterogeneidad como es el medio fracturado, el VER podría no existir o ser mucho
mayor que la escala de observación del problema. Se han planteado discusiones sobre
su existencia en estos medios (Neuman, 1987).
La ecuación diferencial en coordenadas cartesianas que representa el flujo estacionario
en un medio continuo es
∂
[ K xx J x ]+ ∂ [ K yy J y ]+ ∂ [ K zz J z ]= 0
∂x
∂y
∂z
16
(3.3)
Jx, Jy y Jz: Gradiente hidráulico en las tres direcciones.
Un medio homogéneo y anisotrópico se describe con el tensor de conductividad
hidráulica. La medida direccional de la conductividad hidráulica en un dominio permite
encontrar este tensor y el dibujo de las respectivas conductividades direccionales Kxx,
Kyy, Kzz, en coordenadas polares conformarán un elipsoide si el medio es homogéneo y
anisotrópico (Bear, 1972).
El estudio del flujo en el plano de una fractura se basa en las ecuaciones de Navier
Stokes, que representan el comportamiento del flujo entre dos láminas paralelas. En esta
analogía el flujo es considerado laminar con una distribución de velocidades parabólica
en la sección del flujo. En la figura 3.2, se observa esta aproximación.
Figura 3.2. Esquema de flujo en planos de fracturas continuas.
Aplicadas las condiciones de borde, al plano de la fractura, la velocidad y el caudal se
expresan como
V x=
g b 2 ∂p
,
12 ν ∂ x
Q =
gHb
12 ν
3
∂p
,
∂x
(3.4)
Donde b es la abertura de la fractura; H es el ancho de la fractura: g es la gravedad, ν es
la viscosidad cinemática y
∂p
es el gradiente de presión en la dirección del flujo. La
∂x
17
expresión 3.4 es una relación cúbica entre la descarga y la abertura de la fractura que ha
sido ampliamente utilizada en el estudio del flujo en rocas fracturadas (Long, 1982;
Marsily, 1986; Wittke, 1990).
3.2 PLANTEAMIENTO DEL PROBLEMA DEL FLUJO EN MEDIOS
FRACTURADOS
El medio fracturado tiene un carácter heterogéneo mucho más marcado que el medio
poroso, por lo tanto, el hecho de promediar el valor de la conductividad hidráulica sobre
un dominio o el considerar un área transversal homogénea de flujo, nos aleja, en este
caso, mucho más de los comportamientos reales de flujo.
A continuación se expresan algunos de los rasgos observados en escala de la fractura y
en escala de sistemas de fracturas desde observaciones de campo, trazadores y pruebas
hidráulicas. Sin embargo, en cada escala de observación se pueden encontrar nuevos
rasgos. Esto sugiere que la observación de rasgos importantes en el fenómeno de flujo
en rocas fracturadas varía de acuerdo con la escala que se utilice para hacer dicha
observación.
3.2.1 Condiciones de flujo en la escala de la fractura
La fractura está formada por dos superficies onduladas y rugosas que están en contacto
en algunos puntos (abertura cero) y separadas en otros un valor del orden de unas pocas
micras hasta el orden de unos pocos milímetros. Las características físicas observadas
en pequeñas escalas como son la rugosidad de sus paredes, las áreas de contacto, el lleno
en la fractura y otras, generan heterogeneidad en la escala de la fractura que afectan las
condiciones del flujo.
Las áreas de flujo a través de la fractura (secciones abiertas) son muy variables y son
consideradas como canales potenciales de flujo.
18
Figura 3.3. Esquema de canales de flujo en el plano de la fractura
El flujo en la fractura ocurre si la sección abierta se encuentra conectada con otras
secciones abiertas de la misma fractura, que le permita formar un canal de flujo. Bajo
un gradiente hidráulico impuesto, el agua buscará el camino más fácil y se crearán rutas
tortuosas, evitando las secciones de mayor resistencia. Por lo tanto el flujo en un plano
de fractura ocurrirá a lo largo de canales que pueden continuar aislados o pueden
cruzarse con otros como se observa en la figura 3.3.
Las anteriores características han sido observadas a escala de las fracturas y existen
evidencias físicas de estos fenómenos basadas en pruebas de trazadores, en pruebas
hidráulicas y en mediciones de caudales de filtración realizadas en granitos fracturados
en Stripa Mine, Suecia y Fanay-Augere, Francia (Tsang, Neretnieks, 1998).
Una
revisión de este tema se puede encontrar en Gómez (1994).
3.2.2
Condiciones de flujo en la escala de sistemas de fracturas
Las fracturas se presentan en familias que pueden interceptarse y formar especies de
bloques, los cuales debido a la presencia de puentes de roca y áreas de contacto, están
lejos de conformar bloques regulares (los cuales supuestamente constituirían sistemas
con perfecta conexión). Cuando dos planos de fracturas se interceptan pueden presentar
rasgos especiales en la escala de observación. Bajo cierta escala se podrían observar los
siguientes casos mostrados en forma esquemática en la figura 3.4.
19
-
Un canal de flujo en una fractura puede estar conectado a una parte cerrada o
parcialmente cerrada en otra fractura. En este caso la intersección impide o controla
el flujo desde un plano al otro. No existe conexión en la intersección, no existe flujo,
o este ocurre en forma controlada. Si esta condición prevalece, existirán muy pocos
canales que conecten las fracturas. Poca Conexión.
-
Un canal en una fractura puede estar conectado con una parte abierta en la otra
fractura y en este caso el canal continuará sobre más de una fractura. Se presenta
mezcla de canales en las intersecciones y buena conexión.
Figura 3.4. Esquema de intersección de planos de fracturas. El flujo ocurre solamente a
través de canales de flujo que se interconectan entre sí.
En el primer caso existirán muy pocos canales que estén conectados de una fractura a
otra y estos pueden tratarse como un grupo de conductos individuales (existirá poca
conexión). En el segundo caso se tendrá una mezcla de canales en las intersecciones
entre fracturas (existirá una conexión fuerte). Los macizos fracturados poseen una
mayor o menor comunicación hidráulica, que depende de las condiciones de los
esfuerzos tectónicos a que fueron sometidos durante el tiempo geológico y dependen
también de las condiciones físicas de las rocas afectadas. El flujo ocurre a través de una
intrincada red de canales, algunos de ellos no conducen a ninguna parte y forman
clusters que no contribuyen al flujo, solamente una parte de la red de canales contribuirá
20
al flujo en un macizo.
Estas características permitirán expresar el flujo en medio
fracturado como un fenómeno crítico o transición de percolación como se presentará en
el numeral 3.8.
3.2.3 Ecuaciones de flujo en rocas
Para la interpretación de pruebas hidráulicas en formaciones geológicas que constituyen
acuíferos o yacimientos (formaciones que almacenan fluidos) se deben estudiar las
ecuaciones que rigen el flujo hacia o desde un pozo. Estas pruebas se realizan con el
objetivo de conocer el comportamiento hidráulico de dichas formaciones y estimar los
parámetros de flujo.
La interpretación de los procesos de inyección o de bombeo en un pozo utilizando la
ecuación de Theis o la ecuación EI (integral exponencial), se realiza bajo la suposición
de flujo radial en medio homogéneo, isotrópico y dominio infinito. Bajo la hipótesis de
homogeneidad, la ecuación que rige este flujo es la
ecuación de difusión
en
coordenadas cilíndricas y su solución está basada en la suposición de que el pozo es una
fuente lineal instantánea de radio cero. La solución se utiliza para la interpretación de la
prueba, sea ésta de bombeo o de inyección, es decir, no se diferencian los procesos de
entrada o de salida de fluido a la formación geológica, por lo tanto, el bombeo o
inyección son procesos tomados como dos procesos idénticos y se representan con la
misma solución.
La formulación matemática del problema se realiza inicialmente en coordenadas
cartesianas y en régimen transitorio, para lo cual se aplica el principio de conservación
de masa sobre un volumen de control que se escribe como
−
⎡∂
⎢⎣ ∂x ( ρ Vx ) +
∂
∂y
( ρ Vy ) +
∂
∂z
⎤
⎥⎦
( ρ Vz ) =
∂ρ
,
(3.5)
∂t
Se aplica el principio de conservación de momentum (ley de flujo o ley de Darcy)
considerando las componentes cartesianas del potencial de velocidad ϕ, para obtener
21
Vx = − k
∂h
∂x
Vy = − k
,
∂h
Vz = − k
,
∂y
∂h
∂z
.
(3.6 )
Cuando el flujo ocurre bajo condiciones isotérmicas, no se requiere el principio de
conservación de energía.
Si se tiene en cuenta la compresibilidad del acuífero βr α, y la compresibilidad del fluido
βf , el término de la parte derecha de la ecuación (3.5) expresa la parte de fluido
almacenado en la formación y toma la forma
∂ρ
∂t
⎛
⎜
⎝
= mρ⎜β
+
f
β r ⎞⎟ ∂p
m
⎟ ∂t
⎠
,
(3.7 )
donde m es la porosidad, ρ es la densidad del fluido y p es la presión. El coeficiente de
almacenamiento S, se puede expresar en función de las propiedades físicas de la roca
como
⎛
⎝
S = γ m B⎜ β
f
+
βr ⎞
m
⎟
⎠
.
(3.8)
De esta manera el coeficiente de almacenamiento representa la variación del peso del
fluido en el tiempo: ∂W/∂t = Ss B, siendo Ss el coeficiente de almacenamiento específico
(unidades de longitud) y B el espesor del acuífero.
Teniendo en cuenta la ley hidrostática ∂p = γ ∂h, y la densidad constante en el caso de
fluido incompresible, las ecuaciones (3.5), (3.6) y (3.8) se combinan para obtener la
ecuación general de flujo,
⎛ ∂Vx
⎝ ∂x
−⎜
+
∂Vy
∂y
+
∂Vz ⎞
∂z
⎟=
⎠
S ∂h
B ∂t
.
Para el flujo Darciano existe un potencial de velocidad
(3.9 )
en función de la cabeza
hidráulica, ϕ = K∇h y por lo tanto la ecuación diferencial para flujo transitorio en un
22
acuífero confinado de espesor uniforme, isotrópico y homogéneo en coordenadas
cartesianas será
2
2
2
∂ h ∂ h ∂ h
S ∂h
+
+
=
.
2
2
2
kB ∂t
∂x
∂y
∂z
(3.10)
En coordenadas cilíndricas la ecuación anterior se expresa
2
∂ h
∂r
2
+
1 ∂h
r ∂r
=
S ∂h
.
kB ∂t
(3.11)
En el caso de acuífero libre usualmente se ignoran los efectos de la elasticidad del
acuífero considerados pequeños y el coeficiente de almacenamiento se toma como
equivalente a la porosidad efectiva (volumen que puede ser extraído por gravedad desde
un volumen unitario de acuífero saturado).
En un acuífero confinado la cantidad total de agua que se puede extraer, es liberada
debido a los efectos elásticos del acuífero y del fluido como se observa en la expresión
(3.8) que representa el coeficiente de almacenamiento.
Theis (1935) propone una solución a la ecuación de flujo en acuíferos haciendo una
analogía con la ecuación de calor en un plano infinito a partir de una fuente puntual
instantánea estudiada por Carslaw and Jaeger (1959). Debido a que esta es una solución
a una ecuación de difusión clásica se hace la suposición de fuente infinitesimal (fuente
concentrada instantánea) y por lo tanto no se involucra la escala del radio del pozo.
La solución tiene la forma
Q
ho − h =
4π T
∞
r2 S
e −u
∫u u ∂ u , u = 4Tt
(3.12)
Esta solución ha sido la base para la hidráulica de pozos y el tratamiento matemático de
acuíferos durante todo el siglo XX. Teóricamente la solución se aplica a formaciones
geológicas homogéneas de extensión infinita, en las cuales el pozo atraviesa todo el
23
espesor del acuífero, la transmisividad es constante en el tiempo y en el espacio, y el
pozo tiene un diámetro infinitesimal.
3.2.3.1 Analogía entre el flujo de calor en materiales y el flujo en acuíferos
Las ecuaciones anteriores se basan en una analogía total con el flujo de calor en cuerpos
sólidos. El principio de conservación de calor y la Ley de Fourier dan origen a la
ecuación diferencial para conducción de calor, que en una dirección es
∂θ
∂t
=κ
2
∂ θ
,
2
∂x
κ =
λ
ρC
,
(3.13)
Donde κ es la difusividad térmica, λ la conductividad térmica, ρ la densidad y C el calor
especifico, considerado constante e idéntico ya sea para calentamiento o enfriamiento de
cuerpos sólidos.
El análogo en acuíferos es κ = T/S la difusividad hidráulica o coeficiente de difusividad,
siendo T la transmisividad del acuífero (análoga a la conductividad térmica, λ) la cual
describe la habilidad del acuífero para transmitir agua como un todo y es igual a la
conductividad hidráulica por el espesor saturado, kB. Por lo tanto el coeficiente de
almacenamiento S resulta ser análogo al calor específico, C.
Cuando existe una variación de la temperatura el flujo de calor es relajado en forma
instantánea, lo cual resulta ser análogo con el caso de acuíferos confinados puesto que el
agua se descarga instantáneamente desde el almacenamiento debido a la caída de
presión. En acuíferos libres el drenaje ocurre en forma lenta y el tratamiento dado por
Theis (1935) no tiene en cuenta ni el tiempo de demora en la descarga desde el
almacenamiento, ni el volumen de agua retenido en los poros.
En ambos tipos de acuíferos los coeficientes utilizados son considerados constantes e
idénticos en los casos de extracción o inyección; igual consideración se hace para la
interpretación de pruebas hidráulicas en yacimientos que almacenan petróleo.
24
En estos yacimientos, la ecuación que rige la distribución de presión p(r, t) en flujo
radial es la siguiente
∂ 2 p 1∂p 1 ∂p
+
=
∂r 2 r ∂r n ∂t
n=
k
,
mc t μ
(3.14)
en donde n es el coeficiente de difusividad, ct es la compresibilidad total (incluye la roca
y el fluido), μ viscosidad dinámica y m porosidad (Streltsova, 1988). Esta ecuación
tiene la misma forma de la ecuación utilizada en acuíferos.
3.2.3.2 Solución analítica a la ecuación de calor
Ecuaciones del tipo 3.11, 3.13 o 3.14, requieren condiciones iniciales y de borde. La
condición inicial consiste en el valor de la función θ(x, t) en el tiempo inicial. Las
condiciones de borde pueden diferir dependiendo de las condiciones térmicas en las
fronteras, ya se trate de un valor dado en los extremos, valores de derivadas de la
función o alguna relación lineal entre las derivadas y la función (Tijonov y Samarski,
1972).
Una vez se tiene definido el problema de valor de frontera se busca la solución
utilizando diversos métodos como son las transformadas de Laplace y el análisis
dimensional (ver literal 3.5). El método de separación de variables busca la solución de
la ecuación homogénea y aparece una serie que converge. Como la Ley de Fourier en la
ecuación de momentum es lineal se puede usar el principio de superposición
resolviendo se llega a la siguiente expresión
θ ( x, t ) =
l
∫ G ( x, ξ , t ) ϕ (ξ ) ∂ ξ
(3.14.1)
o
Siendo G (x, ξ , t) la llamada función de Green o “función fuente puntual instantánea”
la cual representa la distribución de la temperatura en una barra de longitud, L, en el
tiempo t, siempre y cuando la temperatura en el tiempo inicial sea igual a cero y una
cierta cantidad de calor sea liberada instantáneamente en el punto x = ξ.
25
La búsqueda de soluciones de semejanza permite analizar condiciones importantes de
los problemas y obtener la solución a ecuaciones diferenciales parciales. El tema de
autosemejanza y análisis dimensional se desarrolla en el numeral 3.5 y 3.6.
3.3 PLANTEAMIENTO MATEMÁTICO PARA FLUJO EN MEDIO POROSO
FRACTURADO
En esta investigación se plantean las consideraciones de flujo de medios puramente
fracturados y de medios llamados de doble porosidad. Sin embargo, los aportes finales
de este trabajo no consideran el medio de doble porosidad. Cuando se considera que el
fluido se almacena en el espacio de los poros entre los granos (medio poroso) o matriz
del medio, se tiene el flujo en medio poroso.
Pero si además esa matriz posee
discontinuidades que forman una red intrincada de conductores comunicados entre sí
que poseen algún valor de conductividad hidráulica y permiten el flujo a través de ellos,
se habla entonces del flujo en un medio poroso fracturado o flujo en un medio de doble
porosidad.
En el segundo caso las discontinuidades poseen generalmente dos
dimensiones de varios órdenes de magnitud mayor que la tercera. Mientras en el primero
los poros presentan tres dimensiones del mismo orden de magnitud.
Se debe distinguir entonces entre un medio poroso fracturado en el cual los dos medios
poseen porosidad y permeabilidad. Un ejemplo es la arenisca, en cuyo caso el volumen
almacenado en las fisuras puede ser despreciable con respecto al volumen total de la
matriz. Un medio puramente fracturado se presenta cuando la porosidad de la matriz es
nula y el sistema de fracturas conforma bloques impermeables que no intercambian
flujo, por lo tanto, el flujo ocurre solamente a través del sistema de fracturas; un ejemplo
lo constituyen las rocas ígneas fracturadas como el granito. Estas simplificaciones se
pueden hacer sólo bajo ciertas condiciones.
Debido a que en la mayoría de los casos la conductividad hidráulica del sistema de
fracturas es muchas veces mayor que la conductividad de la matriz porosa, se dice que el
fluido es almacenado en los bloques porosos y luego transportado a través de las
26
fracturas.
Este hecho es importante de considerar cuando se estudian procesos
transitorios debido al efecto del intercambio de flujo entre los dos medios. Por otro lado,
resulta fundamental la consideración de las diversas escalas a considerar en esta
aproximación.
Existen tres escalas fundamentales: La escala del dominio o región de flujo, L, la escala
de los bloques porosos, l, y la escala del poro, d. En este caso el volumen elemental de
referencia debe ser suficientemente grande comparado con el tamaño del bloque, para
poder obtener promedios representativos. Sin embargo, a partir de una caída de presión
en la escala del yacimiento, se generan también flujos locales entre los bloques porosos
y las fracturas como producto de la diferencia de presiones entre ellos.
Bajo la suposición de permeabilidad muy pequeña de la matriz, ocurre un proceso
transitorio mientras las presiones de las fracturas y de los bloques se igualan, ocurriendo
un flujo q, desde ellos hacia la fractura. A partir de la ecuación de continuidad en cada
uno de los medios (se incluye el flujo q y las respectivas porosidades las cuales
aumentan con la presión correspondiente) y la ecuación de momentum se puede obtener
la ecuación de flujo general en el medio poroso fracturado durante el proceso transitorio,
aplicando la ecuación de continuidad en el medio poroso y en la fractura, y la respectiva
ley de flujo
∂ p1 k ∂
( r ∂ p1
=
∂t
r ∂r
∂r
)+ η
∂ ⎡ 1 ∂ ⎛ ∂ p1 ⎞ ⎤
⎜r
⎟ .
∂t ⎢⎣ r ∂r ⎝ ∂ r ⎠ ⎥⎦
(3 .15 )
Donde
η=
k k1l 2
,
A α k2
siendo,
A=
k1
α k2
, y k= 2
μ l m2 ( β 2 + β1 )
μ l m2 ( β 2 + β 1 )
2
Siendo l, la longitud del bloque, k1 y k2 las permeabilidades de la fractura y del medio
poroso respectivamente, m2 es la porosidad de la matriz, β2 el coeficiente de
27
compresibilidad del medio, β1 el coeficiente de compresibilidad del fluido, α es una
constante adimensional que expresa la geometría de la matriz. Esta ecuación representa
la variación de la presión en la fractura p1, mientras ocurre un proceso transitorio que se
detiene cuando las presiones en la matriz y en la fractura se igualan.
En el límite, cuando l→0, η→0 sólo existe medio poroso y la ecuación de flujo en
acuíferos y yacimientos toma la forma 3.11 o 3.14, respectivamente, las cuales
corresponden al flujo de un fluido ligeramente compresible en un medio poroso.
En este caso se considera densidad y viscosidad del fluido como constantes y solamente
el fluido de los poros de la matriz es el que contribuye al soporte de la roca, es decir, la
porosidad de las fracturas no afecta a la ecuación (Barenblat, 1990). En el desarrollo de
este trabajo solamente se tendrá en cuenta el flujo a través de sistemas de fracturas.
3.4
GRUPOS DE RENORMALIZACIÓN, FUNCIONES HOMOGÉNEAS Y
EXPONENTES ANÓMALOS
Una propiedad de la naturaleza es la existencia de escalas características dentro de un
fenómeno. Por ejemplo, las fracturas en una roca se pueden presentar en escalas del
orden de milímetros hasta fracturas del orden de metros y kilómetros, o las ondas en el
océano pueden observarse en diversas escalas de longitud. Por otro lado, los fenómenos
productos de transiciones de fase en el punto crítico, pueden llegar a presentar eventos
en todas las escalas de longitud, como son las transiciones de fase líquido-gas, de
ferrogmanetos o transición de percolación que se presentan en el numeral 3.8. Estudiar
las propiedades de invarianza a la escala de esos fenómenos, permitirá reducir sus
variables y obtener soluciones autosemejantes que cumplen condiciones de
homogeneidad.
28
3.4.1 Autosemejanza, escalamiento y homogeneidad
La autosemejanza de un fenómeno es una propiedad muy importante. Un fenómeno es
auto semejante si la distribución espacial de sus propiedades, en instantes de tiempo
diferentes, se puede obtener mediante una transformación de semejanza de otra
distribución.
Por lo tanto, la distribución espacial de una propiedad permanece
geométricamente similar a sí misma en cada paso de tiempo. En este caso se dice que la
propiedad es invariante bajo una transformación auto semejante. Propiedades como la
presión del fluido en un yacimiento, el frente de una masa de fluido en el suelo o la
velocidad de propagación de una onda, permanecen invariantes en el tiempo siempre y
cuando esas propiedades sean adecuadamente escaladas. La autosemejanza también es
conocida como invarianza bajo escalamiento y es un tipo de simetría que presenta el
fenómeno ante cambios de escala.
La propiedad de autosemejanza simplifica el estudio del fenómeno pues reduce el
número de argumentos en el problema, por lo tanto, la solución de una ecuación
diferencial parcial representativa del fenómeno se reduce a la solución de un conjunto de
ecuaciones diferenciales ordinarias.
Un fenómeno auto semejante puede ser representado con una ley de escalamiento. Se
puede demostrar que la ley de escalamiento es una relación en forma de ley de potencia
de ciertas variables así:
φ ( μ ) = Aμ β
Donde A y β son constantes y β es el exponente de escalamiento.
(3.17)
Las leyes de
escalamiento expresan la propiedad de autosemejanza del fenómeno y hacen parte de las
ecuaciones o formas funcionales homogéneas que representan un fenómeno físico
(propiedad de homogeneidad). Por tanto, las funciones homogéneas se pueden expresar
en función de leyes de escalamiento. Una función ƒ (r) es homogénea si para todos los
valores del parámetro de escala λ, se cumple
29
ƒ(λr) = g(λ) ƒ(r).
(3.18)
Donde g(λ) es una función que toma la forma general (3.17) de escalamiento, g(λ)=λβ
La ecuación 3.18 significa que si se conoce el valor de la función f, en algún punto r0,
entonces es posible conocer la función en todos los puntos r, ya que cualquier número r
se puede escribir como λ r0.
Por lo tanto el valor de la función en r se obtiene
escalando el valor de la función en r0 por g(λ).
Como se puede observar, las leyes de potencia están relacionadas con la autosemejanza
bajo escalamiento de las funciones homogéneas. Los exponentes de estas leyes son
llamados también exponentes críticos en el estudio de transiciones de fase (numeral 3.8)
donde ellos definen singularidades de los parámetros de orden y exponentes anómalos
en fenómenos donde aparecen nuevas escalas a considerar. Es allí donde se puede
observar la equivalencia entre la propiedad de homogeneidad y la llamada invarianza
asintótica con respecto a un grupo de renormalización.
3.4.2 Grupos de renormalización
Los grupos de renormalización han sido utilizados en la mecánica estadística y
fenómenos críticos desde la mitad del siglo pasado. Sin embargo, fue Kenneth Wilson,
en 1979, quien los relacionó directamente con los problemas donde existen múltiples
escalas, característica muy especial de los fenómenos que ocurren en la naturaleza.
Se conoce que los exponentes en los fenómenos críticos presentan valores diferentes a
los obtenidos desde las teorías de campo medio o de Landau. La razón de esto es la
existencia de otra escala de longitud diferente a la longitud de correlación; esa escala de
longitud es la llamada escala microscópica que debe también incluirse en el análisis
dimensional que se hace sobre la función de energía libre. En tal caso se obtiene una
longitud de correlación ξ que toma la forma ξ= t
(-1/2 + α)
l
2α
, en función de la escala
microscópica l y de un exponente anómalo α. Debido a que la renormalización permite
30
introducir nuevas escalas de longitud dentro del problema, resulta ser una herramienta
muy poderosa en la solución de problemas en los cuales se han dejado de considerar
escalas importantes.
La obtención tanto de los exponentes críticos como de los
exponentes anómalos, así como la forma de sus funciones de escalamiento, se hace
escogiendo apropiadamente el grupo de renormalización ante el cual la solución del
problema es invariante.
Los problemas de flujo y transporte en rocas en estado transitorio se representan con
ecuaciones diferenciales que pueden ser solucionadas con la aplicación de grupos de
transformación o renormalización, especialmente cuando se trata de problemas que
deben incluir escalas adicionales. Lo anterior se basa en la teoría matemática del álgebra
de grupos y su aplicación a la solución de ecuaciones diferenciales como se describe a
continuación.
3.4.2.1 Grupos de renormalización de Lie
Los fenómenos de la naturaleza presentan comportamientos en diversas escalas, por lo
tanto, ella exhibe ciertas simetrías, y las ecuaciones matemáticas (ecuaciones
diferenciales) que intentan modelarla deben tomar ventaja de esas simetrías.
La
propiedad fundamental de la invarianza descrita en el numeral 3.4.1 permite realizar
transformaciones que dejan invariante la solución, es decir, permite obtener una solución
autosemejante.
Se deben entonces determinar cuáles son las transformaciones bajo las cuales el
problema es invariante y, luego, usar dichas transformaciones para derivar las
condiciones especiales desde las cuales se puede encontrar la solución o se puede
simplificar el problema. La conceptualización sobre los grupos de Lie se encuentra muy
bien fundamentada y la literatura es extensa; en este trabajo se trata de exponer los
criterios más relevantes para la aplicación de éstos a la solución de la ecuación de
difusión y lograr un mejor entendimiento de la aparición de los exponentes anómalos en
dichas soluciones.
31
3.4.2.2 Conceptos básicos.
El estudio de las ecuaciones diferenciales se puede hacer mediante la utilización de los
Grupos de Lie, los cuales fueron creados por Sophus Lie (Campos, 1995) a partir de la
teoría de grupos algebraicos pero aplicando en este caso transformaciones sistemáticas
infinitesimales.
Se le llama un grupo de transformación de Lie, o grupo de Lie uniparamétrico, a una
familia de transformaciones Tε, de un parámetro ε, que satisface la propiedad de
cercanía local, contiene la identidad y existe inversa para el parámetro, Logan (1987).
Un grupo de transformación corresponde a un cambio de variable.
Los tipos de
transformaciones simétricas bajo las cuales se estudian aquí las propiedades de
invarianza, son transformaciones del plano, Ξ2, y ellas dependen del parámetro ε.
La
representación general es la siguiente,
Tε: t* = ϕ(t, x, ε)
x* = ψ(t, x, ε).
Donde ϕ y ψ son funciones analíticas.
Una transformación particular de la familia se obtiene asignando un valor al parámetro ε.
Para observar la transformación infinitesimal y los respectivos generadores se realiza
una expansión en series de Taylor al grupo de Lie, si ésta converge es llamada serie de
Lie.
Para el tiempo:
t* = t + ε τ (t,x) + o(ε),
Para el espacio:
x*= x + ε ξ(t,x) + o(ε).
Donde τ y ξ son los generadores infinitesimales de la transformación Tε. Los dos
primeros términos de la expansión conforman la transformación infinitesimal del grupo
de Lie. El término generador indica que una transformación finita se obtiene mediante
una aplicación repetitiva de la transformación infinitesimal. El anterior es un método
para obtener las ecuaciones del grupo, si se conocen los generadores infinitesimales, o al
32
contrario, si se conocen las ecuaciones del grupo es posible obtener expresiones para los
generadores infinitesimales.
A continuación, en la Figura 3.5, se observa el esquema de aplicación de una
transformación general Tε sobre una función x=f(t) definida en el intervalo [a,b]
x = f(t)
a
x* = f*(t)
b
t
a
b
t
Figura 3. 5 Esquema de transformación de una función.
La transformación actúa sobre f. Para determinar f* se aplica el mapeo Tε al conjunto
de puntos (t, f(t)) cuando t varía entre a hasta b. Las transformaciones o grupos de Lie
más conocidos son los siguientes:
Grupo de traslación:
t*= t + ε,
x* = x.
Grupo de rotación:
t* = t cosε - x senε,
x* = tsenε+xcosε.
Grupo de estiramiento:
t* = eaε t,
x* = ebε x.
Para cada una de las transformaciones anteriores se definen las expresiones para
designar a cada uno de los generadores de acuerdo a las operaciones algebraicas
requeridas, Campos (1995).
La importancia de los generadores infinitesimales radica
en que son utilizados para la solución de las ecuaciones diferenciales.
Sin la
consideración de los generadores no es posible identificar la simetría completa de una
ecuación diferencial parcial.
33
3.4.2.3 Aplicación de los Grupos de Lie a la ecuación de difusión.
Generalmente, cuando en un problema se reduce el número de variables de una ecuación
diferencial parcial (EDP) a partir de las condiciones de simetría, los invariantes del
grupo se convierten en las nuevas variables.
La variable dependiente ahora será la
variable de semejanza obtenida y al ser reemplazada la nueva forma funcional en la
ecuación diferencial parcial, se obtiene una ecuación diferencial ordinaria (EDO)
fácilmente solucionable. La solución de semejanza obtenida así satisface una ecuación
diferencial parcial auxiliar cuyos coeficientes dependen de los generadores
infinitesimales del grupo. La forma funcional de la solución autosemejante se obtiene
solucionando la ecuación característica anterior.
Bluman y Cole (1974) construyeron la solución de semejanza general de la ecuación de
calor, la cual considerando coeficiente de difusividad unitario se escribe
∂u ∂ 2u
=
∂t ∂x 2
(3.19)
Se trata de encontrar el grupo de transformación que deja invariante la anterior EDP y
los respectivos generadores del grupo, en este caso: η, ξ, τ. Inicialmente se considera un
grupo de Lie uniparamétrico así,
u * = U * ( x, t , u ; ε )
x * = X * ( x, t , u ; ε )
t * = T * ( x, t , u, ε ).
(3.20)
Este grupo deja invariante a S , sistema que incluye la EDP y sus condiciones.
Si v= U* (x, t, θ (x, t ); ε ) satisface a S* cuando u = θ(x,t) es una solución de S.
El sistema S* se obtiene reemplazando u por v en el sistema S. Expandiendo la ecuación
3.20 se obtiene la siguiente forma infinitesimal,
34
x* = x + ε ξ(t,x) + o(ε2)
t* = t + ε τ (t,x) + o(ε2)
(3.21)
u*=u + εη(t,x) + o(ε2),
a partir de la cual se obtiene la ecuación funcional para θ,
θ(x+ ε ξ+ o(ε2), t + ε τ + o(ε2))= θ(x,t)+ εη(t,x, θ)+ o(ε2).
(3.22)
A partir de los términos o(ε) de la forma anterior se obtiene la siguiente EDP que
satisface θ(x,t) dado los generadores η, ξ, τ.
ξ(x,t, θ) θx+ τ(x,t, θ) θt = η(x,t, θ)
(3.23)
Esta ecuación es llamada condición de superficie invariante y su solución se obtiene
utilizando la ecuación característica. La solución de la ecuación 3.23 toma la forma
θ = f(ζ) ƒ(x,t), donde f es una función arbitraria de ζ , siendo ζ y ƒ, funciones
conocidas de x , t.
Si u=θ (x,t) es alguna solución de S, entonces u*=U*(x,t,u; ε) es la solución
correspondiente a S*. A partir de esta solución se obtienen los generadores ηx y ηt para
la primera derivada de u* y para la segunda derivada se obtiene ηxx y ηtt.
Así la
solución invariante u = θ (x,t) satisface:
u ∗xx − ut∗ =θ xx − θ t + ε (η xx − ηt ) + o(ε )
(3.24)
Solucionando e igualando los coeficientes de θ a cero se obtiene el siguiente grupo
35
ξ = κ + δ t + β x + γ xt
τ = α + 2β t + γ t 2
(3.25)
⎛ x2 t ⎞ δ x
+ ⎟⎟ −
+ λ.
⎝ 4 2⎠ 2
η = −γ ⎜⎜
Donde α, β, γ, δ, κ, λ son parámetros arbitrarios. Todos los parámetros, excepto δ,
representan transformaciones triviales del espacio (x,t), así κ representa translación en x;
α representa translación en t; δ significa invarianza galileana; β denota invarianza de
estiramiento y γ representa la invarianza bajo el grupo de un parámetro de
transformaciones proyectivas (homotecia). El grupo (3.25) de seis parámetros es un
grupo de Lie de Transformación que deja invariante la ecuación 3.19. Estos parámetros
son considerados seis elementos que generan un álgebra de Lie de dimensión finita. El
excelente texto de Campos, 1995, los define como v1=κ, v2=α, v3=γ, v4=β, v5=δ, v6=λ,
y llama a la transformación de estiramiento, "escalonamiento".
3.4.2.4 Solución de la ecuación de flujo en medio poroso a partir de Grupos de Lie.
La ecuación de flujo en medio poroso tiene la forma de una ecuación de difusión, por lo
tanto, se utiliza el mismo procedimiento para obtener su solución. Las soluciones de
autosemejanza se obtienen desde las propiedades de invarianza de las EDP bajo
transformaciones de Grupos de Lie. La ecuación de difusión se soluciona considerando
invarianza bajo la transformación particular de estiramiento, lo que significa que las
nuevas variables dependientes e independientes son múltiplos de las anteriores, por lo
tanto la nueva EDP se obtiene en un nuevo sistema de coordenadas que es igual por
algún múltiplo al anterior sistema de coordenadas.
A partir del grupo general (3.25) se toman los siguientes parámetros: α=γ=δ=κ=0. Este
grupo debe dejar invariantes los ejes x=0 y t=0. Por lo tanto, el nuevo sistema se
obtiene tomando estos valores y reemplazando (3.25) en (3.21),
x*= x + ε (β,x)
t* = t + ε (2βt)
u*=u + ε (λu)
36
(3.26)
Ahora se trata de buscar la solución de similaridad, la cual se obtiene desde la invarianza
del grupo, a partir de las siguientes ecuaciones características del grupo
dx
dt
du
=
=
.
βx 2 β t λ u
(3.27)
La solución a esta ecuación diferencial de primer orden nos da la solución de semejanza
ζ =
x
,
2 t1 / 2
(3.28)
la cual lleva a la forma general de la solución de la ecuación lineal de difusión
u ( x, t ) = t 2 / 3 f (ζ ) .
(3.29)
A partir de esta forma general obtenida desde los Grupos de Lie de la ecuación de
difusión, se estudian todas las ecuaciones relacionadas con problemas de difusión.
Los problemas de difusión no lineal (tipo ecuación de filtración) que se originan en la
dinámica de los fluidos en medios porosos, son considerados como fenómenos lejos del
equilibrio, es decir, son procesos dinámicos cuyas soluciones se buscan en el régimen
asintótico intermedio, donde presentan propiedades especiales.
En le caso de difusión no lineal, el coeficiente o parámetro de flujo no es constante sino
que depende del valor de ∂u / ∂t, el cual presenta un cambio de signo dando lugar al
planteamiento de un problema de coeficiente discontinuo.
La explicación de este
comportamiento se mostrará en detalle en el numeral 3.7.5.
La forma general de la
ecuación es entonces
∂u
∂ 2u
= k ( w) 2 2 .
∂t
∂ x
37
(3.30)
La forma de similaridad de esta ecuación igualmente se puede encontrar desde la teoría
de grupos. En este caso ∂u/∂t = 0 sobre ζ=ζc. Haciendo λ / 2β=α(ω), la forma de la
solución a la ecuación no lineal de difusión se puede escribir
u ( x, t ) = t α ( w ) f (ζ , w),
(3.31)
donde α(w) es considerado un exponente anómalo. Por lo tanto el nuevo grupo se puede
reescribir como
x*= x + ε x
t* = t + ε (2t)
(3.32)
u*= u + ε (2α)u.
La idea de encontrar la aproximación para el llamado exponente anómalo α, es que el
comportamiento asintótico de la expansión de perturbación de la solución original,
aunque no es uniforme, incluye una aproximación por perturbación de las propiedades
del grupo.
A partir de estas propiedades se puede obtener una expresión para el
exponente α, Goldenfeld (1992).
3.5 ANÁLISIS DIMENSIONAL Y SOLUCIONES ASINTÓTICAS INTERMEDIAS
Se conoce que las funciones que expresan leyes físicas deben cumplir la propiedad
fundamental de homogeneidad generalizada, la cual permite reducir el número de
argumentos del problema. Esto significa que el sistema tiene algún tipo de simetría, por
lo tanto, el problema debe presentar invarianza ante un grupo de Lie como se estudió en
el numeral anterior.
El Teorema Π de Buckingham, Logan (1987), sobre el cual se basa el análisis
dimensional permite expresar una función de n=k+m variables en una relación entre
cantidades dimensionales que caractericen el fenómeno: a = f (a1,...., ak, b1,.....bm ,ε ).
donde k es el número de parámetros independientes del problema con dimensiones
independientes y m el número de parámetros con dimensiones dependientes.
38
El
fenómeno también depende de un parámetro adimensional ε..
Los productos
adimensionales se forman con los parámetros que tienen dimensiones dependientes y la
forma funcional tendrá un número menor de argumentos
Π = φ (Π1 ,....Π m , ε )
(3.33)
La función φ posee la propiedad de homogeneidad que puede expresarse con la forma
general 3.17, así la transformación de los parámetros con dimensiones independientes,
se puede hallar haciendo simplemente un cambio de sistema de unidades. Los demás
parámetros con dimensiones dependientes variarán de acuerdo a sus dimensiones en la
siguiente forma
b11 = b1
b21 = b2 , ......, bm = Bbm
a11 = a1
a12 = a2 , ......, ak = ak
(3.34)
La anterior transformación forma un grupo simple de estiramiento o grupo de Lie de la
forma (3.25). B es un número positivo arbitrario llamado el parámetro del grupo. Los
productos adimensionales Π, Π1,.....Πm, , no cambian para esa transformación, es decir,
son invariantes respecto a este grupo. El Teorema Π es consecuencia del principio de
covarianza y las relaciones existentes con significado físico en un problema son
invariantes con respecto a un grupo de transformación simple.
Los problema en que estamos interesados se estudian en el régimen de asíntotas
intermedias, donde se obtienen soluciones válidas para tiempos y distancias en los cuales
los detalles de las condiciones iniciales y de borde no tienen influencia, sin embargo, el
sistema se mantiene aún lejos del equilibrio final. Las soluciones a estos problemas son
llamadas asintóticas.
En todo modelo de un fenómeno físico se desprecian aquellos factores que son
considerados no esenciales (ejemplo, una escala del problema) y hacen parte de un
producto adimensional de valor pequeño que en el problema es despreciable, de esta
manera la función φ posee un número menor de argumentos. Si en este caso la función φ
39
tiene un límite finito diferente de cero, existe solución al modelo propuesto y se dice que
el problema posee autosemejanza de primer orden; la solución es invariante respecto a
un grupo simple de estiramiento.
Así, para la forma funcional 3.33 se obtiene solución autosimilar de primer orden
respecto al parámetro Πm,
La función φ puede ser reemplazada por su límite finito
cuando Πm, va a cero, esto para ε=o, en cuyo caso bm desaparece de la relación y
también, Πm.
Sin embargo, existen problemas en los cuales el parámetro bm, no desaparece en el
régimen de asíntotas intermedias, por lo tanto, el límite de una función de la forma
(3.33) cuando Πm va a cero no existe, es decir no se encuentra solución de esta forma.
La función φ puede tener solución aún para valores pequeños de Πm si en lugar de
considerar la forma anterior (3.33) se considera una nueva forma funcional.
Π =φ (
Π1
Π
,....., αmm−−11 , ε )
α1
Πm
Πm
(3.35)
Si en este caso la función φ tiene un límite finito diferente de cero existe solución al
nuevo modelo propuesto y se dice que el problema posee autosemejanza de segundo
orden, en ese caso la solución es invariante respecto a un grupo de Lie más complejo
(β= α ≠1 en (3.17)) que puede tomar la siguiente forma
a ' = Bα m a
b1' = Bα1 b1
b2' = Bα 2 b2 ,
a1' = a1
a2' = a2
, bm' = Bbm
ak' = ak
.
(3.36)
La función φ, que tiene solución asintótica y permanece invariante bajo la siguiente
transformación
Π = Bα m Π ,
Π11 = Bα1 Π1
40
Π12 = Bα 2 Π 2
Π1m = B Π m ,
(3.37)
los exponentes del nuevo grupo α1, α2,.... αm dependen del parámetro ε y se les conoce
con el nombre de exponentes anómalos.
Aunque las formas funcionales (3.33) y (3.35) poseen la propiedad de homogeneidad
para la función φ, existe una diferencia esencial entre las dos formas. En el primer caso
la solución del problema se obtienen mediante un procedimiento simple del análisis
dimensional (asociado a una transformación de estiramiento simple).
En el segundo
caso la homogeneidad de la función debe expresar una propiedad especial del problema.
El exponente anómalo α y la solución completa del problema no se puede obtener desde
el análisis dimensional sino que debe acudirse a otras herramientas para obtener la
solución completa.
En la literatura se han planteado diversas herramientas para la obtención de los
exponentes anómalos, basadas en los conceptos fundamentales de los grupos de
renormalización. Estos exponentes se han obtenido para complementar la solución al
problema no lineal de difusión de una masa de fluido en un suelo no saturado. Cole y
Wagner (1996) realizaron una expansión de perturbación al grupo simple de Lie.
Goldenfeld (1992) realizó una expansión de perturbación a la solución del problema
simple de difusión y renormalización de la masa de fluido. Baremblatt (1987) obtuvo el
exponente anómalo planteando el problema de autovalor o utilizando leyes de
conservación no integrables (Barenblat, 1996). Gómez y Mesa (2000) desarrollan un
algoritmo numérico para resolver el problema de autovalor y estiman tiempos de
propagación y alcance de masa de fluidos.
La aproximación de grupos de renormalización se basa en las ideas de Grupos de Lie y
han sido utilizados en diferentes áreas de la física por algunos autores como Kadanoff
(1966) y Wilson (1979), entre otros, quienes utilizaron los grupos de renormalización
para el estudio de los exponentes críticos en transiciones de fase. Kolmogorov (1941)
aplicó estos conceptos al problema de turbulencia.
41
3.6 LEYES DE CONSERVACIÓN DE MASA
En física es muy importante el estudio de las operaciones ante las cuales los fenómenos
físicos permanecen invariantes (traslación en el espacio, traslación en el tiempo,
rotaciones y estiramientos), las cuales son consideradas por los grupos de
renormalización presentados en el numeral 3.4.2.1.
Existe una relación entre leyes de conservación y simetrías de las leyes físicas: Si las
leyes son simétricas para traslación en el espacio, el momentum se conserva; si las leyes
son simétricas frente a una traslación en el tiempo, la energía se conserva. Si existe
invarianza a una rotación de un ángulo en el espacio, se habla de conservación de
momentum angular. Como se observa la invarianza o simetría conducen a leyes de
conservación.
Barenblat (1996) estudió la ley de conservación de masa para el problema de difusión de
una masa de fluido en un medio poroso. Para ello realiza una perturbación en el
parámetro pequeño ε=(k1/k2)-1.
La perturbación se hace a partir del problema
idealizado, es decir para ε=0. Hasta un término de orden uno, como se observa a
continuación,
φ (ς , ε ) =
1
(1 − ς 2 ) + o(ε )
16
ς0 =
1
+ o(ε ) .
2
(3.38)
En este punto, Barenblat acude a la información cuantitativa que puede ofrecer el
problema, como es una ley de conservación de la masa de agua. Cuando no existe
retención residual en el medio poroso (k1=k2), la masa de agua se mantiene constante
durante toda la evolución del problema, por lo tanto es aplicable una ley de conservación
que es valida en todas las etapas del movimiento. Ésta se obtiene integrando la ecuación
diferencial de flujo.
∂ r1
rh(r , t ) = cte
∂t ∫0
42
(3.39)
La aplicación del análisis dimensional permite involucrar otras variables o parámetros
que afectan a la variable dependiente del problema (la cabeza hidráulica). En el caso de
flujo en acuíferos puede involucrarse otro valor de conductividad y dar origen al
problema de coeficiente discontinuo. Aparece entonces, un nuevo parámetro del
problema, λ=k1 / k2 que debe incluirse en el análisis dimensional de la ecuación. La
consideración de este nuevo parámetro, producto de una física especial del problema, da
lugar a que una variable que siempre ha sido ignorada en la solución tradicional no
pueda omitirse: el radio del pozo, rp.
En este caso, el radio del pozo es una escala
adicional que aparece en el problema, tiene dimensión de longitud y entra a conformar
otro grupo adimensional. Este nuevo grupo hace perder la semejanza de la solución,
puesto que la solución ya no se podrá expresar en términos de una sola variable, sino de
dos variables, la distancia radial y el radio del pozo.
Ahora se trata de definir si ϕ, la función de la forma adimensional con los nuevos
argumentos tiene un límite finito diferente a cero cuando se analiza su comportamiento
en una etapa intermedia de tiempo, es decir, cuando el comportamiento no es ni afectado
por las condiciones iniciales, ni por lo que ocurre en el infinito. Si los nuevos
argumentos pueden simplemente ser ignorados en el problema, entonces el análisis
dimensional es suficiente para obtener la solución y se dice que existe autosimilaridad
completa (de primer orden) respecto a los nuevos parámetros. Por lo tanto, la solución
en este caso ignora el radio del pozo y se tiene la solución equivalente a una fuente
puntual instantánea, como es la que corresponde a la solución presentada por Theis para
flujo en acuíferos o la ecuación integral exponencial en yacimientos.
Si el problema es de coeficiente discontinuo y se ignora la nueva escala (radio del pozo)
se llega a una contradicción en el desarrollo matemático, por lo tanto resultará incorrecta
la suposición de autosemejanza completa respecto a ese nuevo parámetro.
Así, rp
aparecerá en la solución elevado a un exponente anómalo y debe acudirse entonces a
otras técnicas como son la solución de un problema de autovalor o el análisis de un
grupo de renormalización más complejo para obtener el exponente.
43
3.7. NO LINEALIDAD DEL FLUJO EN YACIMIENTOS Y ACUÍFEROS
La realización e interpretación de pruebas hidráulicas en yacimientos y acuíferos deben
considerar las condiciones elásticas y no elásticas de las formaciones geológicas
almacenadoras de fluidos, puesto que tales condiciones dan lugar a la linealidad o no
linealidad de las ecuaciones de flujo, lo que requiere un manejo diferente del problema
desde el punto de vista físico y matemático en cada uno de los casos.
En el primer caso, condición elástica, la solución a la ecuación de flujo se obtiene
aplicando el análisis dimensional y hay autosemejanza de primer orden. En el segundo
caso, condición no elástica, la solución no se puede obtener con el análisis dimensional
solamente, sino que debe acudirse a otra herramienta que permita obtener los exponentes
anómalos a los que da lugar el fenómeno físico. En este último caso se debe aplicar una
suposición de autosemejanza de segundo orden, la cual se encuentra asociada a un
grupo especial de renormalización para obtener la solución de acuerdo a los criterios
expuestos en el numeral 3.5.
3.7.1
Flujo en medios elásticos
Cuando se habla de medios elásticos, se trata de aquellas formaciones geológicas donde
es posible considerar que la porosidad varía linealmente con la presión. Además el
material de la formación se deforma elásticamente, por lo tanto, un proceso de carga es
idéntico a un proceso de descarga: los procesos son reversibles.
En este caso, las ecuaciones de flujo se obtienen considerando un fluido ligeramente
compresible y una formación geológica de gran espesor que presenta un comportamiento
elástico. El coeficiente de difusividad depende de las propiedades físicas del fluido y de
la roca.
Para obtener la expresión del coeficiente de difusividad se considera la
variación lineal de la densidad del fluido con la presión p, y la variación lineal de la
porosidad con el esfuerzo σ (primer invariante del tensor de esfuerzos) que se expresan,
44
ρ
= 1 + β f ( p − p0 )
ρ0
(3.40)
m
= 1 + β r (σ − σ 0 )
m0
Donde, ρo y mo. son valores de la densidad y porosidad de referencia; po y σo son
valores de la presión y esfuerzo de referencia; βf es el coeficiente de compresibilidad
del fluido y βr es el coeficiente de compresibilidad de la roca.
Estas expresiones se reemplazan en la ecuación de continuidad y en la ecuación de
momentum, y se obtiene la expresión para el coeficiente de difusividad
K=
k
μ (m0 β f + β r )
,
(3.41)
La ecuación de flujo en un medio elástico es la ecuación lineal clásica análoga a la
ecuación de conducción de calor.
∂p
∂2 p
=K 2 ,
∂x
∂t
(3.42)
que presenta un coeficiente constante de difusividad el cual expresa las condiciones del
fluido y del medio.
3.7.2
Flujo en medios elastoplásticos
Los estados de esfuerzos actuando sobre los materiales presentes en la naturaleza pueden
llevar a un proceso de deformación irreversible. En este caso, un proceso de carga de la
roca difiere de un proceso de descarga. El material responde en forma diferente, por lo
tanto, la variación de la porosidad dependerá de si la formación geológica se encuentra
en proceso de carga o en proceso de descarga.
45
En procesos donde el esfuerzo efectivo σ aumenta (la presión del fluido disminuye
debido a que el estado de esfuerzos totales del sistema permanece constante) la
porosidad varía según un coeficiente de compresibilidad de la roca βr1 .
En procesos donde el esfuerzo efectivo σ disminuye (la presión del fluido aumenta), la
porosidad varía con un coeficiente de compresibilidad βr2
.
Los coeficientes βr1 y βr2
toman valores diferentes en los dos procesos y por lo tanto, el coeficiente de difusividad
varía para los dos tipos de procesos (Barenblat, 1990).
3.7.2.1 Variaciones de presión a partir de perturbaciones en los yacimientos
Los descensos de presión (caída de presión) en un yacimiento se presentan durante los
procesos de explotación (bombeo o drenaje), fenómenos en los cuales el fluido sale de la
matriz de la roca. Los ascensos de presión se presentan cuando existe una recarga hacia
el yacimiento, ya sea natural o debida a un proceso de inyección, o simplemente cuando
después de un proceso de explotación, se permite la recuperación de los niveles del
fluido. Ambos procesos deben ser considerados como procesos diferentes desde el
punto de vista de las condiciones elásticas de la roca que almacena el fluido. En ambos
casos, la roca y el fluido, comparten el estado de esfuerzos actuando en el yacimiento, el
ascenso o caída de presión ocasiona que el fluido asuma parte de ese estado de esfuerzos
o, por el contrario, que ceda parte de ellos como se describe a continuación.
3.7.2.2 Proceso de caída de presión.
Al disminuir la presión del fluido en un proceso de extracción, la armazón de la roca
soportará mayor peso de las capas superiores y se comprimirá disminuyendo la
porosidad, por lo tanto, por cada descenso unitario de la presión la columna de acuífero
se comprimirá un volumen S2 y cederá ese mismo volumen. En este caso se dice que
el estrato se encuentra bajo carga y el coeficiente de difusividad se expresa como
K2 =
k/μ
,
mβ f + β r 2
(3.43)
donde βr2 es el coeficiente de compresibilidad del medio poros en descarga cuando
existe disminución de presión.
46
3.7.2.3 Proceso de aumento de presión.
En un proceso de inyección, el fluido entra a presión al yacimiento desde el pozo. Al
aumentar la presión del fluido la armazón de la roca soporta menos peso debido a que
parte de los esfuerzos son asumidos por el fluido, la estructura de la roca se expande un
volumen S1 , y ella asume ese mismo volumen de fluido. En este caso se dice que el
estrato se descarga y el coeficiente de difusividad se expresa como
K1 =
k/μ
,
mβ f + β r1
(3.44)
donde βr1 es el coeficiente de compresibilidad del medio poroso cuando existe aumento
de presión; k/μ es la conductividad hidráulica y expresa las condiciones del medio y del
fluido, y tiene un valor igual para ambos casos. La conductividad hidráulica es análoga
a la conductividad térmica.
Los denominadores de las ecuaciones (3.43) y (3.44)
expresan las condiciones de compresibilidad y por tanto de almacenamiento de la roca.
Esta expresión es el análogo al calor específico en la ecuación de conducción de calor.
3.7.3
Coeficiente de difusividad en acuíferos
Los acuíferos son formaciones geológicas que pueden presentar condiciones
elastoplásticas.
Las diversas pruebas hidráulicas o perturbaciones sobre los acuíferos
buscan la obtención de sus parámetros hidráulicos. Los coeficientes de difusividad en
acuíferos dependen de si el proceso o la perturbación corresponde a un aumento de la
presión o si corresponde a una caída de presión y se expresan como
K1 =
T
KB
=
S1
S1
T
KB
=
K2 =
S2
S2
si
∂h
<0
∂t
∂h
>0
si
∂t
(3.45)
donde los parámetros para cada proceso son T, transmisividad (L2 T-1), K, conductividad
hidráulica (L1 T-1), B, espesor del estrato (L), S, coeficiente de almacenamiento
(adimensional). Así, las unidades de la difusividad son (L2 T-1).
47
La conductividad hidráulica es un parámetro tradicionalmente utilizado en acuíferos y se
expresa como: K=k ρ g /μ., donde k es la permeabilidad intrínseca del medio (L2 ), ρ, la
densidad (m/L3), g, gravedad (m2/s) y μ la viscosidad dinámica del fluido (m/LT). En el
caso de acuíferos la variable dependiente es la llamada cabeza hidráulica (energía por
unidad de peso).
El coeficiente de almacenamiento entonces expresa las condiciones de compresibilidad y
variación de la porosidad así: S = mβ +βr, por lo tanto es el parámetro que cambia
cuando el acuífero es considerado elastoplástico y la perturbación realizada genera
ascensos y descensos de presión.
3.7.4
Leyes de autosemejanza en el flujo en pozos
Las ecuaciones que rigen el flujo en yacimientos y acuíferos se solucionan a partir de la
suposición de autosemejanza de acuerdo a criterios de Grupo de Lie presentados en el
numeral 3.4.2.4. Las condiciones físicas del problema presentadas en el numeral 3.7.2
dan lugar a la utilización de uno u otro grupo de renormalización. Aquí se presentan los
dos tipos de autosemejanza vistos: el primero relacionado con un grupo de Lie simple,
de la forma (3.26), y el segundo relacionado con un grupo de Lie más complejo, de la
forma (3.32).
3.7.4.1 Ley de autosemejanza de primer orden
Tradicionalmente se ha considerado que las rocas que almacenan fluidos poseen un
comportamiento elástico, cuyo coeficiente de compresibilidad del medio poroso, βr es
idéntico para cualquier proceso. Por tanto, la porosidad y el coeficiente de
almacenamiento también se mantienen constantes, según lo expresado en el numeral
3.7.3. Para medios porosos en general, se considera que la solución del flujo se obtiene
basados en la autosemejanza de primer orden que se presenta a continuación. Para
obtener la solución autosimilar de la ecuación de flujo (3.11) se requieren las siguientes
condiciones de borde:
48
C.I : h(r ,0) = h0
lim (r
C.B : h(∞, t ) = h0 ,
∂h
∂r
)=
Q .
2π T
(3.46)
r →0
En este caso, la solución es la cabeza hidráulica o energía por unidad de peso (L) que
depende de la distancia desde el pozo, r (L); el tiempo (T); Q, el caudal de explotación
(L3/T) y T, la transmisividad (L2/T). La forma funcional es: h= f (r, t, Q, T). El análisis
dimensional permite obtener dos productos adimensionales y expresar la
forma
funcional así:
r2
π 2 = = u,
Tt
h
π1 =
,
QT −1
h=
Q
f (u ),
T
(3.47)
la función f(u) tiende a un límite finito y satisface la siguiente ecuación diferencial
ordinaria
f " (u ) + (1 / u + 1 / 4) f ' (u ) = 0
(3.48)
1
∂f
=
∂u 4π .
u→∞
(3.49)
con condiciones,
lim u
y solución
ho − h =
Q
4π T
∞
r2 S
e −u
.
u
,
u
∂
=
∫u u
4Tt
(3.50)
En el proceso de bombeo o producción existe un frente de descarga en el cual ocurre
descenso de presión en el acuífero y la solución a la ecuación diferencial ordinaria
utilizando argumentos dimensionales es la misma llamada solución de Theis en
acuíferos o solución Ei en yacimientos.
Si el problema se estudia en medios porosos elastoplásticos las ecuaciones de flujo
consideran coeficientes de difusividad diferentes K1 y K2 con las mismas dimensiones y
no se incluyen constantes adicionales con dimensiones independientes.
El nuevo
argumento dimensional se conforma con los coeficientes de difusividad, λ= K1 / K2.
49
Si se soluciona el problema de Cauchy (condición inicial) en los dos casos anteriores se
puede observar la influencia de las condiciones iniciales durante el proceso,
especialmente para tiempos grandes. Para el segundo caso la solución al problema de
Cauchy utilizando análisis dimensional tendría la siguiente forma
Q
h=
f (u, ε , λ ),
K1
r2
u=
,
K1 t
ε=
l
,
K1 t
λ=
K1 .
K2
(3.51)
En este caso existe un argumento dimensional más, el argumento ε, que incluye la escala
de la condición inicial. Generalmente, se asume que la escala l es muy pequeña y por lo
tanto despreciable, lo cual llevaría a la solución de la forma (3.26). En ambos casos la
solución a un problema de Cauchy olvida los detalles de los datos iniciales y tiende a la
autosemejanza de primer orden. Esto da origen a la validez de la solución tipo fuente
puntual instantánea para la ecuación de difusión clásica presentada en el numeral
3.2.3.2.
En los dos casos anteriores existe una ley integral para el volumen de fluido involucrado
en el problema,
∫
∞
−∞
(h − h0 )dr = Q,
h(r ,0) = h0 .
(3.52)
En este caso se considera que el fluido involucrado en el proceso de bombeo o
disminución de presión es igual al volumen de fluido involucrado en el proceso de
aumento de presión, ya sea por recuperación de presión después de terminar el bombeo o
por un proceso de inyección.
3.7.4.2 Ley de semejanza de segundo orden
Un nuevo enfoque al problema de abatimiento y recuperación de niveles o presión en
yacimientos elastoplásticos plantea los dos procesos de descenso y ascenso de presión
como procesos diferentes. Ellos presentan dos frentes diferentes de descenso y ascenso
de presión que son llamados frente de carga y frente de descarga, en los cuales se
50
involucran parámetros de la roca que cambian según sea el proceso. Para que su
solución tenga en cuenta la diferencia de los procesos debe plantearse una ley de
semejanza de segundo orden, lo que implica un grupo de renormalización especial.
Como se estudió en el numeral 3.6, la ley de semejanza de segundo orden, se origina
cuando no existe una ley de conservación de masa constante, es decir, cuando la masa
de fluido es variable y puede ser normalizada. Se plantea entonces una ley de semejanza
de segundo orden para este problema con una condición inicial que incluye una nueva
escala, la cual ahora será considerada dentro de los argumentos adimensionales; la nueva
escala no desaparece en el régimen de asíntotas intermedias debido a que la forma
funcional es del tipo de ecuación (3.35). En este caso, el grupo de renormalización
asociado a la solución es un grupo más complejo que involucra un exponente anómalo.
La presencia de la nueva escala en la solución hace que la función φ no vaya a un límite
finito cuando la escala l, va a cero, entonces el límite de la función crece asintóticamente
hacia el infinito dependiendo del valor que tome λ.
3.7.5
El problema con coeficiente discontinuo en pruebas de pozos
Existen diversos procesos o perturbaciones que realizados sobre yacimientos
elastoplásticos dan lugar a una componente no lineal y por tanto se consideran como
problemas de coeficiente discontinuo.
El efecto de las condiciones elastoplásticas de las rocas que almacenan fluidos es
fundamental para la solución de las ecuaciones de flujo y posterior interpretación de las
pruebas hidráulicas. Existen dos ejemplos claros de esta situación en el caso de pruebas
hidráulicas: el primero es debido a una perturbación periódica en un pozo, y el segundo
es debido a la explotación y posterior recuperación de niveles en el yacimiento. El
ejemplo del bombeo periódico tiene como propósito ilustrar un problema de coeficiente
discontinuo con condiciones matemáticas manejables.
51
Desde el punto de vista práctico, en explotaciones petroleras la inyección de agua en los
pozos tiene utilidad demostrada. A continuación se expone un problema común en la
explotación de yacimientos o acuíferos, como es el bombeo periódico desde un pozo.
3.7.5.1 Bombeo periódico en pozos
En la práctica un pozo se trabaja a una tasa variable de bombeo, debido a que las
condiciones técnicas de funcionamiento difícilmente pueden ser mantenidas constantes
durante el tiempo que demore el proceso. Cuando se introduce una perturbación en un
yacimiento como es el bombeo a una tasa de producción usualmente periódica, q(t)=qo
senwt o q(t)=qo coswt, se obtiene un frente también periódico que ha sido observado en
sitios cercanos al pozo donde la respuesta puede ser sincrónica con la perturbación.
Lejos del pozo también se observa una respuesta semejante a la anterior, pero algo más
amortiguada y rezagada.
Se requiere, por lo tanto, estudiar el frente de carga y descarga del yacimiento, el cual
presenta un régimen continuo y periódico de ascensos y descensos de la presión debido
al efecto de bombeo periódico. Es importante definir y localizar el frente de descarga, es
decir, la coordenada x=μ(t) que es periódica y depende del tiempo. Este frente define
los dos tipos de procesos, el ascenso y descenso de presión.
Las ecuaciones de flujo en acuíferos en coordenadas cilíndricas considerando coeficiente
discontinuo son las siguientes,
⎛ 1 ∂h ∂ 2 h ⎞
∂h
= K 1 ⎜⎜
+ 2 ⎟⎟
∂t
∂
r
r
∂r ⎠
⎝
⎛ 1 ∂h ∂ 2 h ⎞
∂h
= K 2 ⎜⎜
+ 2 ⎟⎟
∂t
⎝ r ∂r ∂r ⎠
si
∂h
≤ 0,
∂t
(t < π )
si
∂h
≥ 0,
∂t
(π < t < 2π )
Con condiciones iniciales h(r, 0)=ho, y condiciones de borde
h(∝,t)=ho y lim (r
Q cos wt
∂h
) r = rp =
∂r
2π K 1
52
(3.53)
Por lo tanto la variable dependiente se expresa como h = f (r, t , Q , K1 , K2 , rp ).
La forma funcional y los productos adimensionales se expresan como,
rp
K
r2
u=
,, ε =
, λ= 1
K 1t
K 1t
K2
Q cos wt
h=
f (u, ε , λ ),
K1
(3.54)
Donde w es la frecuencia angular y Q la amplitud del caudal de bombeo. Las demás
variables de los argumentos adimensionales han sido identificadas en el numeral 3.7.4.
En la figura 3.5 se puede observar la forma de las oscilaciones del caudal de bombeo,
donde el período T=2π /w, la frecuencia w es igual a 1 y por facilidad se toma un caudal
unitario.
Si se acepta que las oscilaciones de presión son sincrónicas con las oscilaciones del
caudal de bombeo en el pozo, el cambio de fase (cambio entre el descenso y el ascenso
de presión) dependerá del tiempo y estará definido por ro(t) que se expresa a partir del
parámetro adimensional u,
r0 (t ) = u 0 K t 1 / 2
(3.55)
ro(t) define el cambio de signo para ∂h/∂t como se expresa en la ecuación (3.31). Los
cambios de ascensos y descensos de presión se definen así:
Descenso de presión para t<π
Ascenso de presión para π<t<2π donde t, el tiempo se define en segundos.
1.5
Qcoswt
1
0.5
ro
0
-0.5 0
1
2
3
4
5
6
7
8
-1
-1.5
Tiempo, s
Figura 3.5. Oscilaciones de presión de bombeo periódico, forma cosenoidal.
53
Los productos adimensionales u y ε planteados en 3.32, definen las escalas espaciales
del problema y el rango de tiempo en el cual se consideran asíntotas intermedias.
Tradicionalmente la solución a los problemas con caudal de bombeo variable en
yacimientos elásticos se ha obtenido utilizando el principio de superposición lo cual
conlleva a la convolución y posterior aplicación de Transformadas de Laplace. La
solución tipo fuente lineal considera yacimientos elásticos. En este caso, la distribución
de presiones se encuentra independiente del radio del pozo.
3.7.5.2 Ley de autosemejanza de segundo orden
Para que la solución tenga en cuenta la diferencia de los procesos, propuestos en el
numeral anterior, debe plantearse una ley de semejanza de segundo orden. Esta ley se
origina cuando no existe una ley de conservación de masa constante, es decir cuando la
masa de fluido es variable durante el proceso y puede ser normalizada.
Se plantea entonces esta ley con una condición inicial que incluye la escala del radio del
pozo, la cual se considera en un nuevo argumento adimensional. La nueva escala no
desaparece debido a que la forma funcional se puede expresar en función de un
exponente anómalo. Esto hace que la solución se base en un grupo de renormalización
más complejo que aquel utilizado por la ley de autosemejanza de primer orden.
A partir de la forma adimensional expuesta en la ecuación 3.54, se desarrolla el
procedimiento regular para solucionar la ecuación 3.53 vía la suposición de
autosemejanza de segundo orden; para ello se propone un nuevo argumento
adimensional que incluye el radio del pozo.
El hecho de que esta escala no sea despreciable durante todo el proceso hace que la
función f de la ecuación 3.54, tienda a un límite infinito. Debe existir un exponente α
que haga finito el límite de la función cuando rp es muy pequeño, o el tiempo del
proceso es grande. La solución límite podría tomar la siguiente forma
54
h0 − h =
β (Q cos wt ) rpα
K 11+α t α
f (u , λ ) .
(3.56)
Al reemplazar esta forma de solución en la ecuación 3.53, se debe obtener un par de
ecuaciones diferenciales ordinarias, ecuaciones que deben ser diferenciables y continuas
para todas las coordenadas ro de avance del frente de carga y descarga.
Para hallar el exponente anómalo α debe plantearse la ley de conservación no integral ,
cuyo valor debe ser una función de los dos parámetros de flujo, K1 y K2, y la derivada
∂h/∂r evaluada en r = ro.
Como se dijo anteriormente una ley de conservación no
constante es una señal de no linealidad del flujo y presencia de una escala pequeña en la
solución en forma de ley de potencia con exponente anómalo, α. La forma de ley de
potencia ha salido de la consideración de un grupo de renormalización más complejo en
la forma de la solución que incluye al exponente y puede tomar la forma r1 (t ) = β t .
α
En este ejemplo se ha considerado un comportamiento sincrónico entre el bombeo y la
presión dentro del yacimiento.
Debe buscarse el comportamiento más general y
proponer una coordenada de avance del frente para un comportamiento no sincrónico.
Igualmente proponer un nuevo argumento con el término wt, para considerar diferentes
frecuencia. Al aplicar el procedimiento usual aparece una nueva variable independiente
ξ = wt, que es dependiente del tiempo, y en el primer reemplazo en la ecuación
diferencial no es posible obtener la ecuación diferencial ordinaria que lleva a la solución
directamente.
3.7.6
Recuperación de niveles.
A partir del flujo en un pozo en el cual se realiza un bombeo a una tasa constante
durante un tiempo τ, se observa un abatimiento o descenso de la presión en la zona
cercana al pozo. Una vez se detiene el bombeo los niveles en el yacimiento inician un
proceso de recuperación. La observación y datos de recuperación de niveles permiten la
55
aplicación de una metodología específica para la obtención de los parámetros hidráulicos
del yacimiento. La solución clásica en este caso es considerar el yacimiento elástico y
aplicar la solución de Theis o de autosemejanza de primer orden. La solución se
presenta de la siguiente forma,
h − h0 =
Q
[ f (ut +τ ) − f (ut )] ,
4πT
(3.57)
lo cual equivale a decir que la perturbación continúa después de un tiempo τ, pero ahora
considerando que se inyecta un caudal igual a –Q, es decir, se superpone el efecto de un
pozo de caudal negativo para representar la recuperación de niveles.
Durante el proceso de bombeo la presión en el yacimiento disminuye y durante el
proceso de recuperación de niveles, la presión aumenta, observándose dos frentes
definidos de ascensos y descensos de presión.
En yacimientos elastoplásticos, la conductividad hidráulica depende de sí la presión
aumenta o disminuye, lo cual dará lugar a la componente no lineal del problema y por
tanto al planteamiento de un problema de coeficiente discontinuo. En este caso parece
ser mas compleja la localización del frente y la zona de transición de las dos fases.
3.8. TEORÍA DE PERCOLACIÓN Y EL FLUJO EN ROCAS COMO UN
FENOMENO CRÍTICO
La teoría de percolación constituye una herramienta poderosa para estudiar sistemas que
poseen un alto grado de heterogeneidad y que pueden ser representados como una red de
conductores aleatorios. Stauffer (1992), entre otros ejemplos de fenómenos críticos,
presenta una interpretación de un sistema de percolación como un modelo idealizado de
poros en rocas que almacenan un fluido y plantea las leyes fundamentales de la teoría.
En este trabajo se utilizan las leyes y conceptos básicos de la teoría, para estudiar el
medio fracturado como un fenómeno de transición de percolación.
56
El concepto de universalidad se utiliza para estudiar el comportamiento de los
exponentes críticos en transiciones de fase, ya sean las transiciones de fase basadas
en las leyes de la mecánica estadística, o la transición de percolación basada en las
leyes de los grandes números. Estos exponentes críticos se encuentran en las leyes
de potencia que expresan la forma de escalamiento de las propiedades promedias
de los sistemas en la teoría de percolación.
3.8.1
El flujo en rocas como una transición de percolación
Las rocas en la corteza terrestre han sido sometidas durante el tiempo geológico a
campos de esfuerzos los cuales han formado planos de discontinuidad o fracturas dentro
de su matriz granular.
Se considera el caso en el cual el flujo ocurre en forma
preferencial a través de esas fracturas y por tanto el flujo a través de la matriz se puede
considerar despreciable.
El sistema de fracturas estará compuesto por familias de
discontinuidades que se han originado a partir de más de un campo de esfuerzo, los
cuales han afectado la roca en diferentes direcciones y tiempos geológicos.
La experiencia de campo (observaciones y pruebas hidráulicas) ha demostrado
extensivamente que esos sistemas están lejos de formar bloques o planos
cuasihomogéneos sobre los cuales es aplicable el principio tradicional del continuo
homogéneo o volumen elemental de referencia (Hsieh, 1998; Tsang y Neretnieks, 1998).
Es por ello que las teorías clásicas que consideran medios porosos equivalentes son poco
aceptadas para aplicar en estos medios altamente heterogéneos. La teoría de percolación
aplicada inicialmente a los fenómenos críticos en materiales y ferromagnetos ha
permitido interpretar el medio fracturado como el producto de un fenómeno crítico: por
encima del umbral de percolación existe flujo a través del sistema de fracturas y bajo
dicho umbral no hay flujo.
El umbral de percolación existe cuando una ruta formada por fracturas cruza de un lado
al otro el dominio de la roca. Se habla entonces de que se ha formado el racimo
principal.
Esto ha sido producto del proceso de fracturamiento originado por la
acumulación de grandes esfuerzos que ayudan en la formación de la red de fracturas, la
57
cual, en el momento de extenderse y atravesar el dominio produce la liberación o
relajamiento de los esfuerzos. En ese punto el proceso de fracturamiento se detiene, la
fractura no se extiende sobre escalas mayores, y la red de fracturas en la roca se
considera que está muy cerca al umbral de percolación. En este contexto los sistemas
de fracturas se pueden interpretar en términos de percolación como aquellos formados
por grupos de fracturas, siendo el racimo principal un grupo mayor de fracturas a través
del cual se forma la ruta principal de flujo que atraviesa el dominio de un lado al otro.
Existirán además otros grupos o racimos aislados que no llevan a ninguna parte. Todos
estos grupos o racimos de fracturas se presentan en múltiples escalas de longitud.
Cuando el sistema se encuentra bajo el umbral de percolación, el campo de esfuerzo ha
sido insuficiente para formar un grupo de fracturas o ruta de flujo que cruce el dominio,
por lo tanto el medio está desconectado para el flujo.
La teoría de percolación estudia las propiedades físicas de un sistema en función de la
probabilidad de ocupación de sitios (análogo a densidad de fracturas) y del tamaño del
sistema y por tanto, esas propiedades son independientes de la geometría local. El alto
grado de heterogeneidad
del medio fracturado da origen a rasgos en su mayoría
desconocidos o que presentan un alto grado de dificultad para su medición. Los canales
de flujo presentes en el medio forman arreglos intrincados sobre los cuales finalmente la
geometría puede llegar a ser irrelevante.
En este trabajo se estudian las propiedades de un sistema que está muy cerca, pero por
encima del umbral de percolación, lo cual expresa la condición anterior de que los
sistemas de fracturas se encuentran muy cerca de dicho umbral e interesan aquellos
sistemas en los cuales el campo de esfuerzos ha sido suficiente para formar un grupo
principal a través del cual ocurre flujo. Se estudian sistemas conectados de fracturas
muy cerca de la transición de percolación para desarrollar los experimentos numéricos
que se presentan en el capítulo 4.
58
3.8.2 Transición de fase y umbrales de percolación
Los fenómenos físicos, productos de transiciones de fase, son muy comunes en la
naturaleza, y muestran comportamientos en diversas escalas de longitud justo en la
transición de fase o estado crítico. La forma de estudiar el comportamiento de las
magnitudes importantes en ese estado es expresar dicha magnitud, llamada parámetro de
orden, en forma de ley de potencia o ley de escalamiento. El exponente es el llamado
exponente crítico.
Una transición de fase ocurre cuando se presenta una singularidad en la energía libre o
en alguna de sus derivadas. Esto se manifiesta en un cambio brusco de las propiedades
de una sustancia. Por ejemplo, la transición de líquido a gas, de conductores normales a
superconductores o de un material semimagnético a ferromagneto. En las transiciones
de fase de primer orden hay una clara frontera entre las regiones en las cuales cada
estado es estable. Si se cruza la frontera hay un salto a la propiedad correspondiente y
un calor latente asociado. En las transiciones de fase de segundo orden hay un punto
critico, el parámetro de orden se anula, y es posible moverse continuamente de una fase
a otra. En el punto crítico existe divergencia del calor específico.
La expresión parámetro de orden se utiliza para denotar el valor promedio de la variable
que fluctúa, la cual es una señal del orden de la transición o pérdida de simetría en el
sistema.
El parámetro de orden en la transición líquido-gas, es la diferencia de
densidades, y en ferromagnetos el parámetro es la magnetización. Este parámetro de
orden desaparece en el punto crítico.
Las transiciones más conocidas y estudiadas son las llamadas transiciones de fase
térmicas (el parámetro de orden se expresa en función de la temperatura, T, y la
temperatura critica Tc) extensamente estudiadas por la física estadística (Chandler, 1987;
Yeomans, 1992). Se conoce también que la transición de fase magnética de materiales
se puede representar con el modelo Ising. El parámetro de orden, en este caso la
magnetización M, cerca al estado crítico escala de la siguiente forma
59
M ∝ (Tc − T ) β
(3.58)
En el cambio de fase de un fluido desde el estado gaseoso al estado líquido, el parámetro
de orden escala como,
ρ L − ρ c ∝ (Tc − T ) β
(3.59)
Donde β es el llamado exponente crítico que rige el comportamiento en la cercanía al
estado crítico. Existe evidencia sobre la universalidad de estos exponentes; ellos no
dependen de los detalles en la escala microscópica, sino de parámetros más
fundamentales, Así por ejemplo, en transiciones de fase de fluidos se observa un
comportamiento similar, y el exponente que caracteriza el parámetro de orden toma el
valor de β =1/3.
Aunque dos sistemas físicos parezcan diferentes, son idénticos en un nivel más
profundo, poseen un parámetro de orden que fluctúa en todas las escalas en el estado
crítico y comparten el mismo exponente crítico. Esta universalidad permite estudiar las
transiciones de fase con las mismas herramientas debido a que los exponentes críticos
son independientes de la estructura de la malla y sólo dependen de la dimensión del
sistema, lo que permite la utilización de modelos simples para describir el
comportamiento en el estado crítico de cualquier sistema.
K. G. Wilson (1979), explica la existencia de múltiples escalas en el estado crítico de
una transición de fase y propone la importante herramienta de grupos de
renormalización para la obtención de los exponentes críticos, la cual ha sido utilizada en
las transiciones estudiadas por la física matemática y ahora en la transición de
percolación, para estimación de exponentes críticos y estudio de leyes físicas y
estadísticas fundamentales.
La escala que domina el comportamiento crítico de un sistema es la longitud de
correlación, ξ, por lo tanto cualquier ley de escalamiento debe considerar esta escala en
la cercanía del punto crítico. Goldenfeld (1992) estudia la ley escalamiento de la
60
longitud de correlación y hace una comparación entre los resultados obtenidos con la
Teoría de Landau y los resultados experimentales. Las observaciones evidencian la
existencia de una escala adicional en la expresión de la longitud de correlación.
En este trabajo se estudia la analogía entre la transición de percolación regida por leyes
estadísticas y las transiciones térmicas regidas por las leyes de la física estadística, para
aplicarla a la transición en rocas.
3.8.2.1 Transición de Percolación
La transición de percolación es una transición de segundo orden que ocurre entre la fase
conectada y la fase no conectada de un medio determinado. El parámetro de orden es la
llamada resistencia P del sistema, con un estado crítico que se presenta en el llamado
umbral de percolación. Cualquier medio que se encuentre sobre ese umbral se dice que
percola, es decir, que existe por lo menos una ruta que comunica sus extremos. Un
medio que se encuentra por debajo del umbral, no percola, es decir el medio está
incomunicado.
En el umbral de percolación el medio puede poseer una ruta que
comunica sus extremos, y esa ruta se llamará el racimo principal o racimo infinito, que
estará conformado por grupos de sitios vecinos, más cercanos, unidos entre sí.
Las mallas de la Figura 3.6 tienen el mismo número de sitios ocupados, pero sólo la
segunda malla percola, si consideramos que el racimo principal se ha formado cuando el
extremo superior y el extremo inferior se encuentran unidos por una ruta de sitios
vecinos ocupados más cercanos, como se observa en la segunda malla. Se define sitios
vecinos mas cercanos, aquellos que se encuentran a los lados de un sitio, no en diagonal.
La base teórica de la percolación se encuentra en la ley estadística de los grandes
números. El proceso que rige esta teoría es la ocupación aleatoria de sitios con una
probabilidad p sobre un arreglo regular o malla, lo que se llama percolación por sitios.
En las mallas de la Figura 3.6, la probabilidad de ocupación de sitios es p=0.6, llamada
también concentración del sistema (en este caso densidad de elementos o de fracturas).
En el punto crítico de esta transición la probabilidad de ocupación es llamada
61
probabilidad crítica, pc, o umbral de percolación. El parámetro de orden se expresa en
función de (p-pc), lo cual representa la cercanía al punto crítico y escala con este
parámetro en forma de ley de potencia.
P ∝ ( p − pc ) β
(3.60)
En este caso P es la llamada resistencia del racimo infinito (fracción de sitios que
pertenecen al racimo infinito) y β es el exponente crítico correspondiente. En el punto
crítico el parámetro de orden se anula como ocurre con la magnetización en la transición
de fase magnética y con la diferencia de densidades en la transición de fase de fluidos.
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
♦
Figura 3.6 Esquema de una malla de percolación por sitios
La probabilidad crítica de ocupación o estado crítico se define en el límite, por lo tanto,
la teoría de percolación asume sistemas infinitos. Sin embargo, es conocido que los
problemas reales se presentan en escalas finitas y es allí donde se deben escalar las
cantidades importantes. El comportamiento de estas cantidades se estudiará en sistemas
de tamaños finitos.
En el punto crítico, las fuerzas que tienen influencia en escalas pequeñas, son llamadas
fuerzas de corto alcance y han sido generadas lejos del estado crítico. Éstas fuerzas
crecen en la medida en que el sistema se acerca al estado crítico y se empiezan a
desarrollar correlaciones de esas fuerzas en escalas mayores, hasta que, en el punto
crítico, se presentan correlaciones en todas las escalas y se dice que el sistema ha
62
desarrollado fuerzas de largo alcance. En el punto crítico se encuentran racimos en
todas las escalas de longitud.
En el estado crítico de la transición gas-líquido, se observan fluctuaciones de densidad
en todas las escalas, y el punto de opalescencia demuestra la existencia de las fuerzas de
largo alcance. Un poco por encima del umbral el sistema se encuentra conectado, se
dice que la correlación es infinita (permanecen fuerzas de largo alcance) y aparece el
llamado racimo infinito pero aún prevalecen racimos de diferentes tamaños (aunque se
ha desarrollado una correlación de largo alcance, no se eliminan las de corto alcance). La
longitud de correlación diverge en el punto crítico con el exponente ν, y escala de la
siguiente manera,
ξ ∝ ( p − pc ) −ν
(3.61)
3.8.2.2 Leyes de escalamiento en sistemas de tamaño finito.
Las leyes obtenidas en la teoría de percolación son válidas para sistemas muy grandes o
infinitos. Sin embargo, en la realidad, se puede acceder solamente a sistemas finitos y
por ello es necesario estudiar la discrepancia que presentan dichas leyes entre sí.
El
umbral de percolación se ha calculado para sistemas infinitos y toma un valor de 0.5927
(Stauffer, 1990). Éste es la el valor de la densidad crítica en mallas cuadradas y
percolación por sitios. Se utiliza entonces una ley de escalamiento, cerca del punto
crítico que permite estimar el valor del umbral de percolación para diferentes tamaños
de malla. En esta investigación se consideran sistemas cercanos pero por encima, del
umbral de percolación.
El escalamiento de la densidad de la masa del racimo principal, ρ ( ρ = M/LD, D expresa
la dimensión del sistema) se utiliza para plantear la base del escalamiento en tamaño
finito.
Esta densidad no es uniforme para escalas de observación menores que la
longitud de correlación (L<ξ).
Para observaciones bajo esa escala no existe una
distribución homogénea de la masa y las propiedades físicas del sistema deben
considerar la estructura detallada del racimo (propiedades microscópicas). En este caso
la masa sólo depende de la longitud, así M α LD. En el caso de observaciones en
63
longitudes L>ξ, se encuentra que la densidad permanece constante y se observa una
distribución homogénea de la masa, existiendo una dependencia adicional sobre L que
aparece a través de una función de la relación de escala (L/ξ). Por lo tanto la masa
escala como M α LD f(L/ξ). En el régimen homogéneo, la densidad, y por tanto, la
función f, está relacionada con la resistencia P, la cual considera el parámetro de orden
como en la ecuación 3.60. Así se obtiene una ley de escalamiento general para cualquier
cantidad X, en la transición de percolación, de la siguiente forma,
[
X ( L, P) = ( p − pc ) −ν g ( p − pc ) L1 / ν
]
(3.62)
donde g es una función de escalamiento.
Se deduce entonces que existen dos regímenes en los cuales las propiedades asumen un
comportamiento definido.
Por lo tanto, para toda propiedad existe una ley de
escalamiento para tamaños menores que su longitud de correlación, y otra ley para
tamaños mayores que su longitud de correlación. Ver Figura 3.7.
Densidad
p-pc=0.1
p-pc=0.01
p-pc=0.001
p-pc=0.0001
ξ1 ξ2ξ3ξ4
Régimen homogéneo, L>ξ
L
Figura 3.7. Esquema del comportamiento de la densidad de sitios conectados al racimo principal
en función del tamaño de la malla. Este tipo de comportamiento se extiende a otras propiedades.
64
3.8.2 3 Probabilidad de percolación y escalamiento
Sea π la probabilidad de que exista un racimo que percole. Esta probabilidad se conoce
con el nombre de probabilidad de percolación y claramente es función de la probabilidad
de ocupación p. Si el sistema es infinito esta probabilidad toma el valor de uno para
toda p>pc, y toma el valor de cero para toda p<pc, así la probabilidad de percolación
toma la forma de la distribución Delta para dominios infinitos, como se observa en la
figura 3.8.
Considerando que en sistemas infinitos el valor de π es igual a uno, entonces el
exponente crítico es cero y la ley general de la ecuación 3.62 permite obtener la ley de
escalamiento para la probabilidad de percolación así,
π = g [( p − pc ) L1 / ν ]
(3.63)
π
L=∞
L<∞
p
Figura 3.8. Esquema de variación de la probabilidad de percolación π en función de la
densidad de sitios.
Podemos definir la probabilidad de que la malla empiece a percolar como dπ/dp, y el
promedio de la probabilidad en la cual la malla percola por primera vez es
Pc ( L ) = ∫ p
65
dπ
dp,
dp
(3.64)
tal como se calcula el promedio de cualquier función. Si se toma como argumento de la
función g, z = (p-pc)L1/ν , y se reemplaza en las ecuaciones 3.62 y 3.63 se llega a la
siguiente expresión,
Pc ( L ) − Pc ( ∞ ) ∝ L1 / ν ,
(3.65)
donde la constante de proporcionalidad es ∫ z g´(z)dz.
La ecuación 3.65 expresa la diferencia entre la probabilidad crítica para una red de
tamaño finito Pc(L) (se puede obtener con simulaciones numéricas), y la probabilidad
crítica para un tamaño infinito Pc(∞), proporcional, en forma de ley de potencia, con el
tamaño de malla, L y el exponente crítico ν.
El umbral de percolación depende
entonces del tamaño del sistema. Conociendo los dos valores de las probabilidades
críticas se puede estimar también el valor del exponente ν. Esta es una forma de
aproximación al valor crítico del umbral de percolación en sistemas infinitos.
Podremos observar también el ancho de la transición de percolación a través de la
diferencia entre las dos probabilidades críticas y tendremos la desviación de la raíz
media cuadrada c, entre el valor en el umbral de percolación y los datos observados en
simulaciones numéricas, lo cual se puede expresar como
Δ ∝ L1 / ν
(3.66)
donde Δ es la diferencia entre Pc(L) y Pc(∞). En este caso se observa el ancho de
transición de percolación de acuerdo a la probabilidad de ocupación, y éste se interpreta
como la diferencia entre el umbral de percolación en tamaño finito y en tamaño infinito.
Este ancho varía para cada tamaño y tipo de malla, siendo, para una probabilidad de
percolación fija, mayor en sistemas de menor tamaño y menor en sistemas de gran
tamaño. Stauffer (1990).
66
3.8.2.4 Umbrales de percolación en mallas de tamaño finito.
Los criterios a seguir para la estimación de los umbrales de percolación se basan en el
comportamiento y propiedades de la probabilidad de percolación π, vista en el numeral
anterior.
Para tamaños finitos la función se suaviza y presenta dos comportamientos debido a la
ocurrencia de dos eventos: el evento de que la malla no percole o el evento de que la
malla percole; el umbral de percolación se encontrará entonces en el punto de inflexión
de la función. El punto de inflexión debe coincidir con el umbral de percolación y se
encuentra cuando π =1/2 debido a que el experimento es de tipo Bernoulli. El valor
π =1/2 es el valor que divide la distribución en dos partes iguales, por lo tanto nos
basamos en el criterio de la mediana de la distribución para estimar el umbral de
percolación. Una aproximación a la forma funcional entre π y p se puede obtener a
través de experimentos numéricos tipo Montecarlo. Con el criterio anterior se podrán
obtener los valores del umbral de percolación para mallas de diferentes tamaños.
La probabilidad de percolación π así estimada presenta una variabilidad que disminuye
en la medida que aumenta el número de realizaciones y el tamaño de la malla, como se
podrá observar en las simulaciones Montecarlo realizadas en el capítulo 4.
A partir de la ley de potencia de la ecuación 3.65 y de los valores obtenidos del umbral
de percolación en las simulaciones numéricas, es posible estimar el exponente crítico ν
para cualquier tipo de malla.
Los valores obtenidos en la literatura para sistemas de dimensión igual a dos son 0.74
obtenido por Bour and Davy (1997) desde sistemas aleatorios de fracturas, y el valor
teórico de 1/ν, igual a 0.75 reportado por Stauffer (1990), el cual fue calculado desde la
malla Bethe la cual tiene solución analítica para los exponentes críticos.
67
3.8.3 Caracterización de racimos de fracturas
A partir de las observaciones de campo, pruebas de trazadores y ensayos de pozos,
existe evidencia de que las fracturas se encuentran formando pequeños y grandes grupos
(racimos de fracturas), algunas veces totalmente aislados (no existe conexión con otros
grupos) y otras veces se encuentran unidos a una red principal de fracturas. Algunos
grupos que son interceptados por las perforaciones de los pozos muestran mayor
transmisividad que otros grupos, como puede observarse en el esquema realizado con
información obtenida en Mirror Lake, Hsieh (1998). Ver Figura 3.9.
Las pruebas de pozos en estos tipos de formaciones geológicas no pueden interpretarse
con modelos que asumen un medio homogéneo caracterizado por una escala de
observación.
Figura 3.9 Sección vertical de un campo de pozos ilustrando racimos de fracturas.
Sitio Mirror Lake, U.S.A. Tomado de Paul A. Hsieh (1998).
68
3.8.4
Estructura del racimo principal
Un sistema con probabilidad de ocupación por encima del umbral de percolación se
encuentra, en promedio, conectado por un grupo principal formado por fracturas que
atraviesan el sistema desde un extremo al otro de su dominio.
Si se coloca una
condición de flujo en los extremos del sistema, el flujo circula a través de dicho grupo de
fracturas.
Para observar la estructura de los racimos se presenta un sistema de
percolación que representa fracturas localizadas aleatoriamente en un dominio
bidimensional. En la Figura 3.10, en A, puede observarse un sistema por debajo del
umbral de percolación en el cual no se ha formado el racimo principal y se encuentra
desconectado para el flujo. En B se observa el sistema muy cerca, por encima, del
umbral de percolación, y en C su correspondiente racimo principal, extraído de la malla
mostrada en B.
Se observa que dentro del grupo principal existen algunos segmentos de fracturas que no
llevan a ninguna parte, por lo tanto no contribuyen al flujo y se llamaran segmentos
inactivos; la mayor parte de la masa que conforma el grupo principal cerca del umbral
de percolación corresponde a segmentos inactivos. Si se remueven estos segmentos del
grupo principal se obtiene la llamada columna vertebral del sistema (D en la Figura
3.10), la cual, cuando se trata de estudiar la conductividad hidráulica de estos sistemas
resulta ser muy importante ya que conforma la ruta principal de flujo. Al observar
detenidamente la columna vertebral se distinguen dos tipos de segmentos:
a. Aquellos segmentos individuales que unen varios grupos de segmentos, y son
considerados los cuellos de botella del sistema de flujo (transportan todo el flujo).
La remoción de uno de ellos, en este caso un segmento de fractura, desconecta el
sistema que se encuentra muy cerca de la transición de percolación.
b. Aquellos segmentos que se encuentran localizados entre segmentos individuales,
formando grupos en forma de circuitos cerrados, en los cuales se dividen las rutas de
flujo del grupo principal.
69
A
B
C
D
Figura 3.10 Ejemplos de patrones de fracturas con dirección aleatoria. Se observa el
sistema en diferentes estados. A, sistema bajo el umbral de percolación. B sistema
sobre el umbral de percolación. C, racimo infinito. D, columna vertebral del sistema.
Gráfico presentado por Bour and Davy, 1997.
Muy cerca, por encima, del umbral de percolación se observa que en escalas menores a
la longitud de correlación, ξ, se encuentran, en promedio, presentes los dos tipos de
segmentos mencionados en todas las escalas de longitud. Para escalas mayores, el
sistema se puede dividir en cajas de tamaño ξ, y dentro de cada una de ellas se observa
una geometría similar a la de los segmentos descritos anteriormente, conformando la
columna vertebral del grupo principal e igualmente similares en todas las escalas. La
masa de los segmentos que conforman el grupo principal posee igualmente leyes de
escalamiento en forma de ley de potencia como las leyes descritas anteriormente.
70
En el contexto de la transición de percolación, el umbral de percolación de sistemas de
fracturas ha sido considerado en función de la densidad crítica de elementos la cual se
define, como la densidad mínima de fracturas necesaria para formar un grupo principal y
conectar las fronteras del sistema. Esto permite considerar arreglos de fracturas e
involucrar algunos de sus parámetros geométricos, como la longitud y la orientación,
Robinson (1984), Hestir, Long (1990). El parámetro de percolación se define como la
densidad de objetos de longitud especifica por la unidad de área, en arreglos
dimensionales de fracturas. En arreglos de fracturas generados con estos conceptos se
ha comprobado el valor de algunos exponentes universales y se han verificado los
umbrales de percolación, confirmándose la estructura fractal de los sistemas y las leyes
de escalamiento cerca del umbral de percolación. Sin embargo, dada la complejidad de
estos arreglos no se ha considerado el flujo a través de estos sistemas de percolación
(sólo se ha estudiado el flujo sobre estructuras fractales conocidas, ejemplo tipo carpeta
de Sierpinski, Polek (1990), Acuña (1995)), ni se ha estudiado desde simulaciones de
flujo la ley de escalamiento para la conductividad hidráulica en sistemas de fracturas
cerca al umbral de percolación.
La conectividad de los sistemas de fracturas es un parámetro que ha sido considerado
para determinar la mayor o menor conexión hidráulica de un sistema. Sin embargo, la
definición de conectividad no ha sido aceptada por la comunidad científica y solamente
expresa una apariencia subjetiva de una red de fracturas. Se han presentado en la
literatura diversos enfoques de estudio de este parámetro y diversas formas de medición
(Hestir and Long, 1990; Sahimi, 1994; Berkowitz, 1993).
3.8.5
Escalamiento de parámetros de flujo desde la teoría de percolación.
La dependencia de la escala de los parámetros de un sistema es un aspecto importante
que estudia la teoría de percolación. A partir de la caracterización presentada en el
numeral 3.8.3 de los grupos de fracturas y de la presencia de ellas en diferentes escalas,
el fenómeno del flujo en rocas fracturadas se ha presentado como un fenómeno que tiene
71
rasgos importantes en diferentes escalas y es por ello necesario plantear leyes de
escalamiento de sus parámetros de flujo.
La conductividad y difusividad hidráulica de los medios porosos y fracturados, así
como la conducción de sistemas eléctricos, que se expresan como redes de conductores
aleatorios pueden ser estudiados con la teoría de percolación. El escalamiento de los
parámetros se puede expresar con una ley de potencia propuesta por esta teoría. La
difusividad es proporcional a la conductividad de estos sistemas (De Gennes, 1976) lo
cual es una manifestación de la ley de Einstein de la física estadística, que considera a la
difusividad proporcional a la movilidad (velocidad promedia de los electrones
es
proporcional a la corriente eléctrica que ellos producen).
La ley de escalamiento de la conductividad hidráulica debe ser análoga a la
conductividad eléctrica y por lo tanto se expresa aquí la ley de escalamiento para esta
última.
Stauffer (1990) compara la conductividad eléctrica con el parámetro P
(resistencia de la red y parámetro de orden) debido a que en la probabilidad p=1 todo el
material esta conduciendo electricidad y en ese estado ambos parámetros toman el
mayor valor.
En el umbral de percolación ellos se anulan, es decir, dejan de tener
sentido una vez el sistema se desconecta o el racimo principal desaparece.
Para la
permeabilidad el significado físico es que bajo el umbral de percolación no existe una
ruta (conformada por poros o discontinuidades en la roca) a través de la cual pueda
transportarse un fluido. En la Figura 3.11 se observa que las dos cantidades se anulan en
el umbral con diferente pendiente y por tanto diferente exponente crítico.
K
1
P
1/2
k
0
0
1/2
pc
1
p
Figura 3.11
Conductividad K y resistencia P en sistemas de percolación.
Tomado de Stauffer and Aharony, 1992.
72
La razón para que estos dos parámetros se anulen con exponente diferente es que la
mayor parte de la masa del racimo principal P corresponde a extremos o circuitos que
no contribuyen al flujo, por ello la mayor parte de la masa no contribuye a la
conductividad. Cerca al umbral de percolación (al aumentar la probabilidad) la cantidad
total de sitios conectados P crece rápidamente mientras la conductividad aumenta muy
lentamente con pendiente casi nula.
Basados en el comportamiento análogo de la conductividad eléctrica y de la ley general
de escalamiento de la teoría de la percolación, se deduce que la conductividad hidráulica
escala como
para L << ξ
K ∝ L− μ / ν ,
K ∝ ξ − μ /ν ∝ ( p − pc ) μ , para L >> ξ
(3.67)
de acuerdo a la escala de corte ξ. Estas expresiones constituyen la primera y la segunda
ley de escalamiento respectivamente.
El exponente crítico aquí es μ, el cual no ha sido expresado en función de otros
exponentes de la transición de percolación. A través de simulaciones Montecarlo, en
sistemas de conductores eléctricos (se utiliza la ley de Kirchoff y la ecuación de
continuidad), se ha obtenido el exponente μ= 1.3 en d=2, y μ=2.02 en d=3. Estos
exponentes caracterizan a la propiedad del medio que se encuentra relacionada con flujo,
ya sea un medio poroso, un medio fracturado, o un sistema de conductores eléctricos.
Por tanto, a partir de la teoría de percolación se obtiene una ley de escalamiento general
que posee propiedades universales.
73
4. EXPERIMENTOS NUMÉRICOS EN ACUÍFEROS FRACTURADOS
Con el fin de estudiar el comportamiento de los parámetros hidráulicos de acuíferos
fracturados que se expresa en forma de leyes de escalamiento (leyes de potencia), se
diseñó un experimento numérico basado en los conceptos de la teoría de percolación y
en los conceptos de flujo a través de acuíferos.
En el capítulo anterior se planteó la
física del fenómeno que permite estudiar el flujo en rocas fracturadas como un fenómeno
crítico cerca de la transición de fase de percolación.
En este trabajo se estudian las propiedades de un sistema que está muy cerca, pero por
encima, del umbral de percolación e interesan aquellos sistemas en los cuales el campo
de esfuerzos actuando sobre la roca ha sido suficiente para formar un grupo principal de
fracturas a través del cual ocurre flujo en el dominio considerado. Se estudian sistemas
conectados de fracturas muy cerca de la transición de percolación y se simula el flujo a
través de ellos.
Por el hecho de que los sistemas que se simulan son finitos, es
necesario estudiar el efecto del tamaño finito en el valor del umbral de percolación o
probabilidad crítica. Para ello se generan simulaciones Montecarlo y se estiman los
valores promedios de umbrales de percolación en mallas de tamaño finito.
A partir de la teoría de percolación y sus propiedades se construye un medio sintético
que representa arreglos de fracturas, el cual se puede estudiar con las propiedades
universales de la transición de percolación.
La simulación del flujo se realiza en
sistemas de fracturas en forma de arreglos radiales para representar condiciones de
bombeo. En el numeral 4.1 se describe el modelo de simulación de las mallas de
percolación y en el numeral 4.2 se describe el modelo de flujo. Igualmente se plantea la
forma de estimación de los parámetros hidráulicos y se presentan los resultados de las
leyes de escalamiento y exponentes críticos en el numeral 4.3.
74
4.1 ESTIMACIÓN DE UMBRALES DE PERCOLACIÓN EN MALLAS DE
TAMAÑO FINITO
Los criterios a seguir para la estimación de los umbrales de percolación se basan en el
comportamiento y las propiedades de la probabilidad de percolación π presentada en el
numeral 3.8.2.3. Aquí se estiman los valores de los umbrales de percolación de las
mallas estudiadas y su comportamiento o ley de escalamiento
4.1.1. Simulaciones Montecarlo
Para obtener estos umbrales se diseñó el programa PERCOLA que realiza las
simulaciones Montecarlo en mallas cilíndricas (estas mallas pueden representar una
malla radial como se describe mas adelante) y permite obtener la información de si el
sistema percola o no, para una probabilidad de ocupación dada. El sistema percola bajo
las condiciones dadas en 3.8.2.1. El llenado de los sitios ocupados se realiza a partir de
una semilla aleatoria, entre cero y uno, que se reinicia para cada realización: se asigna un
valor de cero si el valor asignado aleatoriamente es menor que el valor de p dado, y un
valor de uno, si es mayor. La información de si el sistema percola permite estimar la
forma funcional entre π, la probabilidad de percolación,
y
ocupación, después de realizar numerosas simulaciones.
Si el sistema percola se
p la probabilidad de
almacena la información de los nodos que pertenecen al racimo principal con el fin de
identificar y trazar la ruta principal de dicho racimo.
Para obtener el umbral de
percolación se estudian mallas cilíndricas de tamaños L=10, 32, 64, 128 y 256. En este
caso, el tamaño de la malla se refiere al número de sitios que tiene, ejemplo la primera
malla tiene diez por diez sitios.
Topológicamente, una malla cilíndrica se puede
construir uniendo y doblando los bordes izquierdo y derecho de una malla cuadrada.
Esta malla cilíndrica se transformará posteriormente formando las mallas radiales
necesarias en las simulaciones de flujo. Ver numeral 4.2.
El programa construido
permite que las simulaciones Montecarlo tengan en cuenta esta topología y así obtener
los valores del umbral de percolación para dichas mallas.
75
4.1.2 Probabilidad de percolación y estimación de umbrales
Para tamaños finitos la función que relaciona π y p, se suaviza y presenta dos
comportamientos debido a la ocurrencia de dos eventos, el evento de que la malla no
percole o el evento de que la malla percole. El umbral de percolación se encontrará
entonces en el punto de inflexión de la función, según criterios expuestos en el numeral
3.8.2.3.
El punto de inflexión debe coincidir con el umbral de percolación y
se
encuentra cuando π =1/2 debido a que el experimento es de tipo Bernoulli. El valor π
=1/2 es el valor que divide la distribución en dos partes iguales, por lo tanto se basa en
el criterio de la mediana de la distribución para estimar el umbral de percolación. Los
trabajos de Hestir y Long (990), Bour y Davy (1997) y Englman (1983) utilizan este
criterio para obtener los umbrales de percolación en sistemas finitos de fracturas.
La probabilidad de percolación π, presenta una variabilidad que disminuye en la medida
que aumenta el número de realizaciones y el tamaño de la malla. Con el fin de disminuir
dicha variabilidad se escogió el número de 384 realizaciones que corresponde a un nivel
de confiabilidad del 95% y un error del 5%. Este fue el criterio de estimación de tamaño
de muestra para un intervalo de confianza. El número de 384 realizaciones disminuyó
sustancialmente la desviación estándar del estimador de π para las mallas de tamaño 32,
64 y 128.
Para la malla de tamaño 256 fue necesario disminuir el número de realizaciones a 269
debido al fuerte incremento en el tiempo de computador de cada simulación. Sin
embargo, se mantuvieron desviaciones estándar del orden de las mallas anteriores. Para
la malla de tamaño 10, el número de realizaciones se elevó a 500 y se obtuvieron valores
de desviación estándar del estimador de π semejantes a las anteriores, del orden de 1.5.
Para la malla mas pequeña se utilizó un número de 500 realizaciones.
Debido a la rapidez de la simulación, en la malla de tamaño 10, se pudo observar que
para diferentes número de realizaciones el valor del umbral cambia en forma no muy
76
sensible, y en vez de hablar de un valor exacto del umbral de percolación, debe pensarse
en un estimador, como se observa en la figura 4.1.
probab.de percolación, %
54
52
384 simulac
500 simulac
1000 simulac
50
48
46
0.564
0.566
0.568
0.57
0.572
0.574
0.576
Probab.de ocupación
Figura 4.1. Umbrales de percolación en función del número de simulaciones para
tamaño de malla 10x10. El valor del umbral de percolación corresponde a una
probabilidad de percolación del 50% que en este caso corresponde a un valor de 0.569
aproximadamente.
Para cada probabilidad de ocupación se realizaron 20 simulaciones de n número de
realizaciones (500 realizaciones para la malla 10 y 384 realizaciones para las mallas 32,
64 y 128, y 269 para la malla 256). El valor del umbral se tomó como la probabilidad de
percolación que corresponde a un valor de π = 50%. Esta aproximación se obtuvo
observando el comportamiento de tres o cuatro puntos muy cercanos al rango del valor
crítico para cada tamaño de malla e interpolando para π = 50%.. Los umbrales de
percolación para las mallas de tamaño finito obtenidos se observan en la Tabla 1.
Tabla 1. Umbrales de percolación estimados para mallas cilíndricas.
Tamaño
10
32
64
128
256
P(L)prom
0.5690
0.5808
0.5857
0.5878
0.5898
P(L) significa el valor del umbral en la malla del tamaño correspondiente L.
77
Stauffer (1992) reporta un valor de umbral de percolación de 0.5927 para mallas
cuadradas de tamaño infinito y percolación por sitios. El valor del umbral cambia de
acuerdo al tipo y dimensión de la malla. En el caso de mallas cilíndricas, se esperaba
que existiera un efecto de tamaño finito ligeramente diferente al de las mallas cuadradas.
Sin embargo se puede observar en la figura 4.2 que los umbrales de percolación para
mallas de tamaño finito pc(L) tienden asintóticamente al valor del umbral en mallas
infinitas.
pc(L)
0.7
0.6
0.5
0
50
100
150
200
250
300
tamaño, L
Figura 4.2. Se observa la variación del valor del umbral de percolación para mallas de
tamaño finito, pc(L), en función del tamaño de la malla L. Los valores de los umbrales
en tamaño finito tienden asintóticamente al valor teórico para mallas cuadradas de
tamaño infinito, pc=0.5927.
Con el fin de observar las formas de la relación π vs p, se realizó el mismo proceso para
un rango de probabilidades p, entre cero y uno para los diferentes tamaños de mallas, lo
que permite comparar el comportamiento de la probabilidad de percolación y el efecto
del tamaño finito sobre dicha probabilidad. Se observa que para las mallas de tamaños
más grandes, la probabilidad de percolación se acerca más a la teórica para tamaño
infinito. Para observar este efecto se dibujaron las curvas de percolación para los cuatro
tamaños de malla analizados que se observan en la figura 4.3.
78
Probab. de percolación
100
90
80
70
60
50
40
30
20
10
0
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Probabilidad de ocupación
Figura 4.3. Efecto del tamaño finito en el umbral de percolación para mallas cilíndricas. Se
tuvieron en cuenta mallas de diferentes tamaños, L=10 se representa con rombos, L=32 se
representa con cuadrados, L=64 se representa con círculos, L=128 se representa con cruces y
L=256 se representa con asteriscos.
4.1.3
Estimación del exponente crítico ν
Para estas mallas, el exponente se obtiene a partir de la ley de potencia dada en la
ecuación 3.66 del capítulo anterior y de los valores obtenidos del umbral de percolación
en las simulaciones numéricas. Se calculan los valores de delta Δ, como la diferencia
entre el umbral de percolación de una malla de tamaño finito y el umbral en tamaños
infinitos, en este caso el valor 0.5927 y se presentan en la tabla 2.
Tabla 2. Valores de Δ para mallas de tamaño finito simuladas
Δ (p(L)-pc(∞)) Tamaño
0.0234
10
0.0114
32
0.007
64
0.0049
128
0.0029
256
79
El exponente que representa este comportamiento es 1/ν que se obtiene para este caso
desde la pendiente de la línea en el espacio log-log. El valor obtenido es de 0.65, como
se observa en la figura 4.4.
0
0
1
3
4
2
Exponente 0.65 R =
-1
log (p-pc(L))
2
-2
-3
-4
log L
Figura 4.4. Estimación en el espacio log-log del exponente ν que corresponde al valor de 1.54,
el valor presentado por la teoría de percolación es de 1.33.
Otros valores del exponente 1/ν obtenidos en la literatura para sistemas en dos
dimensiones son, el valor de 0.74 obtenido por Bour and Davy (1998) desde sistemas
aleatorios de fracturas y el valor teórico de 1/ν igual a 0.75 reportado por Stauffer
(1990), el cual fue calculado desde la malla Bethe, malla que tiene solución analítica
para los exponentes críticos.
El valor obtenido del exponente ν, confirma la universalidad de los exponentes críticos,
en este caso para la transición de percolación en mallas cilíndricas, y permite la
estimación de la longitud de correlación.
4.1.4
Longitud de correlación para sistemas simulados
Para estudiar la longitud de correlación en estos sistemas se plantea la ley general de
escalamiento de acuerdo a los conceptos presentados en el numeral 3.8.2.1. La longitud
de correlación escala como,
ξ ∝ ( p − pc) −ν
80
(4.1)
donde el exponente ν toma el valor de 4/3 en sistemas de dos dimensiones. Este valor
es la escala de corte para plantear las expresiones de escalamiento de los parámetros de
transporte en teoría de percolación.
En el umbral de percolación la longitud de
correlación diverge para las mallas analizadas.
El comportamiento de la longitud de
correlación se observa en la Figura 4.5.
Longitud de correlación
200000
150000
100000
50000
0
0.57
0.575
0.58
0.585
0.59
0.595
0.6
Probabilidad
Figura 4.5 Comportamiento de la longitud de correlación cerca al umbral de percolación para
los tamaños utilizados. El tamaño de la malla 32 se distingue con un círculo, 64 se distingue con
un cuadrado, 128 con un rombo y 256 con asterisco.
Para cada tamaño de malla L, se define la cercanía al umbral de percolación del sistema
como (p-pc), siendo p el valor de la probabilidad en la cual se encuentra el sistema, y pc,
el valor teórico (0.5927). Esta diferencia representa el valor de acercamiento al umbral
de percolación. Se buscaron valores de (p-pc) que cubren hasta cuatro órdenes de
magnitud con el fin de observar el comportamiento de una ley de potencia.
En el numeral 3.8.5 se expuso la forma de la ley de escalamiento de la conductividad de
los sistemas propuesta por la teoría de percolación de acuerdo a la ecuación 3.67. La
forma de la ley de escalamiento a utilizar se define con el criterio de la escala de corte o
longitud de correlación, ξ. Para cada forma se definen diferentes valores de cercanías al
umbral de percolación como se muestra en la Tabla 3. Se compara el valor de L, el
tamaño de la malla, con el valor calculado de la longitud de correlación en unidades de
malla.
81
La primera forma de la ley de escalamiento de la conductividad se puede obtener con
suficiente exactitud debido a que se trata de simular mallas de tamaño menor que la
longitud de correlación. En este caso los parámetros hidráulicos depende del tamaño
del sistema.
La segunda forma de la ley de escalamiento es un poco más difícil de alcanzar con
precisión debido a que se deben cubrir tres órdenes de magnitud y esto puede requerir de
mallas tan grandes que difícilmente pueden ser simuladas en computadores. La forma
de esta ley es válida para mallas de tamaño mayor que la longitud de correlación. La
mayor dificultad se encuentra en la simulación numérica del flujo, debido al extenso
tiempo de computador necesario para correr un número grande de realizaciones.
Tabla 3. Leyes de escalamiento de los parámetros hidráulicos.
p-pc (delta)
0
0.0001
0.001
0.010
0.016
0.06
0.1
Longitud
correlación
Infinito
215443
10000
464.2
248
42.6
21.5
de L < longitud de
Correlación
Todas
Todas
Todas
Todas
L=256>long. Correlación
L=128, 64>long. Correlación
Todas, L >long. Correlación
Forma
escalamiento
1ª Ley
1ª Ley
1ª Ley
1ª Ley
2ª.Ley
2ª Ley
2ª Ley
de
Nota. L es el tamaño de la malla. La longitud de correlación se calcula con la ecuación 3.61.
Las 1ª y 2ª leyes se presentan en la ecuación 3.67.
4.2 MODELO DE FRACTURAS COMO SISTEMAS DE PERCOLACIÓN
El modelo de percolación se aplica a mallas regulares y en este caso se utilizó una malla
radial para representar las condiciones de flujo radial. Sobre esta malla se construye el
arreglo de fracturas que se describe a continuación. Para percolación por sitios como es
el caso desarrollado en este trabajo, cada sitio (punto de intersección de la malla) se
encuentra activo con una probabilidad p. En la medida en que p aumenta, aumenta el
número de sitios activos y por tanto la conexión del sistema. Para cada valor de p se
tiene entonces un sistema con diferentes propiedades de conexión y por lo tanto con
82
diferentes propiedades de flujo. La unión de dos sitios vecinos ocupados forma un
segmento que representa una fractura activa.
La unión de dos sitios vecinos no
ocupados forma un segmento que representa una fractura no activa, es decir se encuentra
cerrada para el flujo. Con estos sistemas de percolación se pueden obtener grupos de
fracturas activas y no activas, de todos los tamaños, incluídos grupos aislados de
fracturas que no contribuyen al flujo total del sistema.
Para este caso de arreglos de fracturas, no se fijan ni orientaciones, ni longitudes y se
maneja un valor de abertura constante. Sin embargo para estos arreglos la longitud de
los segmentos varía de acuerdo a la longitud de dominio que se modele. Esta longitud
representa la distancia de corte entre dos fracturas.
Para crear las mallas de percolación bajo diferente valor de probabilidad p, se realizaron
simulaciones Montecarlo con el programa desarrollado para este trabajo llamado
PERCOLA, el cual permite ajustar un valor de probabilidad y observar sí la malla
percola o no. La malla percola si bajo esa probabilidad se conforma un racimo principal
de sitios ocupados, de lo contrario no percola. En la Figura 4.6 se muestran cuatro
sistemas de arreglos de fracturas en forma radial de tamaño 64 con diferentes
probabilidades o densidades. La primera con una p=0.4, bajo el umbral de percolación,
las dos siguientes con p= 0.6 y 0.7, se encuentran sobre el umbral de percolación, lo cual
significa que son sistemas conectados y se dice que percolan; sobre estos sistemas
conectados se puede simular condiciones de flujo. La malla de p=1 representa un
sistema totalmente conectado.
El programa PERCOLA permite rastrear a través del sistema la ruta formada por el
racimo principal. En la figura 4.7.1 se observa un sistema de fracturas de tamaño 128 y
densidad 0.6, el cual percola. En la Figura 4.7.2 se observa el mismo arreglo con un
trazo del racimo principal. La información del racimo principal es almacenada con el fin
de localizar sobre él los puntos de observación. Un esquema de estos puntos o pozos de
observación se observa en la figura 4.7.3.
83
Sistemas de este tipo se crearon para realizar las simulaciones numéricas de flujo en
diferentes tamaños y estados de conexión como se describe en el siguiente numeral.
1
4
3
Figura 4.6 Mallas radiales de percolación en densidad de ocupación p=0.4, p=0.6, p=0.7
y p=1. El tamaño de la malla es 64 x 64.
84
Figura 4.7.1. Sistema de fracturas sobre el umbral de percolación, con densidad de
ocupación de 0.6 y representado por una malla de percolación de tamaño 128x128.
85
Figura 4.7.2. Racimo principal en malla 128 con densidad 0.6. La línea continua es una
traza del racimo principal.
86
Figura 4.7.3. Localización de los puntos de observación que representan pozos de
observación ubicados sobre el racimo principal del sistema.
87
4.2.1 Modelo de flujo
Sobre las mallas de percolación descritas en el numeral anterior se generan condiciones
de flujo a partir de pruebas numéricas que simulan condiciones de bombeo. Estas
pruebas se realizan colocando condiciones de borde que representan el efecto de
bombeo. El caudal de explotación es una condición tipo Neuman que se impone en el
centro de la malla. El acuífero se extiende hasta el borde exterior de la malla donde se
impone la condición Dirichlet de cabeza constante, la cual representa abatimiento cero
en el límite. Esta última condición intenta expresar el efecto de acuíferos infinitos,
necesario para la interpretación de las pruebas de bombeo.
Se está interesado en obtener las respuestas transitorias de la cabeza hidráulica (nivel
piezométrico en el acuífero) a partir de las pruebas de bombeo generadas sobre el medio.
Debido a que los acuíferos simulados adquieren las características de los sistemas de
percolación, el flujo ocurre a través de aquellos segmentos que pertenecen al racimo
principal, el cual conforma la o las rutas principales de flujo en el acuífero desde el pozo
de bombeo localizado en el centro de la malla, hacia el borde exterior del acuífero. Ver
figuras 4.6 y 4.7.2. Las rutas de flujo generadas adquieren entonces las características
del racimo principal que poseen los sistemas de percolación y se estudian como tal.
Para observar el comportamiento transitorio de la cabeza hidráulica y evaluar los
parámetros de flujo en distintos sitios del acuífero, se localizaron pozos de observación
donde se hace la lectura de cabeza hidráulica.
Estos siete pozos se situaron sobre el
racimo principal del sistema mapeado por el programa PERCOLA, en forma
equidistante desde el pozo de bombeo, como se observa en la figura 4.7.3.
4.2.1.1 Programa de Flujo
Para modelar el flujo en los arreglos de fracturas anteriores y simular las pruebas de
bombeo se utilizó el programa de computador TRINET construido por Kenzi Karasaki
de Lawrence Berkeley Laboratory, California.
Este programa calcula el flujo y
transporte en redes de líneas interconectadas que simulan canales de flujo. Esto permite
88
representar las rutas preferenciales o canales de flujo en fracturas, resultado de la
heterogeneidad estudiada.
El programa de computador soluciona el campo de flujo (Ecuación 3.11) utilizando el
método de elementos finitos de Galerkin. La derivada del tiempo se trata con un
esquema de diferencias finitas.
La distribución de velocidad en la red de fracturas se
calcula a partir de la distribución de presión en cada nodo, solucionando la Ley de
Darcy, simultáneamente, para toda la malla. Esto genera una matriz cuyo ancho de
banda se minimiza reenumerando los nodos en cada ciclo, obteniendo un arreglo de
cabezas hidráulicas en cada nodo y en cada intervalo de tiempo cuando se simulan
procesos transitorios.
Debido a que los sistemas simulados tienen diferentes tamaños se deben calcular las
longitudes de los elementos, arcos y radios, que conforman
el arreglo radial,
manteniendo constante las áreas comprendidas entre cuatro segmentos consecutivos,
Además la malla que representa el arreglo de fracturas es de percolación, y no todos sus
elementos conducen.
Se construyó entonces el programa MALLA para calcular las
longitudes de los elementos y colocarlos en el formato requerido por el programa de
flujo, con las propiedades hidráulicas fijas de cada elemento.
La longitud de los
segmentos también varían con el tamaño de la malla y con el dominio del flujo. La
transmisividad y el coeficiente de almacenamiento de cada segmento se mantienen
constantes en los segmentos que representan fracturas activas y se asignan valores cero a
aquellas fracturas no activas.
Toda la información hidráulica y geométrica de los
elementos se presenta en el archivo elmt. La información de condiciones de borde del
sistema se entra en el archivo node, y la información de localización de pozos de
observación se entra en el archivo npn. Las condiciones de flujo transitorio y tiempos
se entran en el archivo ctrl.
4.2.1.2 Realizaciones de arreglos de fracturas
Con el fin de simular los arreglos de fracturas en diferentes estados de conexión y
evaluar las propiedades hidráulicas promedias de los acuíferos para cada uno de esos
89
estados, se simularon 20 realizaciones para cada estado en diversos tamaños de
sistemas. Se crearon arreglos para cuatro tamaños de malla (L =32, 64, 128 y 256) y
para tres dominios de flujo (50, 1000 y 10000 m).
Como se mencionó en el numeral 4.1, los estados de conexión de los arreglos se
expresan con diversos valores de densidad, definiendo un nivel de cercanía al umbral de
percolación. El nivel de cercanía (p-pc) se expresa como la diferencia entre la densidad
de ocupación del sistema y la densidad crítica, pc, de mallas de percolación en dos
dimensiones.
Se realizaron cerca de 5000 simulaciones de flujo para todas los tamaños
y estados de los sistemas y se leyeron respuestas transitorias de cabeza (en adelante se
denominan, transientes) hidráulica en siete pozos de observación en cada simulación.
4.2.2
Estimación de parámetros hidráulicos
Las propiedades estimadas desde sistemas de percolación equivalen a valores promedios
que permiten obtener la propiedad macroscópica de un continuo. En este caso se
simularon 20 realizaciones por cada parámetro para obtener el parámetro hidráulico
promedio, lo cual equivale a interpretar cerca de 20000 transientes. En este caso los
parámetros hidráulicos son la transmisividad y coeficiente de almacenamiento de
acuíferos fracturados. Para expresar las leyes de escalamiento se evaluaron los dos
parámetros definidos como la conductividad hidráulica (m/s) y la difusividad hidráulica
del acuífero (m2/s) considerada como la relación entre la transmisividad y el coeficiente
de almacenamiento.
Aplicando la metodología de Theis se obtuvo un valor de T y S para cada uno de los
cerca de 20000 transientes generados en las simulaciones de flujo. Dado el número de
información generada en este tipo de pruebas hidráulicas (en campo se obtiene escasa
información, en el espacio y en el tiempo) se desarrolló un programa de computador que
permitiera la evaluación en forma sistemática de los parámetros hidráulicos.
Por lo
tanto, no se utilizó la metodología gráfica de ajuste de Theis, sino que se desarrolló una
metodología que permitiera la evaluación sistemática de las pruebas de bombeo. Esta
90
metodología se basa en la optimización de dos funciones que expresan los abatimientos
a partir de una prueba de bombeo, siendo la primera una función teórica, dada por la
ecuación de flujo, y la segunda,
la función dada por el abatimiento observado en la
simulación. Los parámetros T y S se varían sistemáticamente hasta hallar la pareja de
valores que hace esta función lo más cerca posible a cero, bajo esta condición se
encuentran los parámetros que representan mejor el comportamiento del acuífero.
La
forma de la función a minimizar es,
⎛ C1 × Q ∞ e − u ⎞
delta = ⎜⎜
× ∫ du ⎟⎟ − (sobs ) ,
u u
⎝ T
⎠
(4.2)
siendo delta la diferencia entre el abatimiento teórico dado por los valores T y S y el
abatimiento obtenido dado por la simulación de flujo; C1 y C2 son constantes que
dependen de las unidades utilizadas; Q es el caudal o rata de bombeo; T es la
Transmisividad; S el coeficiente de almacenamiento, t el tiempo desde el inicio del
bombeo; sobs es el abatimiento obtenido de la simulación de flujo y u es
r2 × S
4 × C2 × T × t
que es un argumento adimensional.
El programa desarrollado se llama THEISMOD y permite la evaluación simultánea de
20 simulaciones o más, de flujo. Los parámetros de control son los errores observados,
permitiéndose hasta un error del 20%. Sin embargo, los parámetros obtenidos con las
simulaciones de flujo presentaron valores bastante menores que éste. En las figuras 4.8
y 4.9, se observan las curvas de abatimiento generadas con los parámetros estimados de
transmisividad y coeficiente de almacenamiento de cuatro transientes obtenidos en
diferentes pozos. El ajuste se observa aceptable en los tiempos medios y finales de las
pruebas.
91
Figura 4.8 Ajuste de curvas de abatimiento de la simulación numérica y abatimientos teóricos a
partir de los parámetros de flujo obtenidos con el programa THEISMOD.
92
Figura 4.9 Ajuste de curvas de abatimiento de la simulación numérica y abatimientos
teóricos a partir de los parámetros de flujo obtenidos con el programa THEISMOD.
93
4.3 DISCUSIÓN Y ANÁLISIS DE RESULTADOS
Los experimentos numéricos se realizaron para simular pruebas de bombeo en acuíferos
fracturados. Estos acuíferos son representados con sistemas de percolación que se
encuentran sobre el umbral de la transición de fase.
A partir de las condiciones
hidráulicas impuestas sobre estos sistemas, es posible estudiar comportamientos espaciotemporales que en pruebas de campo es imposible observar. Los parámetros de flujo
estudiados son la conductividad hidráulica y la difusividad (Transmisividad sobre el
coeficiente de almacenamiento) los cuales usualmente se utilizan para el reconocimiento
y evaluación hidráulica de acuíferos.
Los acuíferos estudiados se encuentran en estados diferentes expresados por la densidad
de ocupación o de segmentos abiertos para el flujo. Como se dijo antes, esta densidad
expresa el grado de conexión de los sistemas de fracturas.
Se trata, entonces, de
acuíferos que poseen un medio heterogéneo a través del cual ocurre el flujo. Un sistema
en el umbral corresponde a un valor de p igual a pc (0.5927), la cercanía al umbral se
expresa como p-pc.
El estudio y análisis de los resultados de los parámetros hidráulicos obtenidos desde
numerosas simulaciones numéricas muestran lo siguiente:
1. La localización de puntos que representan pozos de observación en forma equidistante
dentro del acuífero y la evaluación de los parámetros hidráulicos en cada uno de esos
puntos, permite observar la variación espacial de los parámetros hidráulicos entre el sitio
de localización del pozo de bombeo y el limite exterior del acuífero. Este es el llamado
efecto de escala de los parámetros de flujo: el valor de los parámetros hidráulicos
aumenta con la escala de observación. Este comportamiento se observó para todos los
estados de acuíferos simulados y para diversos dominios utilizados, ver figura 4.10.
Los resultados resultan congruentes con los expuestos por Brace (1984) y Clauser
(1990), quienes plantean el efecto de escala desde la medida de la conductividad
94
hidráulica en escalas de laboratorio y de campo. Igualmente resultan congruentes con el
comportamiento observado por Hestir and Long (1990) en simulaciones numéricas a
partir de sistemas de percolación que representan fracturas.
Se debe aclarar que la distancia mas alejada del pozo de bombeo corresponde a un pozo
de observación que debe poseer un fuerte efecto de las condiciones de frontera de los
Conductividad Hidraulica
Conductividad Hidraulica
acuíferos, como es evidente en la figura 4.10.
4.00E-04
3.00E-04
2.00E-04
1.00E-04
0.00E+00
0
2000
4000
6000
8000
10000
3.00E-03
2.00E-03
1.00E-03
0.00E+00
0
Distancia, m
PROB. B
200
400
600
800
Distancia, m
PROB. E
PROB. B
PROB. E
Figura 4.10. Se observa la variación de la difusividad y la conductividad hidráulica con la
distancia desde el pozo de bombeo o efecto de escala. La probabilidad B y E representan el
sistema simulado mas cerca y más lejos del umbral de percolación, respectivamente. Entre estos
dos valores se encuentran los demás estados simulados, los cuales presentan el mismo efecto.
Las simulaciones se realizaron sobre diferentes dominios.
2. Se observa la existencia de dos tipos de comportamiento desde el análisis de los
sistemas más cerca al umbral de percolación y aquellos más alejados.
Estas dos
situaciones dan lugar a dos diferentes leyes de escalamiento de los parámetros
hidráulicos, que se exponen en los siguientes numerales.
Los acuíferos que se
encuentran en estados más cerca al umbral de percolación (p-pc entre 0 y 0.001)
presentan los valores más bajos de conductividad hidráulica y se agrupan entre sí.
Aquellos acuíferos que se encuentran mas alejados del umbral (valores entre p-pc de
0.001 a 0.1), presentan los valores mas altos y se alejan del grupo anterior. Este
95
1000
comportamiento se observa en la figura 4.11 y se encontró en forma generalizada en
todas las simulaciones realizadas.
4.00E-04
Conductividad Hidraulica
Conductividad Hidráulica
8.00E-02
6.00E-02
4.00E-02
2.00E-02
0.00E+00
3.00E-04
2.00E-04
c
1.00E-04
0.00E+00
0
10
20
30
40
50
0
2000
Distancia, m
4000
6000
8000
Distancia,m
Figura 4.11. Se observa que los sistemas que se encuentran mas cerca al umbral de percolación
(valores de p-pc entre 0 y 0.001) presentan los valores de conductividad hidráulica más bajos y
ellos se agrupan en el nivel inferior del gráfico. Los sistemas mas alejados del umbral (valores
de p-pc entre 0.01 y 0.1) presentan los más altos valores y se distribuyen en la parte superior del
gráfico. Los que corresponden a la mayor probabilidad se identifican con asteriscos. Cruces y
círculos llenos, corresponden a sistemas de menor probabilidad respectivamente.
3. Aunque los sistemas simulados con la teoría de percolación en este trabajo,
representan medios heterogéneos conformados por sistemas de fracturas que poseen un
grado mayor o menor de conexión, existe una escala sobre la cual estos sistemas lucen
homogéneos, y es la escala de correlación o escala representativa del sistema. Se
analizaron los sistemas que se encuentran sobre dicha escala y se obtuvo la ley de
escalamiento de sus correspondientes parámetros hidráulicos.
Esta ley que se ha
llamado en este trabajo la segunda ley de escalamiento, relaciona los parámetros
hidráulicos con el valor de cercanía al umbral de percolación (p-pc), lo cual significa un
mayor o menor grado de conexión del sistema de fracturas. En este trabajo se analizan
96
10000
los sistemas de mas alto valor de conductividad, y la ley obtenida corresponde a los
sistemas un poco mas alejados del umbral de percolación (p-pc entre 0.01 y 0.1), pero
lejos de un sistema totalmente conectado (p igual a uno). El exponente μ de la ley de
escalamiento propuesto por la teoría de percolación se obtuvo para diferentes sistemas
de flujo simulados y para cada punto de observación, como se observa en la figura 4.12.
Para valores de p-pc igual o mayor a 0.2 los resultados no se ajustan a esta ley, lo cual
puede significar que a partir de este valor el sistema de fracturas deja de comportarse
como un sistema en transición de percolación.
a
-2
10
-2
-1.5
-1
-0.5
0
-2.5
Log K
Log D
9.5
-3
9
-3.5
8.5
b
-2
-1.5
-1
-0.5
Log(p-pc)
0
-2
10
-2
-1.5
-1
-0.5
9.5
0
Log K
Log D
-2.5
9
-3
8.5
-2
-1.5
-1
-0.5
-3.5
0
Log(p-pc)
Log (p-pc)
Figura 4.12. Comportamiento de los parámetros hidráulicos, difusividad y conductividad
hidráulica, en función del valor de cercanía al umbral de percolación o grado de conexión del
sistema. Este comportamiento define la ley de potencia para la ley de escalamiento de sistemas
que se encuentran en escalas mayores a la escala representativa. En a) se presentan los
resultados para una malla de tamaño 128, y en b) los resultados para una malla de tamaño 256.
Ambos con dominios de flujo de 1000 metros.
97
Se encontró el exponente μ que define la ley de escalamiento de los parámetros
hidráulicos en este régimen homogéneo cuyos valores y coeficiente de correlación se
presentan en la Tabla 4 para la malla 128 y en la Tabla 5 para la malla 256.
Tabla 4. Resultados de simulaciones en una malla 128.
Pozo No.
2
3
4
5
6
DIFUSIVIDAD
Exponente μ
1.23
1.2
1.14
1.18
1.17
CONDUCTIVIDAD HIDRAULICA
R2
Exponente μ
1.35
0.99
1.28
0.99
1.25
1
1.14
1
1.03
0.99
R2
0.99
0.99
0.99
0.99
0.99
Tabla 5. Resultados de simulaciones en una malla 256
Pozo No.
2
3
4
5
6
DIFUSIVIDAD
Exponente μ
0.63
0.71
0.7
0.74
0.77
CONDUCTIVIDAD HIDRAULICA
R2
Exponente μ
1.04
0.99
1.02
0.99
1.00
0.99
0.97
0.99
0.94
0.99
R2
0.99
0.99
0.99
0.99
0.99
La teoría de percolación presenta el exponente universal
μ con un valor de 1.30
(Stauffer, 1990) obtenido mediante simulaciones Montecarlo en sistemas aleatorios de
conductores eléctricos. A partir de simulaciones numéricas en sistemas aleatorios de
fracturas Hestir and Long (1990) presentan valores entre 1 y 2 para el exponente.
En este trabajo se obtiene los exponentes desde simulaciones Montecarlo y simulaciones
de flujo en arreglos radiales de fracturas que simulan acuíferos.
Los valores del
exponente más representativos corresponden a los obtenidos en la malla 128 donde se
hicieron simulaciones con 30 realizaciones. Los valores del exponente en la malla 256
toman valores más bajo y se obtuvieron con 20 realizaciones. Así el exponente μ
98
hallado desde las simulaciones de flujo en sistemas de fracturas de percolación se
encuentra entre 1 y 1.3, valor más cercano al propuesto por Stauffer.
Con este valor obtenido se confirma el comportamiento en forma de ley de potencia de
los parámetros hidráulicos en sistemas que poseen cierto grado de conexión, caso en el
cual es posible encontrar una escala representativa. En este caso el sistema se observa
bajo una resolución constante y no existe dependencia de la resolución del ensayo.
4. Los sistemas o acuíferos fracturados que se observan bajo la escala representativa no
lucen homogéneos, sino que muestran una estructura fractal (los racimos de fracturas
son autosemejantes) y por tanto las leyes que rigen los parámetros hidráulicos obtenidos
bajo esa escala muestran características especiales regidas por otros exponentes. Lo que
imprime estas condiciones a un sistema es la mayor cercanía al umbral de percolación,
es decir, en ese estado (cerca al estado crítico) el flujo ocurre a través de muy pocas rutas
que han sido conformadas por la formación del llamado racimo principal del sistema, al
cual llamaremos el régimen fractal. Las simulaciones realizadas para probabilidades
más cercanas al umbral de percolación presentan menores valores de conductividad
hidráulica, tal como se observa en la figura 4.11.
En este caso los parámetros
hidráulicos dependen del tamaño del sistema en estudio, y se da lugar a la llamada en
este trabajo primera ley de escalamiento.
Los comportamientos de los parámetros hidráulicos obtenidos desde las simulaciones de
acuíferos muy cerca al umbral de percolación, bajo una probabilidad de ocupación fija y
discretización de malla fija, en función de diferentes tamaños de sistemas, se observan
en la figura 4.13.
99
a
10
0
0
1
2
3
Resolucion
3.9
Resolución 3.9
-1
Log D
Log K
9.5
-2
9
-3
8.5
-4
0
Log L
1
2
3
2
3
Log L
10
0
0
1
2
3
Log D
Log K
-1
-2
-3
9
Resolución 1
Resolución 1
8
-4
0
Log L
1
Log L
b
Figura 4.13 Para dos valores de resolución o discretización de malla fija, a) Valor de resolución
3.9, b) valor de resolución 1, se observa el comportamiento de los parámetros hidráulicos, K,
conductividad hidráulica y D, difusividad hidráulica. La resolución se define como la relación
entre el dominio de flujo y el número de discretización de la malla. Se observa el
comportamiento aproximadamente constante para los diferentes pozos de observación (igual
pendiente en el espacio Log-Log).
Los valores de los exponentes obtenidos y sus coeficientes de correlación se presentan
en la Tabla 6 y en la Tabla 7 para dos valores de discretización diferentes.
Tabla 6. Exponentes de la ley de escalamiento para una resolución 3.9.
Pozo No.
2
3
4
5
6
DIFUSIVIDAD
Exponente μ/ν
-0.84
-0.85
-0.91
-0.93
-0.98
CONDUCTIVIDAD HIDRAULICA
R2
Exponente μ/ν
-1.05
0.99
-1.03
0.99
-1.10
0.98
-1.16
0.97
-1.20
0.97
R2
0.93
0.93
0.91
0.90
0.90
100
Tabla 7. Exponentes de la ley de escalamiento para una resolución de 1.
Pozo No.
2
3
4
5
6
DIFUSIVIDAD
Exponente μ/ν
-0.89
-0.88
-0.90
-0.89
-1.10
CONDUCTIVIDAD HIDRAULICA
R2
Exponente μ/ν
-0.97
0.80
-1.08
0.90
-1.25
0.97
-1.30
0.97
-1.20
0.96
R2
0.99
0.99
0.99
0.99
0.99
El exponente obtenido por Stauffer (1992) para el exponente universal en esta ley de
escalamiento es de 0.975, obtenido también desde simulaciones Montecarlo en
conductores aleatorios eléctricos. En este trabajo se obtuvieron exponentes entre 0.84 y
1.1 para un estado del sistema de p-pc igual a 0.001. Para sistemas en estados más
cercanos al umbral de percolación, como por ejemplo, p-pc igual a 0.0001 o 0, se
obtuvieron valores de exponentes entre 0.7 y 0.9.
Los valores de exponentes no
cambian sensiblemente con el sitio de medición de flujo (los sitios se encuentran
localizados sobre la ruta principal del flujo) y esto puede ser reflejo de la estructura
fractal del racimo principal. No se conocen otros trabajos que hayan obtenido este tipo
de ley de escalamiento.
5. Se observa el comportamiento de los parámetros en el régimen fractal, expresando el
valor de los parámetros hidráulicos respecto a discretizaciones del sistema. Se llama
discretización al número de particiones de la malla (tamaño de mallas), para este caso se
utilizan los valores de 32, 64, 128 y 256. La variación de ellos respecto a la resolución
de observación, dejando la escala de observación fija, es una señal del carácter fractal
del medio (Cushman, 1990). Si en el limite, resolución infinita, el parámetro tiende a un
valor de cero, el espacio vacío es considerado fractal. Este comportamiento se observa
en la figura 4.14.
101
El exponente que caracteriza este comportamiento se encontró muy cerca de uno para
todos los valores de difusividad hidráulica, observándose con mayor definición el
comportamiento en los sistemas que se encuentran muy cerca o en el umbral de
percolación.
En el caso de la conductividad hidráulica el comportamiento no se encontró definido,
presentando exponentes con valores entre 0.2 y 0.7, con valores bajos de coeficiente de
correlación. La razón puede ser que en el caso de formaciones acuíferas, la difusividad
expresa mucho mejor las condiciones físicas del sistema.
10
9
9
Log D
Log D
10
8
8
7
0
1
2
3
7
0
4
Log n
1
2
3
4
Logn
Figura 4.14. La difusividad presenta un comportamiento definido en función del valor de
discretización de las mallas utilizadas, n, que toman valores de 32, 64, 128 y 256. En el caso a)
se simuló un dominio de 1000 metros bajo un valor de cercanía al umbral de percolación de p-pc
igual a 0.001. En el caso b) se simulo un dominio de 50 metros justo en el umbral de
percolación. Los exponentes se encontraron muy cercanos a uno.
6. La longitud de correlación que en este caso expresa la escala representativa de los
sistemas utilizados, ha sido considerada como la única escala importante cerca de la
transición de fase en fenómenos críticos, en este caso la transición de percolación. Esta
longitud diverge en el umbral de percolación con un exponente ν que toma un valor de
4/3 como se plantea en el capitulo anterior.
Sin embargo, ha sido considerada
importante otra escala en el estado crítico, considerada microscópica que es relevante, y
102
que es el espaciamiento de las mallas que representan los sistemas, y que se expresa en
función de un exponente anómalo, θ, propuesto por Goldenfeld (1992).
En sistemas de percolación, esta escala representa el espacio entre dos sitios ocupados.
Para el caso de sistemas de fracturas, esta escala representa la longitud de conexión entre
dos fracturas. Bajo esta definición, el escalamiento de la longitud de correlación y
demás parámetros se expresa en función del nuevo exponente.
Con un valor del
exponente μ de 1.3, escala característica de los sistemas utilizados, 3.5 o 7.8, y valor de
difusividad encontrado, se obtiene un exponente θ de 5. Se considera que este valor no
es apropiado y que este exponente debe obtenerse a partir de estimaciones de la longitud
de correlación para lo cual se requiere mayor estudio. Un valor apropiado resulta ser del
orden de 0.8.
7. Las simulaciones de flujo permiten estudiar el comportamiento de la respuesta
transitoria de presión o de cabeza hidráulica a partir de ensayos de bombeo en los
acuíferos simulados.
El nivel de conexión de los sistemas de fracturas dado por la
probabilidad de ocupación del sistema de percolación, influye directamente en el
comportamiento de estos transientes. El comportamiento transitorio de un acuífero que
se encuentra muy alejado del umbral de percolación se estabiliza mas rápidamente que
aquellos que se encuentran cerca al umbral de percolación.
En la figura 4.15 se observa el comportamiento transitorio para diferentes estados de
acuíferos, desde aquellos completamente homogéneos donde todas las fracturas se
encuentran interconectadas, p igual a 1, y que en este trabajo hemos llamado sistemas
euclidianos (Hestir y Long,1990, comparan estos sistemas con el modelo de fracturas
infinitas propuesto por Snow, 1966), hasta aquellos acuíferos que tienen poca conexión
y se encuentran muy cerca al umbral de percolación con un valor de cercanía de p-pc
igual a 0.016.
En aquellos acuíferos que se encuentran más cerca del umbral de percolación existen
menos rutas de flujo y ellas pueden llegar a ser más tortuosas, el sistema apenas se
103
encuentra conectado y la estructura del racimo principal tiene mayor influencia, por lo
tanto la respuesta es más lenta. En los acuíferos alejados del umbral de percolación
existen mas rutas que cooperan al flujo, y por lo tanto la respuesta es más inmediata.
Los valores de la cabeza hidráulica, dependen de la física del flujo dentro de cada uno de
los sistemas.
800
cabeza hidraulica, m
a)
600
400
200
0
0
10
20
30
tiempo, min
p-pc=0.016
p-pc=0.05
p-pc=0.1
b)
cabeza hidraulica, m
400
200
0
0
10
tiempo, min
p-pc=0.016
p-pc=0.05
p-pc=0.1
p=1
Figura 4.15. Comportamiento transitorio de la cabeza hidráulica en acuíferos simulados bajo
diferentes estados de conexión, dados por los valores de p-pc. El acuífero más conectado
(p-pc=0.1) tiene una respuesta de estabilización más rápida, y el acuífero menos conectado (se
encuentra mas cerca del umbral de percolación, p-pc=0.016) demora tiempo en estabilizarse.
Esta respuesta aparece en un pozo de observación localizado a 600 metros del pozo de bombeo.
En a) se observa la respuesta para los sistemas de percolación hasta un tiempo de 30 minutos, y
en la parte b) se incluye la respuesta del sistema totalmente conectado.
104
5. CONCLUSIONES Y RECOMENDACIONES
La heterogeneidad es una característica de los materiales de la naturaleza y depende de
la escala y del lente con que se observe el medio. En la práctica no siempre es posible
observar la heterogeneidad, se tiene limitación del equipo de medida. En este trabajo se
estudia la heterogeneidad de los sistemas de fracturas en rocas, a partir de su
interpretación como el producto de una transición entre la fase no conectada y la fase
conectada. El cambio se presenta en el umbral de percolación, un sistema por debajo del
umbral no tiene rutas de flujo que atraviesen el dominio y un sistema por encima del
umbral posee un grupo de fracturas conectadas de tal manera que permiten el flujo a
través de ellas.
Esta aproximación hace posible
observar las propiedades de los
sistemas en función de la escala y del lente o resolución de observación.
Los parámetros de flujo (conductividad y difusividad hidráulica) de sistemas que
representan acuíferos fracturados se estudiaron a partir de experimentos numéricos
diseñados para simular flujo sobre sistemas de percolación que representan sistemas de
fracturas en diversos estados de conexión. La conexión de los sistemas estudiados
depende de que tan cerca se encuentren del umbral de percolación o punto critico, se
estudiaron sistemas de fracturas sobre el umbral de percolación. La lectura de los
niveles transitorios en puntos que representan pozos de observación, localizados sobre el
racimo principal de fracturas, permitió observar el efecto de escala alrededor del pozo
de bombeo. Los valores de los parámetros de flujo aumentan en la medida en que el
pozo de observación se aleja del pozo de bombeo, para todos los estados de conexión de
los acuíferos.
105
La escala de correlación o escala característica de los sistemas de fracturas fija la escala
de corte entre dos formas de observar los sistemas. Sobre esa escala los sistemas de
percolación lucen homogéneos y bajo esa escala lucen fractales. La forma de expresar
los parámetros en ambos casos es una ley de potencia, cuyo exponente crítico rige su
comportamiento.
Cuando se estudian sistemas en escalas mayores a la escala
característica, sus parámetros hidráulicos dependen del estado de conexión de los
sistemas de fracturas, y es observable el efecto de escala. El exponente crítico se
encontró a partir de la evaluación de los parámetros hidráulicos obtenidos en las
simulaciones de flujo y el valor que toma el exponente es aproximadamente igual al
valor propuesto por la teoría (Stauffer, 1992) y a los valores obtenidos en simulaciones
numéricas que consideran exclusivamente distribuciones geométricas de fracturas
(Hestir and Long, 1990).
Cuando se estudian los sistemas en escalas menores a la escala característica, sus
parámetros hidráulicos muestran condiciones especiales que se rigen con otros
exponentes, que ya no dependen del estado de conexión, sino del tamaño del sistema.
Lo que genera estas condiciones es la mayor cercanía del sistema al estado crítico o
umbral de percolación. En este caso el flujo ocurre a través de muy pocas rutas de flujo
que exhiben un comportamiento fractal, lo cual es posible observar estudiando el sistema
bajo resoluciones variables. Los parámetros de flujo toman valores bastante más bajos
que los que toman en el anterior régimen. El exponente crítico para este caso toma
valores entre 0.8 y 1.1 obtenido desde las simulaciones numéricas de flujo, valores muy
parecidos al exponente propuesto por la teoría de 0.975.
Los valores de los exponentes no cambian sensiblemente con el sitio de medición del
flujo, el cual se localizó sobre el racimo o ruta principal formada por fracturas
interconectadas entre sí. Esto puede interpretarse como un reflejo de la estructura fractal
del medio, y significa que en ese régimen el medio no tiene escala característica, por lo
tanto existe una variación continua en el espacio y en el tiempo de sus propiedades.
106
Una señal del carácter fractal del medio es la variación de los parámetros hidráulicos
con respecto a la resolución del sistema.
En este caso para un dominio de flujo
constante se cambió el tamaño de la malla, lo cual genera resolución variable. Se
observa que los valores de los parámetros disminuyen con un aumento del tamaño de la
malla. Cuando el espacio vacío es considerado fractal, en el límite (resolución infinita),
el parámetro tiende a un valor cero. El exponente que caracteriza este comportamiento
se encontró muy cerca de uno, para todos los valores de difusividad hidráulica. Para la
conductividad hidráulica no se encontró una relación definida de este comportamiento y
se piensa que esto se debe a que la difusividad hidráulica es el parámetro que asume
todas las condiciones físicas de los sistemas acuíferos.
Los conceptos de los fenómenos de transición de fase se incluyeron en el estudio de las
rocas fracturadas debido a que ellos contienen el sentido de función de correlación
(órdenes de corto y largo alcance entre los elementos del sistema), función que expresa,
en forma implícita, el estado de conexión del sistema y a su vez este estado determina
las condiciones físicas del flujo.
El hecho de que los sistemas de fracturas se agrupen en forma de racimos de percolación
y que esos sistemas se encuentren cerca del estado crítico de la transición, se debe al
efecto que causan los campos de esfuerzos sobre la roca en el tiempo geológico y a las
propiedades fisicoquímicas de la roca. Las anteriores condiciones se reflejan en el
estado de conexión de los sistemas de fracturas y deben ser tenidas en cuenta en la
caracterización de sistemas de fracturas.
La longitud de correlación que expresa la escala representativa de los sistemas
estudiados ha sido considerada como la única escala importante cerca de la transición de
fase en fenómenos críticos.
En función de esa escala se expresan las leyes de
escalamiento de los demás parámetros. Sin embargo, en el estado crítico, en sistemas de
tamaño finito, se ha sugerido relevante otra escala, considerada microscópica, y es el
espaciamiento de la malla. Goldenfeld (1992), expresa la longitud de correlación en
función de esta escala y de un exponente crítico θ.
107
En sistemas de percolación la escala microscópica está representada por el espacio entre
dos sitios vecinos inmediatos ocupados. Para el caso de sistemas de fracturas que
conforman acuíferos fracturados, esta escala debe estar representada por la distancia
promedia entre fracturas. Hestir and Long (1990) proponen la longitud promedia entre
sitios en un sistema aleatorio de líneas que representan fracturas, como la representación
de esta escala, sin embargo, ellos no asumen exponente anómalo para la nueva escala.
Las simulaciones de flujo realizadas en los experimentos numéricos permiten estudiar el
comportamiento de las respuestas transitorias a partir de condiciones hidráulica que
expresan bombeo. Estas respuestas se pueden observar para diferentes estados de los
sistemas de fracturas, sistemas con diferentes grados de conexión. Los sistemas más
conectados presentan un menor valor de cabezas hidráulicas y menor tiempo de
estabilización, lo cual sugiere que todos los segmentos de fracturas conducen fluido y
por ello la respuesta es más rápida. Los sistemas menos conectados, cerca al umbral de
percolación, presentan mayores valores de cabeza hidráulica y tiempos mayores de
estabilización. Este comportamiento sugiere que los sistemas que se encuentran más
cerca al umbral de percolación tienen pocas fracturas que conducen, el sistema apenas se
encuentra conectado y el flujo está limitado a las rutas formadas por el racimo principal.
Estas rutas tienen mayor tortuosidad y por tanto la respuesta es más lenta. Los valores
de cabeza hidráulica y de conductividad hidráulica dependen de la física del flujo dentro
de cada sistema, la cual se encuentra relacionada con el estado de conexión del sistema
de fracturas, el cual a su vez esta relacionado con la cercanía al umbral de percolación o
punto crítico. Es decir el comportamiento hidráulico en rocas fracturadas depende de las
condiciones de conexión generadas durante la génesis de los sistemas de fracturas.
Aunque se ha expresado la importancia y el efecto que tiene el grado de conexión sobre
el comportamiento de los parámetros hidráulicos y sus leyes de escalamiento en
acuíferos fracturados, en este trabajo no se ha hecho explícita una forma de medida de
esa conexión. Esta conexión, en este trabajo, sólo se expresa como una cercanía al
umbral de percolación en función de la probabilidad o densidad de fracturas. Se sugiere
108
estudiar con mas detalle el estado de los sistemas de fracturas y su conexión a partir del
origen y formación de dichos sistemas. Esto debe incluir las características y dinámica
de los eventos geológicos, además su efecto sobre la formación de las discontinuidades
en macizos rocosos.
En esta área del conocimiento se requiere mucho estudio y
desarrollo de investigación de tipo teórico y experimental, antes de ser llevada a la
práctica.
Esto con el fin de actuar en forma adecuada frente a la complejidad del
problema de flujo a través de rocas fracturadas.
Todo estudio de reconocimiento en acuíferos y yacimientos fracturados debe incluir el
estudio del estado de su sistema de fracturas y su cercanía a un umbral de flujo, de
manera que la localización de pozos de explotación y de observación sea adecuada y se
realice dentro de las rutas principales del flujo. Por lo tanto, debe estudiarse el carácter
de agrupamiento de fracturas conductoras y no conductoras en la formación. Este
reconocimiento es complejo y costoso, debe incluir un programa completo de ejecución
de pozos de observación, estudios de trazadores, además del reconocimiento geológico
regional y estudio de la dinámica de los campos de esfuerzos. Por otro lado, debe
profundizarse en investigaciones que incluyan observaciones en escalas donde los
sistemas de flujo lucen fractales, situación que posiblemente sea más común de lo que
usualmente se piensa.
Se deben estudiar sistemas de medida y observación desde el
punto de vista de la geología y la hidráulica, que permitan identificar en mejor forma el
carácter heterogéneo de este tipo de formaciones.
En yacimientos fracturados que almacenan petróleo los altos costos de la exploración
obligan aun más al desarrollo de investigaciones profundas de la física del fenómeno, en
función de la dinámica de campos de esfuerzos acoplada al fenómeno de flujo,
investigaciones que deben incluir además el efecto del almacenamiento del fluido en la
matriz de la roca.
Si bien en los sistemas de percolación generados y en las simulaciones de flujo del
experimento numérico se consideraron condiciones radiales, la conexión parcial de los
acuíferos fracturados hacen que el flujo no ocurra exactamente en forma radial, lo cual
109
genera cierta asimetría que no fue evaluada en este trabajo. Igualmente no se tuvieron
condiciones de perdidas cerca al pozo de bombeo, ni efectos de almacenamiento en el
pozo.
La consideración de los yacimientos y acuíferos como formaciones geológicas que no
siguen un comportamiento elástico, da lugar a que los parámetros de flujo sean una
función del tipo de perturbación que el hombre realice sobre ellas. La perturbación
puede crear ascenso o descenso de la presión en yacimientos o de la cabeza hidráulica en
acuíferos.
Formas de perturbaciones son las pruebas hidráulicas ejecutadas para
estimación de los parámetros de flujo y también los procesos continuos de bombeo o
inyección.
En la práctica, estas perturbaciones difícilmente se realizan a una tasa
constante. Más bien se producen variaciones y se originan ascensos y descensos de
presión en forma continua dentro las formaciones, lo cual da lugar a una respuesta no
elástica de la roca. El comportamiento no elástico es un problema que debe asumirse y
tratarse con las leyes físicas y soluciones matemáticas adecuadas, de manera que se
pueda llegar a plantear modelos físicos y matemáticos considerando la no linealidad de
los procesos de flujo.
Se mostró que el efecto del comportamiento no elástico de las rocas puede ser tenido en
cuenta en las ecuaciones de flujo. Ese efecto da lugar a un problema de flujo no lineal y
coeficiente discontinuo, su solución debe buscarse a partir de una ley de autosemejanza
especial. Esta ley se basa en la consideración de un grupo de renormalización más
complejo que el simple utilizado por el análisis dimensional, para lo cual resulta
fundamental el estudio de Grupos de Lie. Este grupo expresa una ley de homogeneidad
en función de una escala adicional que no es irrelevante, el radio del pozo, y un
exponente anómalo que aparece debido al fenómeno físico. El planteamiento de una ley
de conservación no integral, a la cual da origen la física del problema, permite obtener
una expresión para su exponente anómalo y por tanto su solución completa.
En este trabajo se propone el bombeo periódico como un problema de coeficiente
discontinuo que posee las características anteriores.
110
Se hace la suposición de
autosemejanza de segundo orden y la definición de la coordenada de un frente de carga y
descarga con algunas simplificaciones. Se sugieren algunos lineamientos físicos y
matemáticos para considerar la nueva suposición que lleva a obtener la solución del
problema.
Es pertinente continuar con el estudio de este problema siguiendo los
lineamientos propuestos y obtener la expresión del exponente anómalo que rige la
solución asintótica. Vale la pena resaltar una vez mas que escalas no consideradas
importantes hasta el momento hacen parte de la solución cuando se considera la nolinealidad del flujo, en este caso la escala del radio del pozo.
111
REFERENCIAS
Acuña J, Yortsos Y. 1995. Application of fractal geometry to the study of networks of
fractures and their pressure transient. Water Ressources Research. Vol 31,
No.3.pp527-540.
Barenblat, G.I.
1996.
Scaling, self-similarity, and intermediate asymptotics.
Cambridge University Press. 384 p.
Barenblat, G.I, Entov, V, Ryzhik, V. 1990. Theory of Fluid Flows Through Natural
Rocks. Kluwer Academic Publishers, London.
Barenblat, G.I. 1987. Dimensional Analysis. Gordon and Breach science publishers.
Barenblat G.I, Zheltov Iu. P, Kochina I. N. 1960. Basic Concepts in the Theory of
Seepage of Homogeneous Liquids in Fissured Rocks (STRATA). PMM Vol. 24, No.
5. Pp852-864.
Barker J.A. 1988. A Generalized Radial Flow Model For Hidraulic Test in Fractured
Rock. Water Resources Research. Vol 24, No. 10.pp1796-1804.
Beacher, G. B., Lanney, N.A. Einstein, H.H. 1977. Statistical description of rock
properties and sampling. Proceedings, 18th Symp.Rock Mechanics. p5C1, 1-8.
Bear, J. 1972. Dynamics of fluids in porous media. Dover Publication.
Berkowitz, B, Balberg, I. 1993. Percolation Theory and Its Application to Groundwater
Hydrology. Water Resources Research
Berkowitz, B. 1995. Analysis of fracture network connnectivity using percolation
theory. Mathematical Geology. V27. N.4
Billaux , D. Chiles, Long, J. 1989. Three dimensional statistical modelling of a fractured
rock mass. Int. Journal Rock mech. V26. 1989.
Billaux, D. 1990. Hydrogeologie des milieux fractures. Geometrie, connectivite et
comportment hydraulique, Ph.D. Thesis, Ecole des Mines. Paris, 277 p.
Bour O, y Davy. 1997. Connectivity of random fault networks following a power law
fault lenght distribution. Water Resources Research. Vol 33, No. 7. pp1567-1583.
112
Brace, W.F. 1980. Permeability of crystalline and argillaceous rocks. International
journal of rocks mechanics and mining sciences. Vol.17. pp. 241 - 251.
Brace, W.F. Walsh, J.B. 1984. Permeability of crystalline rocks: New in situ
measurements. J. Geophys. Res. 89 (B6), 4327.
Brown, S. 1987. Fluid flow throught rock joints: the effect of surface roughness. Journal
of Geophysical Research. Vol. 92.
Cacas, M.C, Ledoux,E, Marsily,G. 1990. Modeling fractures flow with a stochastic
discrete fracture network.1.The flow model. Water Resources Research, vol 26 (3).
Carrillo, J. 1997. Flujo en medios diaclasados. Universidad Nacional de Medellín.
Carslaw, H.S. Jaeger, C. 1959. Conduction of Heat in Solids. Oxford. Clarendon Press.
Campos, A. 1995. Iniciacion en el analisis de ecuaciones diferenciales mediante grupos
de Lie. Universidad Nacional de Bogota. Preprint, 256p.
Clauser Ch. 1992. Permeability of Crystalline Rocks. Eos, Transacttttions, American
Geophysi. Vol 73, No.21. pp232-237.
Cole J, Wagner B. 1996. On Self-similar solutions of Barenblatt´s nonlinear filtration
equation. Eeuro. Jnl of Applied Mathematics. Vol. 7.pp151-167.
Cook, N.G. 1992. Natural joints in rock: Mechanical, hydraulic and seismic behaviour
and properties under normal stress. Int. Journ. Of Rock Mech. V29. N.3. p.193-197.
Charlaix E, Guyon E, Roux S. 1987. Permeability of a Random Array of Fractures of
Widely Varying Apertures. Transporrt in Porous Media 2.pp31-43.
Cushman, J.H. 1990. Dynaymic of fluids in hierarchical porous media. An introduction
to hierarchical porous media. Academic Press. 505p.
Dagan, G. 1986. Statiscal theory of groundwater flow and transport: pore to laboratory,
laboratory to formation and formation to regional scale. Water Resources Research.
Vol.22. No.9. pp. 120S-134S (14).
Dershowitz, W.S. y Einstein H.H. 1988. Characterizing rock joint geometry with joint
system models. Rock mechanics and rock engineering. Vol. 21. pp. 21 - 51 (31).
Earlougher , R.C. 1977. Advances in Well Test Analysis. Monograph Series. Society
of Petroleum Engineers of AIME. Dallas.
Englman, R. Gur, Y. Jaeger, Z. 1983. Fluid flow through a crack network in rocks.
Journal of applied mechanics. Vol. 50. p707-711.
113
Follin, S. 1990. Effective hydraulic conductivity for fractured media. Procc. Calibration
and reliability in groundwater modelling. The hague. IAHS. No.195.
Gelhar, L. 1986. Stochastic subsurface hydrology from theory to applications. Water
Resources Research. Vol.22 No. 9 pp.135S-145S.
Gentier, S.D. Billaux, D. 1989. Laboratory testing of the voids of a fracture. Rock
Mech. Rock. Eng. 22(2), 149-157.
Goldenfeld, N. 1992. Lectures on phase transitions and the renormalization group.
Addison-Wesley publishing company.
Gomez, S, Mesa, O.J, Cogollo, C.F. Rojas, L. 2000. Difusion de fluidos en suelos no
saturados considerando conductividad hidraulica discontinua. Ingenieria Hidraulica
en Mexico. In press.
Gomez, S, Mesa, O.J. 1997. Flujo en rocas usando grupos de renormalización. Revista
Colombiana de Fisica. V30. No.1.
Gómez, S. 1994. El problema del Flujo en Medio Fracturado.
2º.Congreso
Latinoamericano de Hidrología Subterránea. ALHSUD. V2. P585-595.
Gringarten A. C. 1982. Flow-Test Evaluation of Fractured Reservoirs. Geological
Society of Ameerica. Special Paper 189. pp237-263.
Hestir K, Long C. S. 1990. Analytical Expressions for the Permeability of Random TwoDimensional Poisson Fracture Networks Based on Regulaar Lattice Percolation and
Equivalent Media Theories. Journal of Geophysical Research. Vol 95, No. B13.
Pp21,565-21-,581.
Hestir, K. y Long J. 1990. Analytical expressions for the permeability of random twodimensional poisson fracture networks based on regular lattice percolation and
equivalent media theories. Journal of Geophysical research, Vol. 95, No.B13, pages
21,256 - 21,581.
Hsieh, P. 1998. Scale effects in fluid flow through fractured geologic media. Scale
Dependence and Scale Invariance in Hydrology. Garrison Sposito. Cambridge.
Hsieh, P.A. Neuman S.P. 1985. Fiel determination of the three-dimensional hydraulic
conductivity tensor of anisotropic media. 1. Theory. Water Resources Research. Vol
21. No.11. pp. 1655-1665 (10).
Hsieh, P.A. Neuman S.P. 1985. Field determination of the three- dimensional hydraulic
conductivity tensor of anisotropic media. 2. Methodology and application to fractured
rocks. Water Resources Research. Vol. 21. No.11. pp. 1667-1676.
Kadanoff, L.P. 1966. Scaling laws for Ising model near Tc. Physics 2 (6) 263-272.
114
Karasaki K, Long J, Witheerspoon P. 1998. A new analytic model for fracturedominated reservoirs. SPE Formation Evaluation.pp242-250.
Karasaki, K. 1987. Well Test Analysis in Fractured Meda. Ph.D Thesis. Lawrence
Berkeley Laboratory. University of California.
Kirkpatrick (1983). Science. 220, 671.
Kolmogorov, A. N. 194. The local structure of turbulence in incompressible viscous
fluid for very large Reynolds numbers. Dokl. Akad. Nauk. SSSR30. 299- 303.
Lomize, G.M. 1951. Strömung in klüftigen Gesteinen (Russian). Gosenergoizdat.
Long, J.C. Remer J.S. Wilson, C.R. Witherspoon P.A. 1982. Porous media
equivalentes for networks of discontinuos fractures. Water resources research. Vol
18. No.3 pp. 645-658.
Long, J.C. Witherspoon P.A. 1985. The relationship of the degree of interconnection to
permeability in fracture networks. J. Geophys. Res. 90(B4) 3087- 3097.
Louis, C. 1967. Strömungsvorgänge in klüftigen medien und ihre wirkung auf die
Standsicherheit von bauwerken und Böschungen im fels. Publications of the institut
für Bodenmechanik und Felsmechanik at the TH Karlsruhe, Vol 3.
Marsily, G. 1986. Quantitative hydrogeology: Groudwater hydrology for engineers.
Academic Press. 440p.
Neuman, S.P. 1998. Scale Dependence and Scale Invariance in Hydrology. Garrison
Sposito. Cambridge.
Neuman, S.P. Winter C.L. Newman C.M. (1987). Stochastic theory of field-scale
fickian dispersion in anisotropic porous media. Water Resources Research. Vol.23.
No.3. p.467-479.
Niemi, A. 1994. Modeling flow in fractured medium. Uncertainty analysis with
stochastic continuum approech. Ph.D. Dissertation. Technical research centre of
Finland.
Odling, N. 1991. A conductance mesh approach to the permeability of natural and
simulated fractures patterns. Water resources research. V27 (10).
Neretnieks, J. . 1987. Channeling effects in flow and transport in fractured rocks.
Internationation Symposium Swed. Nucl. Power. Stockholm.
Perez, D. 1995. La explotación del agua subterránea. Un nuevo enfoque. Editorial
Científico-Técnica.
115
Polek, J. 1990. Studies of the hydraulic Behaviour of Hierachically Fractured Rock
Geometries. M.Sc.Thesis. Lawrence Berkeley Laboratory. Univ. Of California.
P72.
Pyrak-Nolte, L. 1987. Hydraulic and mechanical properties of natural fractures in low
permeability rock. 6. Int. Cong. On rock mech.
Raven, K. y Gale, J. 1985. Water flow in natural rock fractures as a funciton of stress
and size sample. Int. J. Rock Mech.
Raven, K. 1988. Interpretation of field tracer tests of a single fracture using a transient
solute storage model. Water Resources Research. V.14.
Robinson , P.C. 1983. Connectiivity of fracture sistems-a percolation theory approach. J.
Phys A: Math. Gen. 16. pp605-614.
Robinson , P.C. 1984. Numerical calculations of critical densities for lines and planes. J.
Phys. A: Math. Gen. 17. pp2823-2830.
Rofail N. 1967. Analysis of Pumping Test in Fractured Rocks. Hidrology of Fractured
Rocks. Vol 1. AIHS-Unesco. pp81-88.
Romm, E.S. 1966. Flow characterstics of fractured rocks. Moscow, Nedra.
Sahimi M, Yortsos Y. 1990. Aplications of Fractal Geometry to Porous Media: A
Review. SPE 20476.pp1-25.
Sahimi, M. 1994. Applications of percolation theory. University of southern California.
Taylor&Francis.
Sahimi, M. 1995. Flow and transport in porous media and fractured rock. From classical
methods to modern approaches. VCH.
Serrano S. 1997. The Theis Solution in Heterogeneous Aquifers. Ground Water. Vol.35,
No. 3.pp463-467.
Sharp, J. C. 1970. Fluid flow through fissured media. Ph.D. Thesis, University of
London.
Snow, D.T. 1969. Anisotropic permeability of fractured media. Water Resorces
Research. Vol. 5. No.6. pp. 1273 - 1289 (16).
Snow, D.T. 1965. A parallel plate model of fractured permeable media. Ph.D.
Dissertation, Berkeley, University of California. 331p.
Stauffer, D. Aharony A. 1992. Introduction to percolation theory. Second Edition.
Taylor&francis.
116
Streltsova T.D. 1988. Well testing in heterogeneous formations. John Wiley&Sons.
Theis, C.V. 1935. The relation between the lowering of the piezometric surface and the
rate of duration of discharge of a well using groundwater storage. Trans. Am.
Geophys. Union. V16. p519-524.
Tsang C.F. y Neretnieks, I. 1998. Flow channeling in heterogeneous fractured rocks.
Reviews of geophysics, 36, 2. pp. 275-298.
Tsang, Y.W. 1987. Channel model of flow through fractured media. Water Resources
Research . Vol 23. N 3. P467-479.
Wilson, K. 1979. Problems in physics with many scales of leght. Scientific American.
241, 2, 158-179.
Warren, J.F. Root, P.J. 1963. The Behaviour of Naturally Fractured Reservoir. Society
of Petroleum Engineers Journal. AIME. V228. P245-255.
Witherspoon, Wang, et al, P.A. 1980. Validity of cubic law for fluid flow in a
deformable rock fracture. Water resources research.
Wittke, W. 1990. Rock Mechanics. Springer-Verlag. 1075p.
Yeomans, J.M. 1992. Statistical Mechanics of Phase Transitions. Clarendon Press
Oxford. 153p.
117