Search

Constitutive laws relating fluid potentials and fluxes in a nonlinear manner are common in several porous media applications, including biological and reactive flows, poromechanics, and fracture deformation. Compared to the standard, linear Darcy’s law, such enhanced flux relations increase both the degree of nonlinearity, and, in the case of multiphysics simulations, coupling strengthbetween processes. While incorporating the nonlinearities into simulation models is thus paramount for computational efficiency, correct linearization, as is needed for incorporation in Newton’s method, is challenging from a practical perspective. The standard approach is therefore to ignore nonlinearities in the permeability during linearization. For finite volume methods, which are popular in porous media applications, complete linearization is feasible only for the simplest flux discretization, namely the two-point flux approximation. We introduce an approximated linearization scheme for finite volume methods that is exact for the two-point scheme and can be applied to more advanced and accurate discretizations, exemplified herein by a multi-point flux stencil. We test the new method for both nonlinear porous media flow and several multiphysics simulations. Our results show that the new linearization consistently outperforms the standard approach. Moreover our scheme achieves asymptotic second order convergence of the Newton iterations, in contrast to the linear convergence obtained with the standard approach.

Discretization of flow in fractured porous media commonly lead to large systems of linear equations that require dedicated solvers. In this work, we develop an efficient linear solver and its practical implementation for mixed-dimensional scalar elliptic problems. We design an effective preconditioner based on approximate block factorization and algebraic multigrid techniques. Numerical results on benchmarks with complex fracture structures demonstrate the effectiveness of the proposed linear solver and its robustness with respect to different physical and discretization parameters.

As a consequence of shearing, wing cracks can emerge from pre-existing fractures. The process involves the interaction of sliding of the existing fracture surfaces and the tensile material failure that creates wing cracks. This work devises a numerical model to investigate how wing cracks emerge, propagate and connect pre-existing fractures under shear processes. A mathematical and numerical model for wing crack propagation based on linear elastic fracture mechanics that also accounts for fracture contact mechanics is presented. Computational efficiency is ensured by an adaptive remeshing technique. The numerical model is verified and validated through a comparison of the analytical and experimental results. Additional numerical examples illustrate the performance of the method for complex test cases where wing-cracks develop for multiple pre-existing and interacting fractures. publishedVersion

Subsurface storage of CO2 is predicted to rise exponentially in response to the increasing levels of CO2 in the atmosphere. Large-scale CO2 injections into the subsurface require understanding of the potential for fluid flow through faults to mitigate risk of leakage. Here, we study how to obtain effective permeability of deformation bands in the damage zone of faults. Deformation bands are relatively small, low permeability features that can have a significant effect on flow dynamics, however, the discrepancy of scales is a challenge for field-scale simulation. A new analytical upscaling model is proposed in order to overcome some of the shortcomings of conventional upscaling approaches for heterogeneous porous media. The new model captures the fine-scale impact of deformation bands on fluid flow in the near-fault region, and can be derived from knowledge of large-scale fault properties. To test the accuracy of the model it is compared to fine-scale numerical simulations that explicitly include individual deformation bands. For a wide range of different stochastically generated deformation bands networks, the upscaling model shows improved estimate of effective permeability compared to conventional upscaling approaches. By applying the upscaling model to a full-field simulation of the Smeaheia storage site in the North Sea, we show that deformation bands with a permeability contrast higher than three orders of magnitude may act as an extra layer of protection from fluid flow through faults. submittedVersion publishedVersion

A standard practice used in the industry to discretizing the gravity term in the two‐phase Darcy flow equations is to apply an upwind strategy. In this paper, we show that this can give a persistent unphysical flux field and an incorrect pressure distribution. As a solution to this problem, we present a new consistent discretization of flow, termed Gravitationally Consistent Multipoint Flux Approximation (GCMPFA), which is valid for both single‐ and two‐phase flows. The discretization is based on the idea that the gravitational term in the flow equations is treated as part of the discrete flux operator and not as a right‐hand side. Here, the traditional formulation representing pressure as a potential is extended to the case including gravity by introducing an additional set of right‐hand side to the local linear system solved in the MPFA construction, thus obtaining an expression of the fluxes in terms of jumps in cell‐center gravities. Numerical examples showing the convergence of the method are provided for both single‐ and two‐phase flows. For two‐phase flow, we show how our new method is capable of eliminating the unphysical fluxes arising when using a standard upwind scheme, thus converging to the correct pressure distribution. publishedVersion

This article reports and describes the data related to the paper “Conforming, non-conforming and non-matching discretization couplings in discrete fracture network simulations” (Fumagalli et al., 2019). The data provided include a set of geometrical input data of Discrete Fracture Networks (DFNs) and a set of simulation results. The geometrical data describe the geometry of fracture networks of increasing complexity. These data also include the geometry of a DFN extruded from a real fracture outcrop in Western Norway. Simulation results are obtained using several different numerical schemes and provide convergence history, plots over line and upscaled output quantities related to the various considered geometries.

Various phenomena in the subsurface are characterised by the interplay between deforming structures such as fractures and coupled thermal, hydraulic and mechanical processes. Simulation of subsurface dynamics can provide valuable phenomenological understanding, but requires models which faithfully represent the dynamics involved; these models therefore are themselves highly complex. This paper presents a mixed-dimensional thermo-hydro-mechanical model designed to capture the process–structure interplay using a discrete–fracture–matrix framework. It incorporates tightly coupled thermo-hydro-mechanical processes based on balance laws for momentum, mass and energy in subdomains representing the matrix and the lower-dimensional fractures and fracture intersections. The deformation of explicitly represented fractures is modelled by contact mechanics relations and a Coulomb friction law, with a novel formulation consistently integrating fracture dilation in the governing equations. The model is discretised using multi-point finite volume methods for the balance equations and a semismooth Newton scheme for the contact conditions and is implemented in the open-source fracture simulation toolbox PorePy. Finally, simulation studies demonstrate the model’s convergence, investigate process–structure coupling effects, explore different fracture dilation models and show an application of the model to stimulation and long-term cooling of a three-dimensional geothermal reservoir. publishedVersion

When modeling fluid flow in fractured reservoirs, it is common to represent the fractures as lower-dimensional inclusions embedded in the host medium. Existing discretizations of flow in porous media with thin inclusions assume that the principal directions of the inclusion permeability tensor are aligned with the inclusion orientation. While this modeling assumption works well with tensile fractures, it may fail in the context of faults, where the damage zone surrounding the main slip surface may introduce anisotropy that is not aligned with the main fault orientation. In this article, we introduce a generalized dimensional reduced model which preserves full-tensor permeability effects also in the out-of-plane direction of the inclusion. The governing equations of flow for the lower-dimensional objects are obtained through vertical averaging. We present a framework for discretization of the resulting mixed-dimensional problem, aimed at easy adaptation of existing simulation tools. We give numerical examples that show the failure of existing formulations when applied to anisotropic faulted porous media, and go on to show the convergence of our method in both two-dimensional and three-dimensional. publishedVersion

In this chapter, we give a unified introduction to the MPFA- and MPSA-type finite volume methods for Darcy flow and poro-elasticity, applicable to general polyhedral grids. This leads to a more systematic perspective of these methods than has been exposed in previous texts, and we therefore refer to this discretization family as the MPxA methods. We apply this MPxA framework to also define a consistent finite-volume discretization of thermo-poro-elasticity. The present chapter introduces the general theory and state-of-the-art of MPFA-type methods, leaving the more technical results to the provided references. We close the chapter by a section containing applications to problems with complex geometries and non-linear physics.

Development of models and dedicated numerical methods for dynamics in fractured rocks is an active research field, with research moving towards increasingly advanced process couplings and complex fracture networks. The inclusion of coupled processes in simulation models is challenged by the high aspect ratio of the fractures, the complex geometry of fracture networks, and the crucial impact of processes that completely change characteristics on the fracture-rock interface. This paper provides a general discussion of design principles for introducing fractures in simulators, and defines a framework for integrated modeling, discretization, and computer implementation. The framework is implemented in the open-source simulation software PorePy, which can serve as a flexible prototyping tool for multiphysics problems in fractured rocks. Based on a representation of the fractures and their intersections as lower-dimensional objects, we discuss data structures for mixed-dimensional grids, formulation of multiphysics problems, and discretizations that utilize existing software. We further present a Python implementation of these concepts in the PorePy open-source software tool, which is aimed at coupled simulation of flow and transport in three-dimensional fractured reservoirs as well as deformation of fractures and the reservoir in general. We present validation by benchmarks for flow, poroelasticity, and fracture deformation in porous media. The flexibility of the framework is then illustrated by simulations of non-linearly coupled flow and transport and of injection-driven deformation of fractures. All results can be reproduced by openly available simulation scripts.

In this work, we consider compressible single-phase flow problems in a porous medium containing a fracture. In the fracture, a nonlinear pressure-velocity relation is prescribed. Using a non-overlapping domain decomposition procedure, we reformulate the global problem into a nonlinear interface problem. We then introduce two new algorithms that are able to efficiently handle the nonlinearity and the coupling between the fracture and the matrix, both based on linearization by the so-called L-scheme. The first algorithm, named MoLDD, uses the L-scheme to resolve for the nonlinearity, requiring at each iteration to solve the dimensional coupling via a domain decomposition approach. The second algorithm, called ItLDD, uses a sequential approach in which the dimensional coupling is part of the linearization iterations. For both algorithms, the computations are reduced only to the fracture by precomputing, in an offline phase, a multiscale flux basis (the linear Robin-to-Neumann codimensional map), that represent the flux exchange between the fracture and the matrix. We present extensive theoretical findings X and in particular, t. The stability and the convergence of both schemes are obtained, where user-given parameters are optimized to minimize the number of iterations. Examples on two important fracture models are computed with the library PorePy and agree with the developed theory. publishedVersion

The hydraulic stimulation of the well RN-34 at the Reykjanes geothermal field in Iceland caused increased seismic activity near the well. Here, we use this as a case study for investigation on how seismic analysis can be combined with physics-based simulation studies to further understand injection-induced fault reactivation. The work presents new analysis of the seismic data combined with application of a recent simulation software for modeling of coupled hydromechanical processes and fault deformation caused by fluid injection. The simulation model incorporates an explicit model of the fault network based on geological characterization combined with insights from seismic analysis. The 3D faulted reservoir model is then calibrated based on injection data. Despite limited data, the work shows how seismic interpretations can be used in developing simulation models and, reciprocally, how the modeling can add to the seismic interpretations in analysis of dynamics.

International audience; Flow in fractured porous media occurs in the earth's subsurface, in biological tissues, and in man-made materials. Fractures have a dominating influence on flow processes, and the last decade has seen an extensive development of models and numerical methods that explicitly account for their presence. To support these developments, we present a portfolio of four benchmark cases for single-phase flow in three-dimensional fractured porous media. The cases are specifically designed to test the methods' capabilities in handling various complexities common to the geometrical structures of fracture networks. Based on an open call for participation, results obtained with 17 numerical methods were collected. This paper presents the underlying mathematical model, an overview of the features of the participating numerical methods, and their performance in solving the benchmark cases.

Convection-driven cooling in porous media influences thermo-poro-mechanical stresses, thereby causing deformation. These processes are strongly influenced by the presence of fractures, which dominate flow and heat transfer. At the same time, the fractures deform and propagate in response to changes in the stress state. Mathematically, the model governing the physics is tightly coupled and must account for the strong discontinuities introduced by the fractures. Over the last decade, and motivated by a number of porous media applications, research into such coupled models has advanced modelling of processes in porous media substantially. Building on this effort, this work presents a novel model that couples fracture flow and heat transfer and deformation and propagation of fractures with flow, heat transfer and thermo-poroelasticity in the matrix. The model is based on explicit representation of fractures in the porous medium and discretised using multi-point finite volume methods. Frictional contact and non-penetration conditions for the fractures are handled through active set methods, while a propagation criterion based on stress intensity factors governs fracture extension. Considering both forced and natural convection processes, numerical results show the intricate nature of thermo-poromechanical fracture deformation and propagation.

This paper presents a model where thermo-hydro-mechanical processes are coupled to a deformation model for preexisting fractures. The model is formulated within a discrete-fracture-matrix framework where the rock matrix and the fractures are considered as individual subdomains, and interaction between them takes place on the matrix-fracture interfaces. A finite volume discretisation implemented in the simulation toolbox PorePy is presented and applied in a simulation showcasing the effects of the different mechanisms on fracture deformation governed by contact mechanics, as well as their different timescales.

Upscaling of flow from pore to Darcy scale is a long-standing research field within flow in porous media. It is well known that non-linearities can occur in near-well regions and high-porosity or fractured media. At the same time, the upscaled non-linear effects associated with high flow rates are hard to quantify a priori with single-scale models. Advances in pore-scale imaging combined with increased computational have made flow simulations in small pore-scale domain feasible, but computations on domains larger than at most a few centimeters are still elusive. In this work, we present a multiscale simulation framework that automatically adapts to non-linear effects as they arise. We formulate a control volume heterogeneous multiscale method (CVHMM) by coupling of a Darcy-scale control volume method with a constitutive relation that is captured based on the fine-scale physics. While the CVHMM formulation works with arbitrary upscaled laws, we emphasize its ability to be applied in fully discrete multiscale context, in particular when a finite element solver is used for solving Navier-Stokes equations on the fine-scale pore geometry. Previous versions of CVHMM are consistent only when the coarse grid is aligned with the upscaled permeability. Herein, we generalize CVHMM by introducing a new coarse solver, thus significantly improving the applicability of the method. The presented method is applied to study flow in near-well regions, as well as media with fractures and irregular grain shapes. The examples show that the method successfully copes with general grids and pore geometries and handles flows with varying degree of non-linearities even outside the domain of applicability of classical upscaled models. In terms of computational efficiency, the method seamlessly localizes computations to regions where non-linear effects are important.

Over the last decades, finite-volume discretisations for flow in porous media have been extended to handle situations where fractures dominate the flow. Successful discretisations have been based on the discrete fracture-matrix models to yield mass conservative methods capable of explicitly incorporating the impact of fractures and their geometry. When combined with a hybrid-dimensional formulation, two central concerns are the restrictions arising from small cell sizes at fracture intersections and the coupling between fractures and matrix. Focusing on these aspects, we demonstrate how finite-volume methods can be efficiently extended to handle fractures, providing generalisations of previous work. We address the finite-volume methods applying a general hierarchical formulation, facilitating implementation with extensive code reuse and providing a natural framework for coupling of different subdomains. Furthermore, we demonstrate how a Schur complement technique may be used to obtain a robust and versatile method for fracture intersection cell elimination. We investigate the accuracy of the proposed elimination method through a series of numerical simulations in 3D and 2D. The simulations, performed on fractured domains containing permeability heterogeneity and anisotropy, also demonstrate the flexibility of the hierarchical framework.

Simulating the deformation of fractured media requires the coupling of different models for the deformation of fractures and the formation surrounding them. We consider a cell-centered finite-volume approach, termed the multi-point stress approximation (MPSA) method, which is developed in order to discretize coupled flow and mechanical deformation in the subsurface. Within the MPSA framework, we consider fractures as co-dimension one inclusions in the domain, with the fracture surfaces represented as line pairs in 2D (face pairs in 3D) that displace relative to each other. Fracture deformation is coupled to that of the surrounding domain through internal boundary conditions. This approach is natural within the finite-volume framework, where tractions are defined on surfaces of the grid. The MPSA method is capable of modeling deformation, considering open and closed fractures with complex and nonlinear relationships governing the displacements and tractions at the fracture surfaces. We validate our proposed approach using both problems, for which analytical solutions are available, and more complex benchmark problems, including comparison with a finite-element discretization.

In simulation of fluid injection in fractured geothermal reservoirs, the characteristics of the physical processes are severely affected by the local occurence of connected fractures. To resolve th ...

Understanding the controlling mechanisms underlying injection-induced seismicity is important for optimizing reservoir productivity and addressing seismicity-related concerns related to hydraulic stimulation in Enhanced Geothermal Systems. Hydraulic stimulation enhances permeability through elevated pressures, which cause normal deformations and the shear slip of preexisting fractures. Previous experiments indicate that fracture deformation in the normal direction reverses as the pressure decreases, e.g., at the end of stimulation. We hypothesize that this normal closure of fractures enhances pressure propagation away from the injection region and significantly increases the potential for postinjection seismicity. To test this hypothesis, hydraulic stimulation is modeled by numerically coupling flow in the fractures and matrix, fracture deformation, and matrix deformation for a synthetic reservoir in which the flow and mechanics are strongly affected by a complex three-dimensional fracture network. The role of the normal closure of fractures is verified by comparing simulations conducted with and without the normal closure effect. publishedVersion

We introduce a new cell-centered finite volume discretization for elasticity with weakly enforced symmetry of the stress tensor. The method is motivated by the need for robust discretization methods for deformation and flow in porous media, and falls in the category of multi-point stress approximations (MPSA). By enforcing symmetry weakly, the resulting method has additional flexibility beyond previous MPSA methods. This allows for a construction of a method which is applicable to all grid types, and in particular the method amends a crucial shortcoming in previous MPSA methods for simplex grids. By formulating the method as a discrete variational problem, we prove convergence of the new method for a wide range of problems, with conditions that can be verified at the time of discretization. We present the first set of comprehensive numerical tests for the MPSA methods in three dimensions, covering Cartesian and simplex grids, with both heterogeneous and nearly incompressible media. The tests show that the new method consistently is second order convergent in displacement, despite being lowest order, with a rate that mostly is between 1 and 2 for stresses. The results further show that the new method is more robust and computationally cheaper than previous MPSA methods.

The effect of heterogeneities induced by highly permeable fracture networks on viscous miscible fingering in porous media is examined using high-resolution numerical simulations. We consider the planar injection of a less viscous fluid into a two-dimensional fractured porous medium which is saturated with a more viscous fluid. This problem contains two sets of fundamentally different preferential flow regimes; the first is caused by the viscous fingering and the second is due to the permeability contrasts between the fractures and rock matrix. We study the transition from a regime where the flow is dominated by the viscous instabilities, to a regime where the heterogeneities induced by the fractures define the flow paths. We find that fractures greatly affect the viscous fingering, even for small permeability differences between the rock matrix and the fractures. The interaction between the viscosity contrast and permeability contrast causes channeling of the less viscous fluid through the fractures and back to the rock. This channeling stabilizes the displacement front in the rock matrix, and the viscous fingering ceases for the higher permeability contrast. Several different fracture geometries are considered, and we observe a complex interplay between the geometries and unstable flow. While we find that the most important dimensionless number determining the effect of the fracture network is a weighted ratio of the permeability of the fractures and the permeability of the rock matrix, the exact point for the cross-over regime is highly dependent on the geometry of the fracture network., Comment: To reproduce simulations, see "R. L. Berge, I. Berre, E. Keilegavlen, and J. M. Nordbotten. Viscous fingering in fractured porous media. Zenodo. doi: 10.5281/zenodo.3249931, June 2019"

We study a heuristic, core-scale model for the transport of polymer particles in a two phase (oil and water) porous medium. We are motivated by recent experimental observations which report increased oil recovery when polymers are injected after the initial waterflood. The recovery mechanism is believed to be microscopic diversion of the flow, where injected particles can accumulate in narrow pore throats and clog it, in a process known as a log-jamming effect. The blockage of the narrow pore channels lead to a microscopic diversion of the water flow, causing a redistribution of the local pressure, which again can lead to the mobilization of trapped oil, enhancing its recovery. Our objective herein is to develop a core-scale model that is consistent with the observed production profiles. We show that previously obtained experimental results can be qualitatively explained by a simple two-phase flow model with an additional transport equation for the polymer particles. A key aspect of the formulation is that the microscopic heterogeneity of the rock and a dynamic altering of the permeability must be taken into account in the rate equations., 20 pages, 9 Figures Submitted to Transport in Porous Media

The last decade has seen a strong increase of research into flows in fractured porous media, mainly related to subsurface processes, but also in materials science and biological applications. Connected fractures totally dominate flow-patterns, and their representation is therefore a critical part in model design. Due to the fracture's characteristics as approximately planar discontinuities with an extreme size to width ratio, they challenge standard macroscale mathematical and numerical modeling of flow based on averaging. Thus, over the last decades, various, and also fundamentally different, approaches have been developed. This paper reviews common conceptual models and discretization approaches for flow in fractured porous media, with an emphasis on the dominating effects the fractures have on flow processes. In this context, the paper discuss the tight connection between physical and mathematical modeling and simulation approaches. Extensions and research challenges related to transport, multi-phase flow and fluid-solid interaction are also commented on.

We study the transport of inertial particles in water flow in porous media. Our interest lies in understanding the accumulation of particles including the possibility of clogging. We propose that accumulation can be a result of hydrodynamic effects: the tortuous paths of the porous medium generate regions of dominating strain, which favour the accumulation of particles. Numerical simulations show that essentially two accumulation regimes are identified: for low and for high flow velocities. When particles accumulate at the entrance of a pore throat (high-velocity region), a clog is formed. This significantly modifies the flow, as the partial blockage of the pore causes a local redistribution of pressure, which diverts the upstream water flow into neighbouring pores. Moreover, we show that accumulation in high velocity regions occurs in heterogeneous media, but not in homogeneous media, where we refer to homogeneity with respect to the distribution of the pore throat diameters. MAEK, EK and KS are supported by Equinor through the Akademia agreement. The work of AM is partially supported via MultiBarr (KK Foundation, project nr. 20180036). FAR was partially supported by the Norwegian Academy of Science and Letters and Equinor through the VISTA AdaSim project #6367. ISP was supported by the Research Foundation-Flanders (FWO) through the Odysseus programme (project GOG1316N) and by Equinor through the Akademia agreement. KK work was partially funded by the Research Council of Norway (NFR) Projects 811716 LAB2FIELD and 810857 MICAP. (FWO) through the Odysseus programme (project GOG1316N); Research Council of Norway (NFR) Projects 811716 LAB2FIELD and 810857 MICAP

Models that involve coupled dynamics in a mixed-dimensional geometry are of increasing interest in several applications. Here, we describe the development of a simulation model for flow in fractured porous media, where the fractures and their intersections form a hierarchy of interacting subdomains. We discuss the implementation of a simulation framework, with an emphasis on reuse of existing discretization tools for mono-dimensional problems. The key ingredients are the representation of the mixed-dimensional geometry as a graph, which allows for convenient discretization and data storage, and a non-intrusive coupling of dimensions via boundary conditions and source terms. This approach is applicable for a wide class of mixed-dimensional problems. We show simulation results for a flow problem in a three-dimensional fracture geometry, applying both finite volume and virtual finite element discretizations. acceptedVersion

Fluid-induced slip of fractures is characterized by strong multiphysics couplings. Three physical processes are considered: Flow, rock deformation and fracture deformation. The fractures are represented as lower-dimensional objects embedded in a three-dimensional domain. Fluid is modeled as slightly compressible, and flow in both fractures and matrix is accounted for. The deformation of rock is inherently different from the deformation of fractures; thus, two different models are needed to describe the mechanical deformation of the rock. The medium surrounding the fractures is modeled as a linear elastic material, while the slip of fractures is modeled as a contact problem, governed by a static-dynamic friction model. We present an iterative scheme for solving the non-linear set of equations that arise from the models, and suggest how the step parameter in this scheme should depend on the shear modulus and mesh size. acceptedVersion

Simulations of fluid flow in naturally fractured rocks have implications for several subsurface applications, including energy storage and extraction, and waste disposal. We are interested in flow in discrete fracture networks, which explicitly represent flow in fracture surfaces, but ignore the impact of the surrounding host rock. Fracture networks, generated from observations or stochastic simulations, will contain intersections of arbitrary length, and intersection lines can further cross, forming a highly complex geometry. As the flow exchange between fractures, thus in the network, takes place in these intersections, an adequate representation of the geometry is critical for simulation accuracy. In practice, the intersection dynamics must be handled by a combination of the simulation grid, which may or may not resolve the intersection lines, and the numerical methods applied on the grid. In this work, we review different classes of numerical approaches proposed in recent years, covering both methods that conform to the grid, and non-matching cases. Specific methods considered herein include finite elements, mixed and virtual elements and control volume methods. We expose our methods to an extensive set of test cases, ranging from artificial geometries designed to test difficult configurations, to a network extruded from a real fracture outcrop. The main outcome is guidances for choice of simulation models and numerical discretization with a trade off on the computational cost and solution accuracy. acceptedVersion

Numerical simulations of brittle fracture using phase-field approaches often employ a discrete approximation framework that applies the same order of interpolation for the displacement and phase-field variables. In particular, the use of linear finite elements to discretize both stress equilibrium and phase-field equations is widespread in the literature. However, the use of P 1 Lagrange shape functions to model the phase-field is not optimal, as the latter contains cusps for fully developed cracks. These should in turn occur at locations corresponding to Gauss points of the associated FE model for the mechanics. Such a feature is challenging to reproduce accurately with low order elements, and element sizes must consequently be made very small relative to the phase-field regularization parameter in order to achieve convergence of results with respect to the mesh. In this paper, we combine a standard linear FE discretization of stress equilibrium with a cell-centered finite volume approximation of the phase-field evolution equation based on the two-point flux approximation constructed over the same simplex mesh. Compared to a pure FE formulation utilizing linear elements, the proposed framework results in looser restrictions on mesh refinement with respect to the phase-field length scale. This ability to employ coarser meshes relative to the traditional implementation allows for significant reductions on computational cost, as demonstrated in several numerical examples.

Porous siliciclastic reservoirs are known to contain structural heterogeneities such as deformation bands, which fall below current seismic resolution and which generally cannot be explicitly represented in reservoir models because of the prohibitively high computational cost. In this study, we built a reservoir model to evaluate fluid flow across a contractionally folded unit containing deformation bands (the Navajo Sandstone in the San Rafael Reef monocline, Utah). Using field data, geometric relationships, and auxiliary computational techniques, we upscale deformation bands to capture flow effects in the large-scale structure, running simulations with variable scenarios of permeability contrast between host rock and deformation bands. Our simulations show that pervasive deformation band arrays (such as the ones present in the monocline) have effects when the contrast of permeability between them and the host rock is of at least three orders of magnitude, delaying water breakthrough and enhancing sweep; in long-term production, this results in larger final produced volumes and higher total recovery. Because of the wide range of deformation band permeabilities used in this study, our findings can be of importance for the prediction of flow and optimization of production strategies in comparable traps and reservoirs. Additionally, auxiliary computational techniques and geometric relationships such as the ones presented in this study can significantly improve the incorporation of small-scale features with strong scale gap into conventional sized reservoirs.

Objective: Chronic kidney disease (CKD) is a serious medical condition characterized by gradual loss of kidney function. Early detection and diagnosis is mandatory for adequate therapy and prognostic improvement. Hence, in the current pilot study we explore the use of image registration methods for detecting renal morphologic changes in patients with CKD. Methods: Ten healthy volunteers and nine patients with presumed CKD underwent dynamic T1 weighted imaging without contrast agent. From real and simulated dynamic time series, kidney deformation fields were estimated using a poroelastic deformation model. From the deformation fields several quantitative parameters reflecting pressure gradients, and volumetric and shear deformations were computed. Eight of the patients also underwent a kidney biopsy as a gold standard. Results: We found that the absolute deformation, normalized volume changes, as well as pressure gradients correlated significantly with arteriosclerosis from biopsy assessments. Furthermore, our results indicate that current image registration methodologies are lacking sensitivity to recover mild changes in tissue stiffness. Conclusion: Image registration applied to dynamic time series correlated with structural renal changes and should be further explored as a tool for invasive measurements of arteriosclerosis. Significance: Under the assumption that the proposed framework can be further developed in terms of sensitivity and specificity, it can provide clinicians with a non-invasive tool of high spatial coverage available for characterization of arteriosclerosis and potentially other pathological changes observed in chronic kidney disease.

In this paper, we introduce a mortar-based approach to discretizing flow in fractured porous media, which we term the mixed-dimensional flux coupling scheme. Our formulation is agnostic to the discretizations used to discretize the fluid flow equations in the porous medium and in the fractures, and as such it represents a unified approach to integrated fractured geometries into any existing discretization framework. In particular, several existing discretization approaches for fractured porous media can be seen as special instances of the approach proposed herein. We provide an abstract stability theory for our approach, which provides explicit guidance into the grids used to discretize the fractures and the porous medium, as dependent on discretization methods chosen for the respective domains. The theoretical results are sustained by numerical examples, wherein we utilize our framework to simulate flow in 2D and 3D fractured media using control volume methods (both two-point and multi-point flux), Lagrangian finite element methods, mixed finite element methods, and virtual element methods. As expected, regardless of the ambient methods chosen, our approach leads to stable and convergent discretizations for the fractured problems considered, within the limits of the discretization schemes.

Solving realistic problems related to flow in porous media to desired accuracy may be prohibitively expensive with available computing resources. Multiscale effects and nonlinearities in the governing equations are among the most important contributors to this situation. Hence, developing methods that handle these features better is essential in order to be able to solve the problems more efficiently. Focus has until recently largely been on preconditioners for linearized problems. This article proposes a two-scale nonlinear preconditioning technique for flow problems in porous media that allows for incorporating physical intuition directly in the preconditioner. By assuming a certain dominant physical process, this technique will resemble upscaling in the equilibrium limit, with the computational benefits that follow. In this study, the method is established as a preconditioner with good scalability properties for challenging problems regardless of dominant physics, thus laying the foundation for further studies with physical information in the preconditioner. publishedVersion

Phase-field approaches to fracture based on energy minimization principles have been rapidly gaining popularity in recent years, and are particularly well-suited for simulating crack initiation and growth in complex fracture networks. In the phase-field framework, the surface energy associated with crack formation is calculated by evaluating a functional defined in terms of a scalar order parameter and its gradients, which in turn describe the fractures in a diffuse sense following a prescribed regularization length scale. Imposing stationarity of the total energy leads to a coupled system of partial differential equations, one enforcing stress equilibrium and another governing phase-field evolution. The two equations are coupled through an energy degradation function that models the loss of stiffness in the bulk material as it undergoes damage. In the present work, we introduce a new parametric family of degradation functions aimed at increasing the accuracy of phase-field models in predicting critical loads associated with crack nucleation as well as the propagation of existing fractures. An additional goal is the preservation of linear elastic response in the bulk material prior to fracture. Through the analysis of several numerical examples, we demonstrate the superiority of the proposed family of functions to the classical quadratic degradation function that is used most often in the literature., 33 pages, 30 figures

The process of initial ice formation in brine is a highly complex problem. In this paper, we propose a mathematical model that captures the dynamics of nucleation and development of ice inclusions in brine. The primary emphasis is on the interaction between ice growth and salt diffusion, subject to external forcing provided by temperature. Within this setting two freezing regimes are identified, depending on the rate of change of the temperature: a slow freezing regime where a continuous ice domain is formed; and a fast freezing regime where recurrent nucleation appears within the fluid domain. The second regime is of primary interest, as it leads to fractal-like ice structures. We analyse the critical threshold between the slow and fast regimes by identifying the explicit rates of external temperature control that lead to self-similar salt-concentration profiles in the fluid domains. Subsequent heuristic analysis provides estimates of the characteristic length scales of the fluid domains depending on the time-variation of the temperature. The analysis is confirmed by numerical simulations. S.A. thanks A. Kvashchuk for useful discussions during preparation of the manuscript and B. Balakin for helpful references. J.M.N. visited the Isaac Newton programme 'Melt in the Mantle' at Cambridge University while working on this manuscript. I.S.P. acknowledges support from the Akademia grant of Statoil and from the Research Foundation - Flanders FWO through the the Odysseus programme grant G0G1316N. The authors also acknowledge the feedback from the editor and reviewers that helped to improve the presentation of the paper.

Conventional vertically integrated models for CO2 storage usually adopt a vertical equilibrium (VE) assumption, which states that due to strong buoyancy, CO2 and brine segregate quickly, so that the fluids can be assumed to have essentially hydrostatic pressure distributions in the vertical direction. However, the VE assumption is inappropriate when the time scale of fluid segregation is not small relative to the simulation time. By casting the vertically integrated equations into a multiscale framework, a new vertically integrated model can be developed that relaxes the VE assumption, thereby allowing vertical dynamics to be modeled explicitly. The model maintains much of the computational efficiency of vertical integration while allowing a much wider range of problems to be modeled. Numerical tests of the new model, using injection scenarios with typical parameter sets, show excellent behavior of the new approach for homogeneous geologic formations.

Discrete fracture models are an attractive alternative to upscaled models for flow in fractured media, as they provide a more accurate representation of the flow characteristics. A major challenge in discrete fracture simulation is to overcome the large computational cost associated with resolving the individual fractures in large-scale simulations. In this work, two characteristics of the fractured porous media are utilized to construct efficient preconditioners for the discretized flow equations. First, the preconditioners are tailored to the fracture geometry and presumed flow properties so that the dominant features are well represented there. This assures good scalability of the preconditioners in terms of problem size and permeability contrast. For fracture dominated problems, numerical examples show that such geometric preconditioners are comparable or preferable when compared to state-of-the-art algebraic multigrid preconditioners. The robustness of the physics-based preconditioner for less favorable fracture conditions is further demonstrated by a systematic degradation of the fracture hierarchy. Second, the preconditioners are physics preserving in the sense that conservative fluxes can be computed even for an inexact pressure solutions. This facilitates a scheme where accuracy in the linear solver can be traded for efficiency by terminating the iterative solvers based on error estimates, and without sacrificing basic physical modeling principles. With the combination of these two properties a novel preconditioner is obtained which bridges the gap between multiscale approximations and iterative linear solvers.

The accurate description of fluid flow and transport in fractured porous media is of paramount importance to capture the macroscopic behaviour of an oil reservoir, a geothermal system, or a CO2 sequestration site, to name few applications. The construction of accurate simulation model for flow in fractures is challenging due to the high ratios between a fracture's length and width, which makes modeling by lower-dimensional manifolds a natural option. In this paper we present a mixed-dimensional Darcy problem able to describe pressure and Darcy velocity in all the dimensions, i.e. in the rock matrix, in the fractures, and in their intersections. Moreover, we present a mixed-dimensional transport problem which, given the Darcy velocity, describes coupled advection and diffusion of a passive scalar into the fractured porous media. The approach can handle both conducting and blocking fractures. Our computational grids are created by coarsening of simplex tessellations that conform to the fractures surfaces. An accurate choice of the discrete approximation of the previous model, by virtual finite element and finite volume, allows us to simulate complex problem with a good balance in term of accuracy and computational cost. We illustrate the performance of our method by comparing to benchmark studies for two-dimensional fractured porous media, as well as a complex three-dimensional fracture geometry.

Cataclastic deformation bands, which are common in sandstone reservoirs and which may negatively affect fluid flow, are generally associated with notable thickness variations. It has been suggested previously that such thickness variations represent an important control on how deformation bands affect fluid flow. The effects of such thickness variations are tested in this study though statistical analysis and fluid flow simulation of an array of cataclastic deformation bands in Cretaceous sandstones in in the Bassin de Sud-Est in Provence, France. Spatial outcrop data are statistically analyzed for incorporation in flow simulation models, and numerical simulations are used to investigate the effects of notable thickness variations on how the deformation bands influence effective permeability and flow dynamics. A suite of simulations is performed using a combination of fine-scale and coarse-scale grids, revealing that the effective permeability of the simulated reservoir is reduced by a factor of 15–25. More interestingly, the simulations further demonstrated that, as compared to the overall effect of the deformation band array on fluid flow, thickness variations along the bands proved to have negligible effects only. Thus, our simulations indicate that the configuration and connectivity of the deformation bands, together with the permeability contrast between the bands and the host rock and the mean band thickness, are the most important controls on the effective permeability. Our findings represent new insight into the influence of deformation bands on fluid flow in subsurface aquifers and reservoirs, indicating that thickness variations of individual deformation bands are of less significance than previously thought.

The multiscale control-volume methods for solving problems involving flow in porous media have gained much interest during the last decade. Recasting these methods in an algebraic framework allows one to consider them as preconditioners for iterative solvers. Despite intense research on the 2D formulation, few results have been shown for 3D, where indeed the performance of multiscale methods deteriorates. The interpretation of multiscale methods as vertex based domain decomposition methods, which are non-scalable for 3D domain decomposition problems, allows us to understand this loss of performance. We propose a generalized framework based on auxiliary variables on the coarse scale. These are enrichments the coarse scale, which can be selected to improve the interpolation onto the fine scale. Where the existing coarse scale basis functions are designed to capture local sub-scale heterogeneities, the auxiliary variables are aimed at better capturing non-local effects resulting from non-linear behavior of the pressure field. The auxiliary coarse nodes fits into the framework of mass-conservative domain-decomposition (MCDD) preconditioners, allowing us to construct, as special cases, both the traditional (vertex based) multiscale methods as well as their wire basket generalization. submittedVersion

Domain decomposition (DD) methods, such as the additive Schwarz method, are almost exclusively applied to linearized equations. In the context of nonlinear problems, these linear systems appear as part of a Newton iteration. However, applying DD methods directly to the original nonlinear problem has some attractive features, most notably that the Newton iterations now solve local problems, and thus are expected to be simpler. Furthermore, strong, local nonlinearities may to a less extent affect the numerical algorithm. For linear problems, DD can be applied both as an iterative solver or as a preconditioner. For nonlinear problems, it has until recently only been understood how to use DD as a solver. This article offers a systematic study of domain decomposition strategies in the context of nonlinear porous-medium flow problems. The study thus compares four different approaches, which represents DD applied both as a solver and preconditioner, to both the linearized and nonlinear equations. Our model equations are those obtained from a fully implicit discretization of immiscible two-phase flow in heterogeneous porous media. In particular we emphasize the case of nonlinear preconditioning, an algorithm that to our knowledge so far has not been studied nor implemented for flow in porous media. Our results show that the novel algorithm is up to 75% faster than the standard algorithm for the most challenging problems for a moderate number of subdomains.

Discrete fracture networks is a key ingredient in the simulation of physical processes which involve fluid flow in the underground, when the surrounding rock matrix is considered impervious. In this paper we present two different models to compute the pressure field and Darcy velocity in the system. The first allows a normal flow out of a fracture at the intersections, while the second grants also a tangential flow along the intersections. For the numerical discretization, we use the mixed virtual finite element method as it is known to handle grid elements of, almost, any arbitrary shape. The flexibility of the discretization allows us to loosen the requirements on grid construction, and thus significantly simplify the flow discretization compared to traditional discrete fracture network models. A coarsening algorithm, from the algebraic multigrid literature, is also considered to further speed up the computation. The performance of the method is validated by numerical experiments.