Numerical models for the large-scale simulation of fault and fracture mechanics

The possible activation of pre-existing faults and the generation of new fractures in the subsurface may play a critical role in several fields of great social interest, such as the management and the exploitation of groundwater resources, especially in arid areas, the hydrocarbon recovery and storage, and the monitoring of the seismic activity in the Earth’s crust. The sliding and/or opening of a fault can create preferential leakage paths for the pore fluid escape, causing a matter of great concern in the process of storing fluids and hydrocarbons underground. The most challenging effect connected to a fault activation is the possible earthquake triggering. Many earthquakes associated with the production and injection of fluids have been recently reported. Similar issues arise also in the development of unconventional hydrocarbon reservoirs, that has recently experienced a dramatic increase thanks to the deployment of the “fracking” technology, which is based on the massive generation of fractures through the injection of fluids at high pressures. The use of this technique in densely populated areas has raised a large scientific debate on the possible connected environmental risks. The over-exploitation of fresh aquifers in arid regions has caused the generation of significant ground fissures. In this thesis, a novel formulation based on the use of Lagrange multipliers has been developed for the stable and robust numerical modeling of fault mechanics. A fault or fracture is simulated as a pair of inner surfaces included in a 3D geological formation where Lagrange multipliers are used to prescribe the contact constraints. The standard variational formulation of the contact problem with Lagrange multipliers is modified to take into account the energy dissipated by the frictional work along the activated fault portion. This term is computed by making use of the principle of maximum plastic dissipation, whose application defines the direction of the limiting shear stress vector. The novel approach has been verified against analytical solutions and applied in a number of real-world problems. In particular, we test the novel approach in four cases: (i) mechanics of two adjacent blocks, to investigate the numerical properties of the algorithm; (ii-iii) ground fractures due to groundwater withdrawal, with different geometries; (iv) fault reactivation in an underground reservoir subject to primary production and Underground Gas Storage cycles. The results are analyzed and commented. In the fourth case, the possible magnitude of the seismic events triggered by fault reactivation is computed, in order to evaluate whether underground human activities may generate seismicity. The application of the fault model to large-scale problems gives rise to a set of sparse discrete systems of linearized equations with a generalized non-symmetric saddle point structure. The second part of this thesis is devoted to the development of efficient algorithms for the iterative solution of this kind of system. We focus on a preconditioning technique, denoted as “constraint preconditioning”, which exploits the native block structure of the Jacobian. The quality and performance of the preconditioner relies on two steps: (i) the preconditioning of the leading block and (ii) the Schur complement computation. In this work, novel preconditioning techniques for the leading block based on a multilevel framework are developed and tested. The main idea behind the multilevel preconditioner is to improve the quality of the factorized approximate inverses borrowing the scheme of incomplete factorizations, thus introducing some sequentially in perfectly parallelizable algorithms. The proposed approach is robust, from a theoretical point of view, and very efficient in parallel environment. As to the latter point, i.e. the Schur complement computation, it can be done with the aid of different approximations. The main difference is whether the Jacobian is symmetrized or not. The computation can be founded on the FSAI approximation of the leading block inverse or on a physically-based block diagonal block algorithm. The Schur complement must be inverted, thus other possibilities come in. The approximate Schur complement can be inverted through FSAI, if symmetric, or an incomplete factorization, if non-symmetric, but it can also be solved exactly, thanks to a direct solver. The performances of the proposed algorithms are finally investigated and discussed in a set of real-world numerical examples.

L’eventuale attivazione di faglie preesistenti e la generazione di nuove fratture nel sottosuolo può svolgere un ruolo cruciale in diversi campi di grande interesse sociale, come la gestione e lo sfruttamento delle risorse idriche sotterranee, in particolare in aree aride, l’estrazione e lo stoccaggio di idrocarburi e il monitoraggio dell’attività sismica nella crosta terrestre. Lo scorrimento e/o l’apertura di una faglia può creare percorsi di dispersione preferenziali per la fuoriuscita del fluido interstiziale causando gravi problemi nel processo di stoccaggio sotterraneo di fluidi e idrocarburi. L’effetto più problematico connesso ad un’attivazione di faglia è l’eventuale innesco di un terremoto. Recentemente, sono stati riportati molti terremoti associati all’estrazione e all’iniezione di fluidi. Problemi simili si verificano anche nella coltivazione di giacimenti di idrocarburi non convenzionali, che recentemente hanno registrato un drastico aumento grazie alla diffusione della tecnologia chiamata “fracking”, che si basa sulla creazione di numerose fratture attraverso l’iniezione di fluidi ad alta pressione. L’uso di questa tecnica in aree densamente popolate ha sollevato un ampio dibattito scientifico sui possibili rischi ambientali ad essa connessi. L’eccessivo sfruttamento di acquiferi in regioni aride ha generato notevoli fratture nel terreno. In questa tesi, è stata sviluppata una nuova formulazione basata sui moltiplicatori di Lagrange, in grado di offrire una robusta modellazione numerica della meccanica delle faglie. Una faglia o una frattura viene simulata come una coppia di superfici interne inserite in una formazione geologica tridimensionale, in cui i moltiplicatori di Lagrange vengono utilizzati per imporre le condizioni di contatto. La classica formulazione variazionale del problema di contatto con i moltiplicatori di Lagrange è modificata per includere l’energia dissipata dal lavoro delle forze d’attrito nella parte di faglia attivata. Questo quantità viene calcolata utilizzando il principio della massima dissipazione plastica, la cui applicazione definisce la direzione del vettore di tensione tangenziale in caso di scorrimento. Il nuovo approccio è stato validato grazie a soluzioni analitiche e utilizzato in diversi problemi derivanti da applicazioni reali. In particolare, il nuovo approccio è stato testato in quattro casi: (i) la meccanica di due blocchi adiacenti, per analizzare il comportamento dell’algoritmo da un punto di vista numerico; (ii-iii) fratture del sottosuolo causate da estrazione d’acqua subsuperficiale, con differenti geometrie; (iv) riattivazione di faglie in un giacimento soggetto a produzione primaria e cicli di stoccaggio di gas. I risultati sono stati analizzati e commentati. Nel quarto caso, si è determinata la possibile magnitudo relative ad eventi provocati dalla riattivazione delle faglie, per capire se le attività umane nel sottosuolo possano provocare fenomeni sismici. L’applicazione del modello di faglia a problemi di grande scala genera un insieme di sistemi discreti e sparsi di equazioni linearizzate, con una struttura tipo punto sella non simmetrica generalizzata. La seconda parte di questa tesi è rivolta allo sviluppo di algoritmi efficienti per la soluzione iterativa di sistemi di questo tipo. In particolare, si utilizza la tecnica detta “constraint preconditioning”, che sfrutta la suddivisione in blocchi dello jacobiano. La qualità del precondizionatore si basa su due fasi: (i) il precondizionamento del blocco principale e (ii) il calcolo del complemento di Schur. In questo lavoro, sono state sviluppate e testate nuove tecniche di precondizionamento per il blocco principale, basate su un framework multilevello. L’idea di base di questo approccio è il miglioramento delle inverse approssimate grazie allo schema delle fattorizzate incomplete, introducendo, quindi, una sequenzialità in algoritmi perfettamente paralleli. La tecnica presentata è robusta da un punto di vista teorico e molto efficiente in ambiente parallelo. Per quanto riguarda l’ultimo punto, ossia il calcolo del complemento di Schur, sono state proposte diverse approssimazioni. La differenza principale risiede nel considerare o meno lo jacobiano simmetrizzato. Il calcolo può essere basato sulla FSAI, per approssimare l’inversa del blocco principale, o su un approccio diagonale a blocchi, che è legato alla caratteristiche fisiche del problema. Il complemento di Schur deve essere invertito, quindi altre scelte entrano in gioco. Si può utilizzare la FSAI, se la matrice è simmetrica, o una fattorizzata incompleta, se non simmetrica, ma è possibile anche risolvere esattamente il complemento di Schur, grazie ad un solutore diretto. Infine, le performance degli algoritmi proposti sono analizzate e discusse su un insieme di problemi numerici derivanti da applicazioni reali.