An adaptive multiresolution scheme is proposed for the numerical solution of a spatially two-dimensional model of sedimentation of suspensions of small solid particles dispersed in a viscous fluid. This model consists in a version of the Stokes equations for incompressible fluid flow coupled with a hyperbolic conservation law for the local solids concentration. We study the process in an inclined, rectangular closed vessel, a configuration that gives rise a well-known increase of settling rates (compared with a vertical vessel) known as the "Boycott effect". Sharp fronts and discontinuities in the concentration field are typical features of sedimentation phenomena. This solution behavior calls for locally refined meshes to concentrate computational effort on zones of strong variation. The spatial discretization presented herein is naturally based on a finite volume (FV) formulation for the Stokes problem including a pressure stabilization technique, while a Godunov-type scheme endowed with a fully adaptive multiresolution (MR) technique is applied to capture the evolution of the concentration field, which in addition induces an important speed-up of CPU time and savings in memory requirements. Numerical simulations illustrate that the proposed scheme is robust and allows for substantial reductions in computational effort while the computations remain accurate and stable.

This paper examines two-phase flow in porous media with heterogeneous capillary pressure functions. This problem has received very little attention in the literature, and constitutes a challenge for numerical discretization, since saturation discontinuities arise at the interface between the different homogeneous regions in the domain. As a motivation we first consider a one-dimensional model problem, for which a semi-analytical solution is known, and examine some different finite-volume approximations. A standard scheme based on harmonic averaging of the absolute permeability, and which possesses the important property of being pressure continuous at the discrete level, is found to converge and gives the best numerical results. In order to investigate two-dimensional flow phenomena by a robust and accurate numerical scheme, a recent multi point flux approximation scheme, which is also pressure continuous at the discrete level, is then extended to account for two-phase flow, and is used to discretize the two-phase flow pressure equation in a fractional flow formulation well suited for capillary heterogenity. The corresponding saturation equation is discretized by a second-order central upwind scheme. Some numerical examples are presented in order to illustrate the significance of capillary pressure heterogeneity in two-dimensional two-phase flow, using both structured quadrilateral and unstructured triangular grids.

This paper presents numerical results on the development of compositional fluid mixing simulators in porous media. These simulators integrate geological processes (source rock maturation, hydrocarbon generation, migration, charge/filling, etc.) and reservoir processes (fluid mixing through Darcy's flow, advection, and diffusion, gravity segregation, etc.). The model governing equations are written with a proper choice of solution variables so that numerical mass conservation is preserved for all chemical components. The approximation procedure uses the finite volume method for space discretization, the backward Euler scheme in time, and an adaptive time stepping technique. The traditional simulator for solving the isothermal gravity/chemical equilibrium problem is deduced as a special example of the simulators presented here. Extensive numerical experiments are given to show segregation and instability effects for multiple components.

We consider the application of a perfectly matched layer (PML) technique applied in Cartesian geometry to approximate solutions of the electromagnetic wave (Maxwell) scattering problem in the frequency domain. The PML is viewed as a complex coordinate shift ("stretching") and leads to a variable complex coefficient equation for the electric field posed on an infinite domain, the complement of a bounded scatterer. The use of Cartesian geometry leads to a PML operator with simple coefficients, although, still complex symmetric (non-Hermitian). The PML reformulation results in a problem which preserves the original solution inside the PML layer while decaying exponentially outside. The rapid decay of the PML solution suggests truncation to a bounded domain with a convenient outer boundary condition and subsequent finite element approximation (for the truncated problem).
For suitably defined Cartesian PML layers, we prove existence and uniqueness of the solutions to the infinite domain and truncated domain PML equations provided that the truncated domain is sufficiently large. We show that the PML reformulation preserves the solution in the layer while decaying exponentially outside of the layer. Our approach is to develop variational stability for the infinite domain electromagnetic wave scattering PML problem from that for the acoustic wave (Helmholtz) scattering PML problem given in [12]. The stability and exponential convergence of the truncated PML problem is then proved using the decay properties of solutions of the infinite domain problem. Although, we do not provide a complete analysis of the resulting finite element approximation, we believe that such an analysis should be possible using the techniques in [6].

We consider the flow of two-phases in a porous medium and propose a modified version of the fractional flow model for incompressible, two-phase flow based on a Helmholtz regularization of the Darcy phase velocities. We show the existence of global-in-time entropy solutions for this model with suitable assumptions on the boundary conditions. Numerical experiments demonstrating the approximation of the classical two-phase flow equations with the new model are presented.

The analysis of the Multi Point Flux Approximation (MPFA) method has so far relied on the possibility of seeing it as a mixed finite element method for which the convergence is then established. This type of analysis has been successfully applied to triangles and quadrilaterals, also in the case of rough meshes. The MPFA method has however much in common with another well known conservative method: the mimetic finite difference method. We propose to formulate the MPFA O-method in a mimetic finite difference framework, in order to extend the proof of convergence to polyhedral meshes. The formulation is useful to see the close relationship between the two different methods and to see how the differences lead to different strengths. We pay special attention to the assumption needed for proving convergence by examining various cases in the section dedicated to numerical tests.

The authors formulated in [32] a multipoint flux mixed finite element method that reduces to a cell-centered pressure system on general quadrilaterals and hexahedra for elliptic equations arising in subsurface flow problems. In addition they showed that a special quadrature rule yields $\mathcal{O}(h)$ convergence for face fluxes on distorted hexahedra. Here a first order local velocity postprocessing procedure using these face fluxes is developed and analyzed. The algorithm involves solving a 3$\times$3 system on each element and utilizes an enhanced mixed finite element space introduced by Falk, Gatto, and Monk [18]. Computational results verifying the theory are demonstrated.

Coupled flow- and rock mechanics simulations are necessary to achieve sufficient understanding and reliable production forecasts in many reservoirs, especially those containing weak or moderate strength rock. Unfortunately these runs are in general significantly more demanding with respect to computing times than stand-alone flow simulations. A scheme is presented whereby the number of needed rock mechanics simulations in such a setting can be reduced to a minimum. The scheme is based on constructing optimal input for flow simulations from a few rock mechanics runs. Results obtained with the scheme are at least as accurate as traditional coupled runs, but computations are considerably faster, often as much as two orders of magnitude.

The pumping function of the heart is driven by an electrical wave traversing the cardiac muscle in a well-organized manner. Perturbations to this wave are referred to as arrhythmias. Such arrhythmias can, under unfortunate circumstances, turn into fibrillation, which is often lethal. The only known therapy for fibrillation is a strong electrical shock. This process, referred to as defibrillation, is routinely used in clinical practice. Despite the importance of this procedure and the fact that it is used frequently, the reasons for defibrillation's effectiveness are not fully understood. For instance, theoretical estimates of the shock strength needed to defibrillate are much higher than what is actually used in practice. Several authors have pointed out that, in theoretical models, the strength of the shock can be decreased if the cardiac tissue is modeled as a heterogeneous substrate. In this paper, we address this issue using the bidomain model and the Courtemanche ionic model; we also consider a linear approximation of the Courtemanche model here. We present analytical considerations showing that for the linear model, the necessary shock strength needed to achieve defibrillation (defined in terms of a sufficiently strong change of the resting state) decreases as a function of an increasing perturbation of the intracellular conductivities. Qualitatively, these theoretical results compare well with computations based on the Courtemanche model. The analysis is based on an energy estimate of the difference between the linear solution of the bidomain system and the equilibrium solution. The estimate states that the difference between the linear solution and the equilibrium solution is bounded in terms of the shock strength. Since defibrillation can be defined in terms of a certain deviation from equilibrium, we use the energy estimate to derive a necessary condition for the shock strength.

Fractional diffusion equations provide an adequate and accurate description of transport processes that exhibit anomalous diffusion that cannot be modeled accurately by classical second-order diffusion equations. However, numerical discretizations of fractional diffusion equations yield full coefficient matrices, which require a computational operation of $O(N^3)$ per time step and a memory of $O(N^2)$ for a problem of size $N$. In this paper we develop a fast second-order finite difference method for space-fractional diffusion equations, which only requires memory of $O(N)$ and computational work of $O(N log^2 N)$. Numerical experiments show the utility of the method.

Natural fractures reside in various subsurface formations and are at various length scales with different intensities. Fluid flow in fractures, in matrix and between matrix and fractures are following different flow physics. It is thus a great challenge for efficiently modeling and simulation of fluid flow in fractured media due to the multi-scale and multi-physics nature of the flow processes.
Traditional dual-porosity and dual-permeability approach represents fractures and matrix as different continuum. The transfer functions or shape factors are derived to couple the fluid flow in matrix and fractures. The dual-porosity and dual-permeability model can be viewed as a multi-scale method and the transfer functions are used to propagate fine-scale information to the coarse-scale reservoir simulation. In this paper, we perform a detailed study to better understand the optimal way to propagate the fracture information to the coarse-scale model based on the detailed fracture characterization at fine-scale.
The Discrete Fracture Modeling (DFM) approach is used to represent each fracture individually and explicitly. The multiple sub-region (MSR) method is previously used for upscaling calculations based on fine-scale flow solution by finite volume method on the DFM. The MSR method is the most appropriate upscaling procedure for connected fracture network but not for disconnected fractures. In this paper, we propose an adaptive hybrid multi-scale approach that combines MSR and DFM adaptively for upscaling calculation for complex fractured subsurface formations that usually involve both connected fracture network and disconnected fractures. The numerical results suggest that adaptive hybrid multi-scale approach can provide accurate upscaling results for flow in a complicated geological system.

We study the dynamics of two-phase flow with gravity and point out three different transport mechanisms: non-cyclic advection, solenoidal advection, and gravity segregation. Each term has specific mathematical properties that can be exploited by specialized numerical methods. We argue that to develop effective operator splitting methods, one needs to understand the interplay between these three mechanisms for the problem at hand.

For multi-phase flow in porous media, the fluids will experience reduced conductivity due to fluid-fluid interactions. This effect is modeled by the so-called relative permeability. The relative permeability is commonly assumed to be a scalar quantity, even though justifications for this modeling choice seldom are presented. In this paper, we show that the relative permeability can yield preferential flow directions, and that furthermore, these directions may vary as a function of the fluids present. To model these effects properly, the relative permeability should be considered a tensor. As shown in the paper, standard numerical methods cannot simulate tensor relative permeabilities. We therefore propose two new methods to remedy the situation. One scheme can be seen as an attempt to amend traditional methods to handle tensor relative permeabilities. The second scheme is a consistent approach to simulate the tensors. We also present numerical simulations of cases where upscaling leads to tensor relative permeabilities. The results show that the new methods are indeed capable of capturing tensor effects, moreover, they are in very good agreement with fine scale reference solutions.

Multi-phase flow in porous media is most commonly modeled by adding a saturation-dependent, scalar relative permeability into the Darcy equation. However, in the general case anisotropically structured heterogeneities result in anisotropy of upscaled parameters, which not only depends on the solid structure but also on fluid-fluid or fluid-fluid-solid interaction. We present a method for modeling of incompressible, isothermal, immiscible two-phase flow, which accounts for anisotropic absolute as well as relative permeabilities. It combines multipoint flux approximation (MPFA) with an appropriate upwinding strategy in the framework of a sequential solution algorithm. Different tests demonstrate the capabilities of the method and motivate the relevance of anisotropic relative permeabilities. Therefore, a porous medium is chosen, which is heterogeneous but isotropic on a fine scale and for which averaged homogeneous but anisotropic parameters are known. Comparison shows that the anisotropy in the large-scale parameters is well accounted for by the method and agrees with the anisotropic distribution behavior of the fine-scale solution. This is demonstrated for both the advection dominated as well as the diffusion dominated case. Further, it is shown that off-diagonal entries in the relative permeability tensor can have a significant influence on the fluid distribution.

Equations describing flow in porous media averaged to allow for lateral spatial variability but integrated over the vertical dimension are derived from pore scale equations. Under conditions of vertical equilibrium, the equations are simplified and employed to describe migration of $\rm{CO}_2$ injected into an aquifer of variable thickness. The numerical model based on the vertical equilibrium equations is shown to agree well with a fully three-dimensional model. Trapping of $\rm{CO}_2$ in undulations at the top of the aquifer is shown to retard $\rm{CO}_2$ migration.