Download Linear Poroelasticity

Numerical approximation of reservoir fault stability
with linear poro-elasticity
Aproximación de la estabilidad de falla de
yacimientos con poroelásticidad lineal
Duran Triana O. & Devloo P. R. B.
Universidade Estadual de Campinas,
Campinas, Brasil
[email protected], [email protected]
Abstract
En este documento se resumen los resultados del trabajo [1], conclusiones, recomendaciones para
futuras investigaciones y abarca rápidamente los siguientes puntos:
• Revisión teórica de las ecuaciones constitutivas para el modelo poroelástico monofásico y sus
principales hipótesis siendo consideradas.
• Forma adimensional de la formulación, conceptos teóricos y una descripción matemática del método
semi-analítico usado para la reactivación de falla utilizando teoría de inclusiones.
• La formulacion débil, discretización por elementos finitos continuos, implementación computacional
y validación con algunas soluciones analíticas.
• El efecto de diferentes estrategias de producción, así como el contraste de rigidez del yacimiento
- rocas circundantes, en el inicio de la reactivación, es estudiado.
I. Motivación
Modelos existentes para reactivación de falla
Existen un gran número de casos de gran interés en el estudio de reactivación de fallas,
subsidencia, deformación del yacimiento y cambios en el perfil de esfuerzo debido a la extracción de fluidos e.g [2], [3], [4], [5], [6], [7], [8]. En esta sección modelos existentes para estudiar
el efecto de la depleción del yacimiento son revisados.
En la mayoría de los casos, el reservorio es tratado como una inclusión poroelástica rodeada
por un material puramente elástico. De forma convencional en [9]; [10]; [11]; [12], dos abordages
son implementadas
• El yacimiento es modelado como una inclusión poroelástica con propiedades similares a
las de las rocas circundantes, pero con presiones de poro variables en el espacio [2]; [3];
[13]; [14].
• El yacimiento es modelado como una inclusión elipsoidal con propiedades diferentes a las
de las rocas circundantes, pero con pressiones de poro uniformemente distribuidas sobre
la inclusión [15] [7], [16].
La figua 1 muestra una respresentación esquemática de una inclusión poroelástica que
modela el yacimiento con una variación de presión ∆P, induciendo alteraciones en el perfil
de esfuerzos en las rocas adyacentes, que pueden reactivar parte de las fallas (Lineas rojas)
pre-existentes (Lineas verdes).
Fig. 1: Problema de reactivación de falla divido en dos partes. Izquierda) Problema inicial.
Derecha) problema de perturbación. Inclinación de la falla θ , zonas verdes respresentan partes
reactivadas en los planos de falla. Fuente: redibujado de [17].
Método de Galerkin
El uso apropiado de espacios de aproximación para problemas de consolidacion, es atribuido
a [18], en su estudio numérico en condiciones de deformación plana, se presentan espacios
de aproximación que satisfacen la condición Ladyzhenskaya-Babuska-Brezzi (LBB) [19] para la
aproximación de los desplazamientos y la aproximación de la presión de poro. Sin embargo esta
propuesta no es muy popular, debido a que en situaciones pacticas su implementación requiere
de programación compleja. En este trabajo es usado una librería neopz para el desarrollo de
algoritmos que se basan en el método de elementos finitos, en donde diferentes ordenes de
aproximación pueden ser elegidas de forma arbitraria, tanto para los desplazamientos y la
presión de poro, teniendo como ventaja que la condicion LBB es facilmente satisfecha. [20]
muestra que para pequeños pasos de tiempo, oscilaciones no físicas pueden ser introducidas.
no obstante desde que la condición LBB es satisfecha una aproximación razonable es usar
elementos cubicos - cuadráticos, cuando los desplazamientos y la presión de poro son las
variables de estado.
II. Objetivos
El objetivo principal es el desarrollo de modelo un computacional simple y robusto que permita
realizar un rápido análisis parametrico para identificar los puntos mas críticos de la reactivación
de falla causados por la producción o injección de fluidos del yacimiento. Idealmente cualquier
modelo debe ser calibrado con datos reales para aplicaciones practicas, esta tarea es dejada
para estudios futuros.
El objetivo científico es alcanzado mediante los siguientes hitos:
• Entendimiento de los mecanismos fundamentales de la reactivación de falla, inducida por
producción/injección de fluidos a través de efectos poroelásticos.
• Desarrollo de un modelo computacional, simple y robusto que permite cuantificar la depleción del yacmiento para tener encuenta ambas abordages, presión de poro no uniforme
junto con constraste de materiales entre el yacimiento y las rocas circundantes
• Analizar el inicio de la reactivación causado por la perturbación de un estado en equilibrio,
ver figura 1.
• Identificar los parámetros más importantes que afectan la reactivación de falla.
III. Metodología
A. Consideraciones
Generalmente existe interacción entre el fluido y la roca, el fluido fluye a través de la roca,
como consecuancia de los ciclos de descarga/carga, consecuentemente la presión de poro
varía con el espacio y el tiempo interferiendo con el proccesso de deformación. Para analizar
esta situación, son introducidas las siguientes consideraciones:
• No existen efectos incerciales, i.e. un proceso cuasi-estático.
• Régimen de pequeñas deformaciones, i.e. elásticidad linear.
• Flujo isotérmico.
• Flujo dominado por difusión de masa.
• Sólido isotropico completamente saturado.
• El tensor de permeabilidad es diagonal.
• Formación rocosa compresible. La velocidad del fluido es caracterizada por la lei de Darcy
(la derivación relacionada con Navier-stokes es dada en [21]). Algunas consideraciones
son intrínsecamente relacionadas con la lei de Darcy:
• El fluido es homogéneo y newtoniano.
• No ocurre reacción química, precipitación ni absorción.
• No hay efectos electro-cinéticos.
B. Formulación fuerte y débil
Usando el procedimiento en [22], el sistema completamente acoplado es descrito y transformado mediante el uso de los grupos adimensionales presentados en [1], como también su
correspondiente direcretización espacio-temporal.
C. Discretización espacio-temporal
Siguiendo [1] el sistema semi-discreto es dado por:
u Kelasticity −Qc
ū
F
=
T
p̄
Fp
Qc
S+H
(1)
Una vez la discretización espacial ha sido realizada, el sistema 1 respresenta un conjunto
de ecuaciones ordinarias en el tiempo. Por conveniencia las ecuaciones son escritas en la
siguiente forma:
u d
0 0
ū
F
Kelasticity −Qc
ū
=
(2)
+
Fp
0
H
p̄
QTc S dtD p̄
Diferencias finitas en el tiempo son usadas para tratar el problema de valor inicial resultalte.
Considerando la siguiente ecuación:
dx
+Cx = F
(3)
dt
La discretización en el tiempo es realizada por la regla trapezoidal generalizada [23]. Esta
representa un método explícito, para la aproxcimación de:
(xn+1 −xn )
dx
dt n+ξ =
∆tD
(4)
xn+ξ = (1 − ξ ) xn + ξ xn+1
H
donde ∆tD es el paso de tiempo, xn y xn+1 , son los vectores de estado en los tiempos
tDn y tD n+1 , y ξ es el parámetro que determina el nivel de implicitud, el cual es acotado por
0 ≤ ξ ≤ 1. el valor de ξ puede ser obtenido desde las tabla I, que describe diferentes esquemas
temporales.
La aplicación de este procedimiento a la ecuación 2 y multiplicando por −1 el conjunto de
ecuaciones correspondiente al problema de difusión, se recupera la simetría obteniendo el
siguiente sistema:
onKelasticity
−Qc
ū
T
−Qc
−S − ξ ∆tD H n+ξ p̄ n+1
(5)
0
0
ū
Fu
=
+
∆tD F p n+ξ
QTc S − (1 − ξ ) ∆tD H n+ξ p̄ n
TABLE I: Esquemas temporales caracterizados por diferentes valores de ξ . Fuente: [23].
ξ
ξ
ξ
ξ
ξ value
0
1
0.5
0.6667
=
=
=
=
ξ=
1+
1
∆tD
Esquema Temporal
Método de Euler progresivo
Método de Euler regresivo
Método de Crank-Nicolson
Método de Zienkiewicz Galerkin
1
− ln(1+∆t
Método logaritmico de Sanhu
D)
!
ξ=
1+
tk
tk+1 −tk
−
t1
k−1−t
ln 1+ t k
Método logaritmico de Hwang
k
[23] prueba que ξ ≥ 12 corresponde a métodos incondicionalmente estables en el tiempo.
El problema inicial es dado por la solución de:
Kelasticidad
−Qc
ū
Fu
=
(6)
p̄
∆tD F p n+ξ
QTc
S + ξ ∆tD H
El problema de perturbación en la figura 1 es solucionado por el sistema (5), detalles sobre
el procedimiento de solción son dados en [1]
Kelasticidad
−Qc
Rigidez
−QTc
−S − ξ ∆tD H n+ξ
(7)
0
0
Matriz Masa
QTc S − (1 − ξ ) ∆tD H n+ξ
0
0
ū
Fu
Vector cargan+ξ = −
+
(8)
∆tD F p n+ξ
QTc S − (1 − ξ ) ∆tD H n+ξ p̄ n
D. Implementación
En resumén, la implementación de la formulación fuerte usando la librería neopz se puede
resumir a un punto clave, el cual es la traducción de la formulación debil in lineas de codigo
neopz.
La figura 2 muestra el flujo de trabajo para generar los modelos computacionales en la
herramienta desarrollada LPA (Linear Poroelastic Analysis). La información es procesada desde
la geometría, creada por un programa enmallador (GID) e importada en LPA, con base en
la geometría y las propiedades de material, condiciones iniciales y de contorno, es creado el
modelo computacional multifísico. Finalmente los resultados son vizualidados usando Paraview.
Fig. 2: Flujo de trabajo en LPA.
IV. Validación y Resultados
De manéra resumida en este documento, se presentan algunos resultados obtenidos en la
investigación. A continuación se muestran una serie de resultados que verifican y que ilustran
como el problema siendo de reactivación de falla tratado.
A. Modelo de yacimiento Rectangular
Reservoir Material
Hosting Material
Fig. 3: Yacimiento rectangular representando la inclusión y condiciones de contorno.
Considerando una inclusión rectangular presentada en la figura 3, donde Ra , Rb , SBRa , SBRb ,
representa respectivamente: espesor del yacimiento, extensión del yacimiento, dimensión de las
rocas circundantes. RD es la profundidad medida de la superficie hasta el centro del yacimiento.
La inclusión es depletada ±∆Pex (e.g. ∆Pex < 0 representa producción) para todo x ∈ Ω (e.g. una
distribución de presión de poro homogénea. Las condiciones de contorno son: esfuerzo normal
cero en la superficie σn = 0, en las laterales los desplazamientos son cero en la dirección x,
en el basamento los desplazamientos son cero en la dirección y. El problema es encontrar los
cambios de esfuerzos, dadas las condiciones de contorno y la distribución de presión de poro.
Las figuras 4, 5 y 6 muestran λFR para fallas inclinadas en 60 grados antihorarios medidos
desde la dirección positiva del eje x, en diferentes offset (offset reprsenta la distancia del
afloramiento de la falla hasta el origen del eje x, ver figura 3).
Fig. 4: Factor de reactivación para un yacimiento rectangular, Izquierda) fallas con inclinación
de 60 grados. Derecha) λFR a lo largo de los planos de falla. Relación de arqueo (Arching ratio)
normalizados por (1−2ν()
(1−ν) .
Fig. 5: λFR a lo largo de planos de falla. Izquierda) Offset 0. Derecha) Offset 1. λFR normalizados
por (1−2ν)
(1−ν) , (λFR < 0 tendencia a reactivar).
Fig. 6: λFR a lo largo de la extensión del yacimiento. λFR normalized by
a reactivar).
(1−2ν)
(1−ν) ,
(λFR < 0 tendencia
El análisis de reactivacion es basado en la localización del contorno λFR = 0, desde que este
representa un cambio de estado de compresión a extensión. Usando este criterio cualquier
plano de falla con inclinación de 60, que es intersectado por λFR = 0 tiene tendencia a reactivar,
donde λFR es negativo y tendencia a estabilizar donde λFR es positivo. El analísis anterior aplica
cuando un régimen de fallamiento normal es considerado en la región de interés, en el caso
de régimen inverso los signos son invertidos.
B. Threshold de reactivación
Sismos inducidos con modulo de ∆CFS mayor que 0.1 MPa han sido encontrados en condiciones donde ∆CFS es cerca de zero [13], [24], [7], [12], [8].
Si este valor es asociado con la consistencia del material que compone la falla, sería posible
0.1
t
definir un valor límite (Threshold) λFRt = ∆CFS
α∆pex = α∆pex , todos los valores negativos estrictamente
representan una fuerte tendencia a reactivar.
Asi, usando α = 0.8 y ∆Pex = −10 MPa, es equivalente a λFRt = −0.0125, esto puede ser
interpretado en los resultados anteriores como:
• Figura 4 muestra que las fallas con mayores offsets, a presar de tener valores negativos,
tienen una tendencia de reactivación nula.
• Figura 5 muestra que la falla en offset 0 exibe una tendencia de reactivación a profundidades localizadas en el flanco izquierdo del yacimiento.
• Figura 6 muestra que cualquier plano con inclinación de 60 que pasa por el yacimiento
reactiva a lo largo de la extension del mismo, como el caso de la falla con offset 1.
C. Variación del coeficiente de frición
El coeficiente de fricción en la mayoría de las rocas puede variar desde 0.4 a 0.8 [8], [25].
La figura 7 muestra el efecto de la variación de µs en el contorno nullo para fallas incinadas
60 grados.
Fig. 7: Efecto de µs en el contorno λFR = 0. Factor de reactivación normalizado por
(1−2ν)
(1−ν) .
D. Contraste de materiales
El contraste de materiales es un factor importante. En relación, al caso inicial con módulos
elásticos similares (caso base), se presentan dos casos en relación a la rigidéz que compone
el yacimiento:
• Caso 1: Roca rígida, un decimo del caso base case (e.g. λcase1 = 0.1 ∗ λbase and µcase1 =
0.1 ∗ µbase ).
• Caso 2: Roca dúctil, diez veces mayor que el caso base (e.g. λcase2 = 10.0 ∗ λbase , µcase2 =
10.0 ∗ µbase ).
La figura 8 muestra que con el mismo conjunto de fallas, el contraste del material tiene un
efecto drama tico en el factor de reactivacion. Cuando la roca del reservorio es mas dúticl se
evidancia la prescencia de diferentes regiones con signos negativos.
Fig. 8: efecto de la rígidez en λFR . Izquierda) Yacimiento rígido. Derecha) Yacimiento dúctil.
Valores normalizados por (1−2ν)
(1−ν) .
Fig. 9: Efecto del buzamiento del yacimiento con (45 grados) en λFR . Valores normalizados
por (1−2ν)
(1−ν) .
Fig. 10: Yacimiento anticlinal, efecto de la geometria en λFR . Valores normalizados por
(1−2ν)
(1−ν) .
E. Buzamiento del yacimiento
La inclinacion del yacimiento tienen un gran efecto en λFR . La figura 9 muestra que las fallas
cerca de los flancos del yacimiento tienen grande tendencia a reactivar comparado con el caso
de un yacimento con buzamiento cero.
F. Geometría del yacimiento
Teniendo encuenta geometrias mas realistas, este ejemplo muesta el efecto de un yacimiento
con forma de anticlinal. La figura 10 muestra que las fallas cerca de los flancos tienen una
tendecia mas fuerte a reactivación en comparación con el yacimiento rectangular.
G. Yacimeinto con distribución de poro no homogénea
Todos los effectos anteriores representan situaciones con distribución de poro uniforme
dentro de la inclusión, situación que es facilmente tratable con el enfoque analítico de la teoría
de inclusiones ([8], [16]), estos casos que fueron reproducidos con la herramienta desarrollada
a manera de validación. No obstante cuando la distribución de presión de poro varía en el
espacio y el tiempo, en geometrias mas realisticas con medios heterogeneos, este analítico
enfoque es insuficiente, poniendo en evidencia la ventaja de una abordage numérica.
Mateniendo la brevedad de este documento el lector interesado en los casos con distribucion
de poro no homogénea y estrategias de producción, debe referirse a [1].
V. Conclusiones
De este trabajo es posible concluir que:
• Usando un modelamiento numérico fue mostrado que el potencial de reactivación de falla
depende fundamentalmente de: el coeficiente de fricción, la geometría del reservorio, el
angulo de buzamiento y el constraste de materiales. Asi, es realmente importante tener
encuenta estas variables en un analísis de reactivación, independientemente de cual sea
el enfoque usado.
• Diferentes estrategias de depleción, affectan la reacivación de falla en formas muy diferentes. La dirección de depleción o stress path tiene un importante rol en la evaluación de
las zonas propensas a reactivación.
• En la literatura es bastante común encontrar modelos hidrodinámicos acoplados con deformación, que no consideran un tamaño optimo para evitar los efectos cel contorno. En
este trabajo es presentado una serie de gráficos adimensionales que permiten seleccionar
las dimensiones apropiadas de las rocas circundates para evitar el efecto del contorno. Es
importante recalcar que estas figuras pueden ser usadas en simulaciónes no lineares.
• La reactivacion de falla presentada en este trabajo aplica para flujo multifásico desde el
que los terminos de acoplamiento son función de la presión ponderada por la saturación
de cada fase.
• La librería neopz permite el desarrollo de complejos algoritmos de elementos finitos que
usan el paradigma de orientación a objetos. Por tanto permite una implementación rápida
y de facil mantenimiento.
• Este trabajo muestra que LPA genera resultados más precisos, comparados con la teoría
de inclusiones y STARS CMG, software comercial de uso común en la industria del Petróleo
para problemas geomecánicos.
VI. Recomendaciones
El enfoque usado provee una base robusta para un simulador de yacimientos teniendo en
cuenta efectos poroelásticos. Las siguientes modificaciones para este trabajo son recomendadas:
• LPA no ha sido usando en un analísis inverso en contraste con un caso real.
• Extensión a la tercera dimension, flujo multifásico y acoplamiento con cambios de permeabilidad.
• Modelar la interacción mecanica entre bloques fallados y el yacimiento depletado. El
modelo actual asume una respuesta completamente reversible. Sin embargo en rocas
poco consolidadas la relación esfuerzo-deformación puede ser muy diferente, todo esto
implica analisis no lineal.
References
[1] O. Durán, “Numerical approximation of reservoir fault stability with linear poro-elasticity,” Master’s thesis, State University
of Campinas, Brazil, 2013.
[2] J. Geertsma, “Problems of rock mechanics in petroleum production engineering,” International Society for
Rock Mechanics, vol. 1, p. 10, 1966. [Online]. Available: http://www.onepetro.org/mslib/servlet/onepetropreview?id=
ISRM-1CONGRESS-1966-099
[3] P. Segall, “Stress and subsidence resulting from subsurface fluid withdrawal in the epicentral region of the 1983 coalinga
earthquake,” Journal of geophysical research, vol. 90, p. 6801, 1987.
[4] L. W. Teufel, “Effect of reservoir depletion and pore pressure drawdown on in situ stress and deformation in
the ekofisk field, north sea,” American Rock Mechanics Association, vol. 1, p. 10, 1991. [Online]. Available:
http://www.onepetro.org/mslib/servlet/onepetropreview?id=ARMA-91-063
[5] N. Morita, “Rock-property changes during reservoir compaction,” SPE Formation Evaluation (Society of Petroleum
Engineers), vol. 7, pp. 197–205, 1992. [Online]. Available: http://www.onepetro.org/mslib/servlet/onepetropreview?id=
00013099
[6] J. Grasso, “Mechanics of seismic instabilities induced by the recovery of hydrocarbons,” pure and applied geophysics,
vol. 139, pp. 507–534, 1992. [Online]. Available: http://dx.doi.org/10.1007/BF00879949
[7] J. W. Rudnicki, “Alteration of regional stress by reservoirs and other inhomogeneities: stabilizing or destabilizing?”
International Congress on Rock Mechanics, vol. 3, pp. 1629– 1637, 1999.
[8] C. D. Hawkes, “Assessing fault reactivation tendency within and surrounding porous reservoirs during fluid production
or injection,” International Journal of Rock Mechanics and Mining Sciences, vol. 46, no. 1, pp. 1 – 7, 2009. [Online].
Available: http://www.sciencedirect.com/science/article/pii/S1365160908000610
[9] D. Kosloff and R. Scott, “Finite element simulation of wilmington oil field subsidence: I. linear modelling,”
Tectonophysics, vol. 65, no. 3–4, pp. 339 – 368, 1980. [Online]. Available: http://www.sciencedirect.com/science/article/
pii/0040195180900827
[10] B. Maillot, “Numerical simulation of seismicity due to £uid injection in a brittle poroelastic medium,” International Journal
of Geophys, vol. 1, p. 10, 1999.
[11] J. Roest, “Geomechanical analysis of small earthquakes at the eleveld gas reservoir,” Rock Mechanics in Petroleum
Engineering, vol. 1, p. 8, 1994. [Online]. Available: http://www.onepetro.org/mslib/servlet/onepetropreview?id=00028097
[12] M. Ferronato, “Numerical modelling of regional faults in land subsidence prediction above gas/oil reservoirs,” International
Journal for Numerical and Analytical Methods in Geomechanics, vol. 32, no. 6, pp. 633–657, 2008. [Online]. Available:
http://dx.doi.org/10.1002/nag.640
[13] P. Segall, “Induced stresses due to fluid extraction from axisymmetric reservoirs,” pure and applied geophysics, vol. 139,
pp. 535–560, 1992. [Online]. Available: http://dx.doi.org/10.1007/BF00879950
[14] ——, “Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France,” J. Geophys. Res.,
vol. 99, pp. 15 423–15 438, 1994.
[15] M. Addis, “The influence of the reservoir stress-depletion response on the lifetime considerations of well completion
design,” SPE/ISRM Rock Mechanics in Petroleum Engineering, vol. 1, p. 9, 1998.
[16] J. D. Eshelby, “The determination of the elastic field of an ellipsoidal inclusion, and related problems,” Proceedings of
the Royal Society of London. Series A, Mathematical and Physical Sciences, vol. 241, no. 1226, pp. pp. 376–396, 1957.
[Online]. Available: http://www.jstor.org/stable/100095
[17] H. F. Wang, Theory of linear poroelasticity: with applications to geomechanics and hydrogeology, ser. Princenton seires
in geophysics. Princeton University Press, 2000, no. II.
[18] C. T. Hwang, “On solutions of plane strain consolidation problems by finite element methods,” Canadian Geotechnical
Journal, vol. 8, no. 1, pp. 109–118, 1971. [Online]. Available: http://www.nrcresearchpress.com/doi/abs/10.1139/t71-009
[19] M. Brezzi, F.; Fortin, Mixed and hybrid finite elements methods, ser. Springer series in computational mathematics.
Springer-Verlag, 1991. [Online]. Available: http://books.google.com.br/books?id=yYwZAQAAIAAJ
[20] P. A. Vermeer and A. Verruijt, “An accuracy condition for consolidation by finite elements,” International Journal
for Numerical and Analytical Methods in Geomechanics, vol. 5, no. 1, pp. 1–14, 1981. [Online]. Available:
http://dx.doi.org/10.1002/nag.1610050103
[21] J. Bear, Dynamics of fluids in porous media, ser. Environmental science series. American Elsevier Pub. Co., 1972, no.
v. 1. [Online]. Available: http://books.google.com.br/books?id=yrU-AQAAIAAJ
[22] S. C. H. Wong, “Computerized, interactive calculation of dimensionless variables,” Computers & Education, vol. 14,
no. 6, pp. 525 – 536, 1990. [Online]. Available: http://www.sciencedirect.com/science/article/pii/036013159090111J
[23] L. R. Wynne, The finite element method in the static and dynamics deformation and consolidation of porous media, 2nd ed.
John Wiley, 2000.
[24] F. Guyoton and P. Volant, “Interrelation between induced seismic instabilities and complex geological structure,”
Geophysical Research Letters, vol. 19, no. 7, pp. 705–708, 1992. [Online]. Available: http://dx.doi.org/10.1029/92GL00359
[25] L. N. Germanovich, “Fault stability inside and near a depleting petroleum reservoir,” American Rock Mechanics Association,
vol. 90, p. 12, 2004.