A posteriori error estimation and modeling of unsaturated flow in fractured porous media

This doctoral thesis focuses on three topics: (1) modeling of unsaturated flow in fractured porous media, (2) a posteriori error estimation for mixed-dimensional elliptic equations, and (3) contributions to open-source software for complex multiphysics processes in porous media. In our first contribution, following a Discrete-Fracture Matrix (DFM) approach, we propose a model where Richards' equation governs the water flow in the matrix, whereas fractures are represented as lower-dimensional open channels, naturally providing a capillary barrier to the water flow. Therefore, water in the matrix is only allowed to imbibe the fracture if the capillary barrier is overcome. When this occurs, we assume that the water inside the fracture flows downwards without resistance and, therefore, is instantaneously at hydrostatic equilibrium. This assumption can be justifiable for fractures with sufficiently large apertures, where capillary forces play no role. Mathematically, our model can be classified as a coupled PDE-ODE system of equations with variational inequalities, in which each fracture is considered a potential seepage face. Our second contribution deals with error estimation for mixed-dimensional (mD) elliptic equations, which, in particular, model single-phase flow in fractured porous media. Here, based on the theory of functional a posteriori error estimates, we derive guaranteed upper bounds for the mD primal and mD dual variables, and two-sided bounds for the mD primal-dual pair. Moreover, we improve the standard results of the functional approach by proposing four ways of estimating the residual errors based on the conservation properties of the approximations, that is, (1) no conservation, (2) subdomain conservation, (3) local conservation, and (4) pointwise conservation. This results in sharper and fully-computable bounds when mass is conserved either locally or exactly. To our knowledge, to date, no error estimates have been available for fracture networks, including fracture intersections and floating subdomains. Our last contribution is related to the development of open-source software. First, we present the implementation of a new multipoint finite-volume-based module for unsaturated poroelasticity, compatible with the Matlab Reservoir Simulation Toolbox (MRST). Second, we present a new Python-based simulation framework for multiphysics processes in fractured porous media, named PorePy. PorePy, by design, is particularly well-suited for handling mixed-dimensional geometries, and thus optimal for DFM models. The first two contributions discussed above were implemented in PorePy.Denne avhandlingen tar for seg tre emner: (1) modellering av flyt i umettet porøst medium med sprekker, (2) a posteriori feilestimater for blandet-dimensjonale elliptiske ligninger, og (3) bidrag til åpen kildekode for komplekse multifysikk-prosesser i porøse medier. I det første bidraget anvender vi en Discrete-Fracture Matrix (DFM) (Diskret-Sprekk Matrise) metode til å sette opp en modell hvor Richard's ligning modellerer vann-flyt i matrisen, og sprekkene representeres som lavere-dimensjonale åpne kanaler, som naturlig virker som kapillærbarrierer til vann-flyten. Derfor vil vann i matrisen kun få tilgang til sprekken når kapillærbarrieren blir brutt. Når det inntreffer, antar vi at vannet i sprekken flyter nedover uten motstand, og at hydrostatisk ekvilibrium derfor inntreffer øyeblikkelig. Slike antakelser kan rettferdiggjøres for sprekker med tilstrekkelig stor apertur (åpning), hvor kapillærkrefter ikke har noen innvirkning. Fra et matematisk standpunkt kan modellen klassifiseres som en sammenkoblet PDE-ODE med variasjonelle ulikheter hvor hver sprekk behandles som en filtreringsfase. Det andre bidraget tar for seg feilestimater for blandet-dimensjonale elliptiske ligninger, som modellerer en-fase flyt i porøse medier med sprekker. Her anvender vi teorien for "funksjonal a posteriori feilestimater" til å finne øvre skranker for primær og dual variablene, samt øvre og nedre skranker for primær-dual paret. Dessuten viser vi at vi kan forbedre standardresultatene fra "funksjonal a posteriori feilestimater" ved å foreslå fire måte å estimere residualfeilen basert på bevaringsegenskapene til diskretiseringen. De fire forskjellige bevaringsegenskapene er; ingen bevaringsegenskap, under- domene bevaring, lokal bevaring og punktvis bevaring. Dette fører til skarpere skranker som er mulige å beregne når masse er bevart enten lokalt, eller eksakt. Vi kjenner ikke til andre tilgjengelige feilestimater for sprekknettverk som inkluderer snitt av sprekker og sprekkrender som ligger innenfor domenets rand. Det siste bidraget omhandler utvikling av åpen kildekode. Først presenterer vi imple- menteringen av en multipunktfluks-basert modul for flyt i umettet deformerbart porøst medium som er kompatibelt med "Matlab Reservoir Simulation Toolbox"(MRST). I tillegg presenterer vi et nytt Python-basert rammeverk for simulering av multifysikkprosesser i porøse medier med sprekker, som heter PorePy. Dette rammeverket er designet for å håndtere geometrier med blandede dimensjoner og er derfor optimalt for DFM modeller. De to første bidragene i avhandlingen (nevnt over) er implementert i PorePy.Doktorgradsavhandlin

Development of robust and efficient solution strategies for coupled problems

Det er mange modeller i moderne vitenskap hvor sammenkoblingen mellom forskjellige fysiske prosesser er svært viktig. Disse finner man for eksempel i forbindelse med lagring av karbondioksid i undervannsreservoarer, flyt i kroppsvev, kreftsvulstvekst og geotermisk energiutvinning. Denne avhandlingen har to fokusområder som er knyttet til sammenkoblede modeller. Det første er å utvikle pålitelige og effektive tilnærmingsmetoder, og det andre er utviklingen av en ny modell som tar for seg flyt i et porøst medium som består av to forskjellige materialer. For tilnærmingsmetodene har det vært et spesielt fokus på splittemetoder. Dette er metoder hvor hver av de sammenkoblede modellene håndteres separat, og så itererer man mellom dem. Dette gjøres i hovedsak fordi man kan utnytte tilgjengelig teori og programvare for å løse hver undermodell svært effektivt. Ulempen er at man kan ende opp med løsningsalgoritmer for den sammenkoblede modellen som er trege, eller ikke kommer frem til noen løsning i det hele tatt. I denne avhandlingen har tre forskjellige metoder for å forbedre splittemetoder blitt utviklet for tre forskjellige sammenkoblede modeller. Den første modellen beskriver flyt gjennom deformerbart porøst medium og er kjent som Biot ligningene. For å anvende en splittemetode på denne modellen har et stabiliseringsledd blitt tilført. Dette sikrer at metoden konvergerer (kommer frem til en løsning), men dersom man ikke skalerer stabiliseringsleddet riktig kan det ta veldig lang tid. Derfor har et intervall hvor den optimale skaleringen av stabiliseringsleddet befinner seg blitt identifisert, og utfra dette presenteres det en måte å praktisk velge den riktige skaleringen på. Den andre modellen er en fasefeltmodell for sprekkpropagering. Denne modellen løses vanligvis med en splittemetode som er veldig treg, men konvergent. For å forbedre dette har en ny akselerasjonsmetode har blitt utviklet. Denne anvendes som et postprosesseringssteg til den klassiske splittemetoden, og utnytter både overrelaksering og Anderson akselerasjon. Disse to forskjellige akselerasjonsmetodene har kompatible styrker i at overrelaksering akselererer når man er langt fra løsningen (som er tilfellet når sprekken propagerer), og Anderson akselerasjon fungerer bra når man er nærme løsningen. For å veksle mellom de to metodene har et kriterium basert på residualfeilen blitt brukt. Resultatet er en pålitelig akselerasjonsmetode som alltid akselererer og ofte er svært effektiv. Det siste modellen kalles Cahn-Larché ligningene og er også en fasefeltmodell, men denne beskriver elastisitet i et medium bestående av to elastiske materialer som kan bevege seg basert på overflatespenningen mellom dem. Dette problemet er spesielt utfordrende å løse da det verken er lineært eller konvekst. For å håndtere dette har en ny måte å behandle tidsavhengigheten til det underliggende koblede problemet på blitt utviklet. Dette leder til et diskret system som er ekvivalent med et konvekst minimeringsproblem, som derfor er velegnet til å løses med de fleste numeriske optimeringsmetoder, også splittemetoder. Den nye modellen som har blitt utviklet er en utvidelse av Cahn-Larché ligningene og har fått navnet Cahn-Hilliard-Biot. Dette er fordi ligningene utgjør en fasefelt modell som beskriver flyt i et deformerbart porøst medium med to poroelastiske materialer. Disse kan forflytte seg basert på overflatespenning, elastisk spenning, og poretrykk, og det er tenkt at modellen kan anvendes i forbindelse med kreftsvulstmodellering.There are many applications where the study of coupled physical processes is of great importance. These range from the life sciences with flow in deformable human tissue to structural engineering with fracture propagation in elastic solids. In this doctoral dissertation, there is a twofold focus on coupled problems. Firstly, robust and efficient solution strategies, with a focus on iterative decoupling methods, have been applied to several coupled systems of equations. Secondly, a new thermodynamically consistent coupled system of equations is proposed. Solution strategies are developed for three different coupled problems; the quasi-static linearized Biot equations that couples flow through porous materials and elastic deformation of the solid medium, variational phase-field models for brittle fracture that couple a phase-field equation for fracture evolution with linearized elasticity, and the Cahn-Larché equations that model elastic effects in a two-phase elastic material and couples an extended Cahn-Hilliard phase-field equation and linearized elasticity. Finally, the new system of equations that is proposed models flow through a two-phase deformable porous material where the solid phase evolution is governed by interfacial forces as well as effects from both the fluid and elastic properties of the material. In the work that concerns the quasi-static linearized Biot equations, the focus is on the fixed-stress splitting scheme, which is a popular method for sequentially solving the flow and elasticity subsystems of the full model. Using such a method is beneficial as it allows for the use of readily available solvers for the subproblems; however, a stabilizing term is required for the scheme to converge. It is well known that the convergence properties of the method strongly depend on how this term is chosen, and here, the optimal choice of it is addressed both theoretically and practically. An interval where the optimal stabilization parameter lies is provided, depending on the material parameters. In addition, two different ways of optimizing the parameter are proposed. The first is a brute-force method that relies on the mesh independence of the scheme's optimal stabilization parameter, and the second is valid for low-permeable media and utilizes an equivalence between the fixed-stress splitting scheme and the modified Richardson iteration. Regarding the variational phase-field model for brittle fracture propagation, the focus is on improving the convergence properties of the most commonly used solution strategy with an acceleration method. This solution strategy relies on a staggered scheme that alternates between solving the elasticity and phase-field subproblems in an iterative way. This is known to be a robust method compared to the monolithic Newton method. However, the staggered scheme often requires many iterations to converge to satisfactory precision. The contribution of this work is to accelerate the solver through a new acceleration method that combines Anderson acceleration and over-relaxation, dynamically switching back and forth between them depending on a criterion that takes the residual evolution into account. The acceleration scheme takes advantage of the strengths of both Anderson acceleration and over-relaxation, and the fact that they are complementary when applied to this problem, resulting in a significant speed-up of the convergence. Moreover, the method is applied as a post-processing technique to the increments of the solver, and can thus be implemented with minor modifications to readily available software. The final contribution toward solution strategies for coupled problems focuses on the Cahn-Larché equations. This is a model for linearized elasticity in a medium with two elastic phases that evolve with respect to interfacial forces and elastic effects. The system couples linearized elasticity and an extended Cahn-Hilliard phase-field equation. There are several challenging features with regards to solution strategies for this system including nonlinear coupling terms, and the fourth-order term that comes from the Cahn-Hilliard subsystem. Moreover, the system is nonlinear and non-convex with respect to both the phase-field and the displacement. In this work, a new semi-implicit time discretization that extends the standard convex-concave splitting method applied to the double-well potential from the Cahn-Hilliard subsystem is proposed. The extension includes special treatment for the elastic energy, and it is shown that the resulting discrete system is equivalent to a convex minimization problem. Furthermore, an alternating minimization solver is proposed for the fully discrete system, together with a convergence proof that includes convergence rates. Through numerical experiments, it becomes evident that the newly proposed discretization method leads to a system that is far better conditioned for linearization methods than standard time discretizations. Finally, a new model for flow through a two-phase deformable porous material is proposed. The two poroelastic phases have distinct material properties, and their interface evolves according to a generalized Ginzburg–Landau energy functional. As a result, a model that extends the Cahn-Larché equations to poroelasticity is proposed, and essential coupling terms for several applications are highlighted. These include solid tumor growth, biogrout, and wood growth. Moreover, the coupled set of equations is shown to be a generalized gradient flow. This implies that the system is thermodynamically consistent and makes a toolbox of analysis and solvers available for further study of the model.Doktorgradsavhandlin

Iterative schemes for surfactant transport in porous media

In this work, we consider the transport of a surfactant in variably saturated porous media. The water flow is modelled by the Richards equations and it is fully coupled with the transport equation for the surfactant. Three linearization techniques are discussed: the Newton method, the modified Picard, and the L-scheme. Based on these, monolithic and splitting schemes are proposed and their convergence is analyzed. The performance of these schemes is illustrated on five numerical examples. For these examples, the number of iterations and the condition numbers of the linear systems emerging in each iteration are presented.publishedVersio

Towards large-scale modelling of fluid flow in fractured porous media

To date, the complexity of fractured porous media still precludes the direct incorporation of small-scale features into field-scale modelling. These features, however, can be instrumental in shaping and triggering coarsening instabilities and other forms of emergent behaviour which need to be considered on the field-scale. Here we develop numerical simulation methods for this purpose and demonstrate their improved performance in single-and two-phase flow simulations with models of fractured porous media. Material discontinuities in fractured porous media strongly influence single-and multi-phase fluid flow. When continuum methods are used to model transport across such interfaces, they smear out jump discontinuities of concentration or saturation. To overcome this drawback, we “explode” hybrid finite-element node-centred finite-volume models along these introducing complementary finite-volumes along the material interfaces. With this embedded discontinuity discretization we develop a transport scheme that realistically represents the dependent variable discontinuities arising at these interfaces. The main advantage of this new scheme is its ability to honour the flow effects that we know that these discontinuities have in physical experiments. We have also developed a new time-stepping control scheme for the transport equation. It allows the user to specify the volume fraction of the model in which he/she is prepared to relax the CFL condition. This scheme is applied in a study of the impact of fracture pattern development on solute transport. These two-dimensional simulations quantify the effect of the fractures on macro-scale dispersion in geomechanically generated fracture geometries, as opposed to stochastically generated ones. Among other insights, the results indicate that fracture density, fracture spacing, and the fracture-matrix flux ratio control anomalous mass transport in such media. We also find that it is crucial to embed discontinuities into large-scale models of heterogeneous porous media

A proof of convergence of a finite volume scheme for modified steady Richards’ equation describing transport processes in the pressing section of a paper machine

A number of water flow problems in porous media are modelled by Richards’ equation [1]. There exist a lot of different applications of this model. We are concerned with the simulation of the pressing section of a paper machine. This part of the industrial process provides the dewatering of the paper layer by the use of clothings, i.e. press felts, which absorb the water during pressing [2]. A system of nips are formed in the simplest case by rolls, which increase sheet dryness by pressing against each other (see Figure 1). A lot of theoretical studies were done for Richards’ equation (see [3], [4] and references therein). Most articles consider the case of x-independent coefficients. This simplifies the system considerably since, after Kirchhoff’s transformation of the problem, the elliptic operator becomes linear. In our case this condition is not satisfied and we have to consider nonlinear operator of second order. Moreover, all these articles are concerned with the nonstationary problem, while we are interested in the stationary case. Due to complexity of the physical process our problem has a specific feature. An additional convective term appears in our model because the porous media moves with the constant velocity through the pressing rolls. This term is zero in immobile porous media. We are not aware of papers, which deal with such kind of modified steady Richards’ problem. The goal of this paper is to obtain the stability results, to show the existence of a solution to the discrete problem, to prove the convergence of the approximate solution to the weak solution of the modified steady Richards’ equation, which describes the transport processes in the pressing section. In Section 2 we present the model which we consider. In Section 3 a numerical scheme obtained by the finite volume method is given. The main part of this paper is theoretical studies, which are given in Section 4. Section 5 presents a numerical experiment. The conclusion of this work is given in Section 6

A stable mesh-independent approach for numerical modelling of structured soils at large strains

We describe the large strain implementation of an elasto-plastic model for structured soils into G-PFEM, a code developed for geotechnical simulations using the Particle Finite Element Method. The constitutive model is appropriate for naturally structured clays, cement-improved soils and soft rocks. Structure may result in brittle behavior even in contractive paths; as a result, localized failure modes are expected in most applications. To avoid the pathological mesh-dependence that may accompany strain localization, a nonlocal reformulation of the model is employed. The resulting constitutive model is incorporated into a numerical code by means of a local explicit stress integration technique. To ensure computability this is hosted within a more general Implicit-Explicit integration scheme (IMPLEX). The good performance of these techniques is illustrated by means of element tests and boundary value problems.Peer ReviewedPostprint (author's final draft

Simulazione numerica del flusso e trasporto di contaminanti in mezzi porosi a saturazione e densità variabile

La legislazione riguardante la salvaguardia e la tutela delle risorse idriche, e tra queste le acque sotterranee, è in continua crescita in tutti i paesi industrializzati. La protezione delle acque di falda dal sovrasfruttamento e dalla contaminazione di origine diversa (rifiuti urbani e industriali, pesticidi e fertilizzanti, scorie nucleari, ecc...) richiede la previsione degli effetti indotti dalle attività umane sulla quantità e qualità delle risorse sotterranee, previsione che si può conseguire solo attraverso l'impiego di idonei modelli matematico-numerici. Un problema di stringente attualità in tutti i paesi che si affacciano sul Mediterraneo è l'inquinamento degli acquiferi costieri per intrusione di acqua di mare. La simulazione della penetrazione del cuneo salino comporta lo sviluppo di modelli accoppiati di flusso e trasporto che possono essere accuratamente ed efficientemente risolti col metodo degli elementi finiti (FEM) che viene qui implementato in un mezzo poroso tridimensionale a saturazione variabile, e che è quindi in grado di simulare sia la zona insatura (suoli superficiali) che quella satura (falde in pressione). Le non linearità che scaturiscono dall'accoppiamento e dalle leggi costitutive della permeabilità e del coefficiente di immagazzinamento nella zona insatura sono risolte con le tecniche di Picard e di Newton parziale. I modelli discreti finali linearizzati sono trattati col metodo dei gradienti coniugati opportunamente precondizionati per le matrici simmetriche di flusso (PGC) e quelle non simmetriche di trasporto (GMRES, Bi-CGSTAB, TFQMR). Le procedure descritte sono implementate nel codice FEM CODESA-3D (COupled variable DEnsity and SAturation) di cui è offerto un esempio applicativo.In the industrialized countries subsurface water resources are increasingly subject to regulations for protection from over-exploitation and from contamination arising from urban, industrial, nuclear, military, and agricultural activities. Prediction of the effects of anthropogenic impacts on water quantity and quality is an important part of proper aquifer management, and can be achieved through the use of mathematical models. As an example, seawater intrusion in coastal aquifers represents a serious environmental problem, especially in the countries of the Mediterranean basin, and can be simulated using coupled models of water flow and solute transport. Sophisticated groundwater models such as these can be accurately and effciently solved numerically via finite element discretizations of the three-dimensional porous medium. Both saturated (groundwater) and unsaturated (soil water) zones can be represented, and nonlinearities arising from storage-pressure head and conductivity-pressure head dependencies in the unsaturated zone and from coupling of the two equations can be resolved using Picard, Newton, or partial Newton methods. The resulting linearized systems of equations can be solved using a variety of preconditioned conjugate gradient-like methods applicable to symmetric and non-symmetric systems. The mathematical formulation and numerical procedures to be described form the basis of the CODESA-3D (COupled variable DEnsity and SAturation) model

Flux Postprocessing and Semi-Analytical Particle Tracking for Finite-Element-Type Models of Variably Saturated Flow in Porous Media

Aquifers are predominantly conceptualized as porous media. In real-world applications, the material properties and boundary conditions are heterogeneous requiring numerical models for an accurate description of flow and transport processes. Such numerical models yield discretized hydraulic heads and velocities on nodes or within elements. Particle tracking is a computationally advantageous and fast scheme to determine trajectories and travel times. Accurate particle tracking relies on conforming velocity fields that ensure local mass conservation in elements, and a continuous normal velocity component on element boundaries. While cell-centered finite-volume and mixed finite-element methods result in conforming velocity fields by definition, this is not the case for continuous Galerkin methods, such as the standard finite element method (FEM), and some finite-difference discretizations. Nonetheless standard FEM and also finite difference methods (FDM), formulated in finite-element terms, are often used for subsurface flow modeling because they yield a continuous approximations of hydraulic heads, and easily handle unstructured grids and material anisotropy. Acknowledging these advantages and the wide-spread use of finite-element-type simulations, the aim of this thesis is to present a novel framework for computing conforming velocity fields and accurate particle trajectories for finite-element type primal solutions of variably saturated flow in porous media. In this thesis, two different postprocessing methods to compute conforming velocity fields based on non-conforming primal solutions and semi-analytical particle tracking techniques for triangles, tetrahedra, and triangular prisms are presented. The first postprocessing method is a projection mapping a non-conforming, element-wise given velocity field onto a conforming velocity field in lowest-order Raviart-Thomas-Nédélec (RTN0) space, which meets the requirements of accurate particle tracking. The projection is based on minimizing the difference in the hydraulic gradients at the element centroids between a finite-element-type primal solution and the hydraulic gradients consistent with the RTN0 velocity field imposing element-wise mass conservation for variably saturated flow in porous media. The results of the RTN0-projection are close to those of a cell-centered finite volume method (FVM) defined for comparison and the finite-element-type primal solution. Consistency and convergence of the RTN0-projection are empirically shown for saturated flow based on a test case including hydraulic anisotropy. However, the RTN0-projection requires solving a large indefinite system of linear equations, which strongly reduces the choice of appropriate solvers, and occasionally shows numerical artifacts like velocities of wrong magnitude and unphysical direction in sub-regions of the domain in special cases. The second postprocessing technique reconstructs a cell-centered finite-volume solution from a finite-element-type primal solution of variably saturated flow in porous media to obtain conforming, mass-conservative fluxes in RTN0-space. The method is exemplified for triangular prisms because this is one of the most common elements used for catchment-scale subsurface discretization. The finite-volume flux reconstruction only solves a linear elliptic problem whose size is on the order of the number of elements, which is computationally much faster than solving the RTN0-projection or the initial Richards equation describing non-linear transient and variably saturated flow. Compared to other postprocessing schemes, the finite-volume flux reconstruction is numerically stable, fast to compute, and does not induce severe numerical artifacts when applied to heterogeneous domains with strongly varying velocities. Its advantage lies in the application to non-linear flow laws and transient problems, because it is assumed that the non-linearities are already solved by the primal solution and the transient term can be treated as a change in storage. It is shown that the results of the finite-volume flux reconstruction are close to the finite-element-type primal solution for variably saturated three-dimensional flow with heterogeneous material properties and boundary conditions. Semi-analytical particle tracking is based on element-wise analytical solutions for particle trajectories and associated travel times given a numerically approximated velocity field. This facilitates the direct computation of the spatial coordinates where the particle exits an element from its entry point and the attributes of the element-wise velocity field. In this thesis, element-wise analytical solutions are given for triangles, tetrahedra, and triangular prisms using the linear average velocity field derived from the fluid fluxes in RTN0-space. The contribution of this thesis is to provide a computational framework for approximating conforming, mass-conservative velocity fields and, based on this, semi-analytical particle tracking routines applicable to finite-element-type models of variably saturated flow in porous media on unstructured grids