Search

Physics-informed neural networks (PINNs) have received significant attention as a unified framework for forward, inverse, and surrogate modeling of problems governed by partial differential equations (PDEs). Training PINNs for forward problems, however, pose significant challenges, mainly because of the complex non-convex and multi-objective loss function. In this work, we present a PINN approach to solving the equations of coupled flow and deformation in porous media for both single-phase and multiphase flow. To this end, we construct the solution space using multi-layer neural networks. Due to the dynamics of the problem, we find that incorporating multiple differential relations into the loss function results in an unstable optimization problem, meaning that sometimes it converges to the trivial null solution, other times it moves very far from the expected solution. We report a dimensionless form of the coupled governing equations that we find most favourable to the optimizer. Additionally, we propose a sequential training approach based on the stress-split algorithms of poromechanics. Notably, we find that sequential training based on stress-split performs well for different problems, while the classical strain-split algorithm shows an unstable behaviour similar to what is reported in the context of finite element solvers. We use the approach to solve benchmark problems of poroelasticity, including Mandel's consolidation problem, Barry-Mercer's injection-production problem, and a reference two-phase drainage problem. The Python-SciANN codes reproducing the results reported in this manuscript will be made publicly available at https://github.com/sciann/sciann-applications.

The injection and storage of freshwater in saline aquifers for the purpose of managed aquifer recharge is an important technology that can help ensure sustainable water resources. As a result of the density difference between the injected freshwater and ambient saline groundwater, the pressure field is coupled to the spatial salinity distribution, and therefore experiences transient changes. The effect of variable density can be quantified by the mixed convection ratio, which is a ratio between the strength of two convection processes: free convection due to the density differences and forced convection due to hydraulic gradients. We combine a density-dependent flow and transport simulator with an ensemble Kalman filter (EnKF) to analyze the effects of freshwater injection rates on the value-of-information of transient pressure data for saline aquifer characterization. The EnKF is applied to sequentially estimate heterogeneous aquifer permeability fields using real-time pressure data. The performance of the permeability estimation is analyzed in terms of the accuracy and the uncertainty of the estimated permeability fields as well as the predictability of breakthrough curve arrival times in a realistic push-pull setting. This study demonstrates that injecting fluids at a rate that balances the two characteristic convections can maximize the value of pressure data for saline aquifer characterization. Keywords: Managed aquifer recharge; Density-dependent flow; Inverse modeling; Ensemble Kalman filter; Value of information; Permeability estimation

We investigate tracer transport on random discrete fracture networks that are characterized by the statistics of the fracture geometry and hydraulic conductivity. While it is well known that tracer transport through fractured media can be anomalous and particle injection modes can have major impact on dispersion, the incorporation of injection modes into effective transport modeling has remained an open issue. The fundamental reason behind this challenge is that-even if the Eulerian fluid velocity is steady-the Lagrangian velocity distribution experienced by tracer particles evolves with time from its initial distribution, which is dictated by the injection mode, to a stationary velocity distribution. We quantify this evolution by a Markov model for particle velocities that are equidistantly sampled along trajectories. This stochastic approach allows for the systematic incorporation of the initial velocity distribution and quantifies the interplay between velocity distribution and spatial and temporal correlation. The proposed spatial Markov model is characterized by the initial velocity distribution, which is determined by the particle injection mode, the stationary Lagrangian velocity distribution, which is derived from the Eulerian velocity distribution, and the spatial velocity correlation length, which is related to the characteristic fracture length. This effective model leads to a time-domain random walk for the evolution of particle positions and velocities, whose joint distribution follows a Boltzmann equation. Finally, we demonstrate that the proposed model can successfully predict anomalous transport through discrete fracture networks with different levels of heterogeneity and arbitrary tracer injection modes. © 2017 Elsevier Ltd., PKK and SL acknowledge a grant (16AWMP- B066761-04) from the AWMP Program funded by the Ministry of Land, Infrastructure and Transport of the Korean government and the support from Future Research Program (2E27030) funded by the Korea Institute of Science and Technology (KIST). PKK and RJ acknowledge a MISTI Global Seed Funds award. MD acknowledges the support of the European Research Council (ERC) through the project MHetScale (617511). TLB acknowledges the support of European Research Council (ERC) through the project Re- activeFronts (648377). RJ acknowledges the support of the US Department of Energy through a DOE Early Career Award (grant DE-SC0009286). The data to reproduce the work can be obtained from the corresponding author.

It is well documented that cells can migrate in response to gradients in stiffness (durotaxis) and gradients in strain (tensotaxis) in the underlying substrate. Understanding the potential physical mechanisms at play during this motion has motivated recent efforts to unravel the role of surface tension in the interaction between droplets and soft solids. Here, we present a multiphysics phase-field model of fluid–solid interaction, which allows us to isolate the effects of strain gradients—something difficult to achieve in experiments. Our high-fidelity numerical simulations in two and three dimensions elucidate the physics of tensotaxis, and show how localized forces in a soft substrate can be used to move and merge droplets deposited on it.

We present the application of a class of deep learning, known as Physics Informed Neural Networks (PINN), to inversion and surrogate modeling in solid mechanics. We explain how to incorporate the momentum balance and constitutive relations into PINN, and explore in detail the application to linear elasticity, and illustrate its extension to nonlinear problems through an example that showcases von Mises elastoplasticity. While common PINN algorithms are based on training one deep neural network (DNN), we propose a multi-network model that results in more accurate representation of the field variables. To validate the model, we test the framework on synthetic data generated from analytical and numerical reference solutions. We study convergence of the PINN model, and show that Isogeometric Analysis (IGA) results in superior accuracy and convergence characteristics compared with classic low-order Finite Element Method (FEM). We also show the applicability of the framework for transfer learning, and find vastly accelerated convergence during network re-training. Finally, we find that honoring the physics leads to improved robustness: when trained only on a few parameters, we find that the PINN model can accurately predict the solution for a wide range of parameters new to the network—thus pointing to an important application of this framework to sensitivity analysis and surrogate modeling.

We report the emergence of anomalous (non-Fickian) transport through a rough-walled fracture as a result of increasing normal stress on the fracture. We show that the origin of this anomalous transport behavior can be traced to the emergence of a heterogeneous flow field dominated by preferential channels and stagnation zones, as a result of the larger number of contacts in a highly stressed fracture. We show that the velocity distribution determines the late-time scaling of particle spreading, and velocity correlation determines the magnitude of spreading and the transition time from the initial ballistic regime to the asymptotic anomalous behavior. We also propose a spatial Markov model that reproduces the transport behavior at the scale of the entire fracture with only three physical parameters. Our results point to a heretofore unrecognized link between geomechanics and particle transport in fractured media.

We propose a thermodynamic approach to modeling unsaturated flow in porous media, where the liquid saturation is understood as the state variable. The free energy functional is designed as a symmetric expansion of the traditional capillary energy density in Richards equation, therefore removing ambiguities on the interpretation of the higher-order term in the model equation. The proposed definition renders a formulation that leads naturally to an entropy function of the system, and we show that the model describes an entropy-increasing process for an isolated system. The new formulation reproduces gravity fingering during infiltration in soil. We show that the nonlinear and singular structure of the capillary pinning function in the fourth-order term plays a fundamental role in the behavior and stability of infiltration fronts, promoting front pinning and the persistence of fingered infiltration at relatively large flux ratios.

The density driven convection phenomenon is expected to have a significant and positive role in CO 2 geological storage capacity and safety. But predictions on reservoir time and space scales are difficult to validate because data are generally sparse and will only be useful for a small part of the relevant time period. Laboratory scale data are valuable to validate the numerical models. In this paper we focus on the comparison of experimental and numerical determination of CO 2 mass transfer in a laboratory experiment. We developed an experimental protocol for the determination of density-driven mass transfer of CO 2 in water-saturated Hele-Shaw cells with different apertures. We used a CCD camera to capture images of the initiation of density-driven convection caused by dissolution of CO 2 in water and the subsequent development of convective fingers. The visualization of the phenomenon allowed consistently stopping the experiment when dissolved CO 2 first reached the bottom of the cell. We determined the total mass of dissolved CO 2 during the experiment using a catalytic combustion-based total carbon analyzer (TC-analyzer). This experimental procedure was repeated several times for uncertainty analysis. Thus a combination of quantitative and qualitative experimental results for the same Hele-Shaw cell configuration was obtained for validation of corresponding numerical simulation results. A numerical simulation of the phenomenon was carried out using the STOMP-WCS simulator. We found that in order to accurately simulate numerically the phenomenon occurring in the Hele-Shaw cell, existent variations in the cell apertures should be taken into account. Thus we observed a good agreement between the experimental and numerical results in terms of total dissolved CO 2 mass, timescale of mass transport and morphology of the convection fingers. In addition correlations are obtained between total dissolved CO 2 mass, arrival time of dissolved CO 2 to bottom of the cell, and the Rayleigh number.

Estimating the carbon dioxide (CO[subscript 2]) storage capacity of deep saline aquifers is important for identifying those most suitable for sequestration, and for planning the future development of CO[subscript 2] storage projects. Currently, capacity estimates are highly uncertain due in part to uncertainty in the dominant constraint on capacity: both the pressure buildup from injection and the space available to trap CO[subscript 2] have been identified as constraints, but have not been rigorously compared to determine the conditions under which each is more limiting. In this study, we evaluate their relative importance in an idealized aquifer using simple, but dynamic models of how pressure rises during injection and how CO[subscript 2] becomes trapped in the pore space. We show that there exists a crossover injection duration, T[subscript c], below which pressure constraints dominate, but above which the CO[subscript 2] migration becomes the more limiting constraint. We illustrate this behavior by applying the models to the Fox Hills Sandstone., United States. Dept. of Energy (Grant DE-FE0002041), MIT Energy Initiative, Reed Research, MIT Martin Family Society of Fellows for Sustainability

We experimentally studied the displacement of a viscous wetting fluid (water) by an inviscid non-wetting fluid (air) injected at the bottom of a vertical Hele-Shaw cell filled with glass microbeads. In order to cover a wide parameter space, the permeability of the porous medium was varied by using different bead size ranges and diverse air flow rates were generated by means of a syringe pump. A LED light table was used to back illuminate the experimental cell, allowing a high speed camera to capture images of the drainage process at equal time intervals. The invasion occurred in intermittent bursts. Image processing of the bursts and fractal analysis showed successive transitions from capillary invasion to viscous fingering to fracturing during the same experiment, dependent on the medium permeability, the air injection flow rate, and the vertical position in the cell. The interplay between the capillary, viscous and gravity forces determines the nature of the invasion pattern and the transitions, from capillary invasion to viscous fingering with decreasing fluid pressure on one hand and from viscous fingering to fracturing with decreasing effective overburden pressure on the other hand.

Coupling between fluid flow and mechanical deformation in porous media plays a critical role in geologic storage of CO 2 One of the key issues in simulation of CO 2 sequestration is the ability to describe the mechanical and hydraulic behavior of faults, and the influence of the stress tensor and change in pressure on fault slip. Here, we present a new computational approach to model coupled multiphase flow and geomechanics of faulted reservoirs. We represent faults as surfaces embedded in a three-dimensional medium by using zero-thickness interface elements to accurately model fault slip under dynamically evolving fluid pressure and fault strength. We incorporate the effect of fluid pressures from multiphase flow in the mechanical stability of faults by defining a fault pressure. We employ a rigorous formulation of nonlinear multiphase geomechanics based on the increment in mass of fluid phases, instead of the change in porosity. Our nonlinear formulation is required to properly model systems with high compressibility or strong capillarity, as can be the case for geologic CO 2 sequestration. To account for the effect of surface stresses along fluid-fluid interfaces, we use the equivalent pore pressure in the definition of multiphase effective stress. We develop a numerical simulation tool by coupling a multiphase flow simulator with a mechanics simulator, using the unconditionally stable fixed-stress scheme for a computationally efficient sequential solution of two-way coupling between flow and geomechanics. We validate our modeling approach using several synthetic test cases that illustrate the onset and evolution of earthquakes from fluid injection.

During geologic storage of carbon dioxide (CO[subscript 2]), trapping of the buoyant CO[subscript 2] after injection is essential in order to minimize the risk of leakage into shallower formations through a fracture or abandoned well. Models for the subsurface behavior of the CO[subscript 2] are useful for the design, implementation, and long-term monitoring of injection sites, but traditional reservoir-simulation tools are currently unable to resolve the impact of small-scale trapping processes on fluid flow at the scale of a geologic basin. Here, we study the impact of solubility trapping from convective dissolution on the up-dip migration of a buoyant gravity current in a sloping aquifer. To do so, we conduct high-resolution numerical simulations of the gravity current that forms from a pair of miscible analogue fluids. Our simulations fully resolve the dense, sinking fingers that drive the convective dissolution process. We analyze the dynamics of the dissolution flux along the moving CO[subscript 2]–brine interface, including its decay as dissolved buoyant fluid accumulates beneath the buoyant current. We show that the dynamics of the dissolution flux and the macroscopic features of the migrating current can be captured with an upscaled sharp-interface model., Seventh Framework Programme (European Commission) (CO2-MATE Project PIOF-GA-2009-253678), Seventh Framework Programme (European Commission) (Project PANACEA Grant Agreement 282900), United States. Dept. of Energy (DE-FE0009738)

Infiltration of water in dry porous media is subject to a powerful gravity-driven instability. Although the phenomenon of unstable infiltration is well known, its description using continuum mathematical models has posed a significant challenge for several decades. The classical model of water flow in the unsaturated flow, the Richards equation, is unable to reproduce the instability. Here, we present a computational study of a model of unsaturated flow in porous media that extends the Richards equation and is capable of predicting the instability and captures the key features of gravity fingering quantitatively. The extended model is based on a phase-field formulation and is fourth-order in space. The new model poses a set of challenges for numerical discretizations, such as resolution of evolving interfaces, stiffness in space and time, treatment of singularly perturbed equations, and discretization of higher-order spatial partial–differential operators. We develop a numerical algorithm based on Isogeometric Analysis, a generalization of the finite element method that permits the use of globally-smooth basis functions, leading to a simple and efficient discretization of higher-order spatial operators in variational form. We illustrate the accuracy, efficiency and robustness of our method with several examples in two and three dimensions in both homogeneous and strongly heterogeneous media. We simulate, for the first time, unstable gravity-driven infiltration in three dimensions, and confirm that the new theory reproduces the fundamental features of water infiltration into a porous medium. Our results are consistent with classical experimental observations that demonstrate a transition from stable to unstable fronts depending on the infiltration flux., United States. Dept. of Energy (Early Career Award Grant DE-SC0003907)

We perform a stability and convergence analysis of sequential methods for coupled flow and geomechanics, in which the mechanics sub-problem is solved first. We consider slow deformations, so that inertia is negligible and the mechanical problem is governed by an elliptic equation. We use Biot’s self-consistent theory to obtain the classical parabolic-type flow problem. We use a generalized midpoint rule (parameter α between 0 and 1) time discretization, and consider two classical sequential methods: the drained and undrained splits. The von Neumann method provides sharp stability estimates for the linear poroelasticity problem. The drained split with backward Euler time discretization ( α = 1) is conditionally stable, and its stability depends only on the coupling strength, and it is independent of time step size. The drained split with the midpoint rule ( α = 0.5) is unconditionally unstable. The mixed time discretization, with α = 1.0 for mechanics and α = 0.5 for flow, has the same stability properties as the backward Euler scheme. The von Neumann method indicates that the undrained split is unconditionally stable when α ⩾ 0.5. We extend the stability analysis to the nonlinear regime (poro-elastoplasticity) via the energy method. It is well known that the drained split does not inherit the contractivity property of the continuum problem, thereby precluding unconditional stability. For the undrained split we show that it is B-stable (therefore unconditionally stable at the algorithmic level) when α ⩾ 0.5. We also analyze convergence of the drained and undrained splits, and derive the a priori error estimates from matrix algebra and spectral analysis. We show that the drained split with a fixed number of iterations is not convergent even when it is stable. The undrained split with a fixed number of iterations is convergent for a compressible system (i.e., finite Biot modulus). For a nearly-incompressible system (i.e., very large Biot modulus), the undrained split loses first-order accuracy, and becomes non-convergent in time. We also study the rate of convergence of both splits when they are used in a fully-iterated sequential scheme. When the medium permeability is high or the time step size is large, which corresponds to a high diffusion of pressure, the error amplification of the drained split is lower and therefore converges faster than the undrained split. The situation is reversed in the case of low permeability and small time step size. We provide numerical experiments supporting all the stability and convergence estimates of the drained and undrained splits, in the linear and nonlinear regimes. We also show that our spatial discretization (finite volumes for flow and finite elements for mechanics) removes the well-documented spurious instability in consolidation problems at early times.

We analyze stability and convergence of sequential implicit methods for coupled flow and geomechanics, in which the flow problem is solved first. We employ the von Neumann and energy methods for linear and nonlinear problems, respectively. We consider two sequential methods with the generalized midpoint rule for tn+α, where α is the parameter of time discretization: namely, the fixed-strain and fixed-stress splits. The von Neumann method indicates that the fixed-strain split is only conditionally stable, and that its stability limit is a coupling strength less than unity if α ⩾ 0.5. On the other hand, the fixed-stress split is unconditionally stable when α ⩾ 0.5, the amplification factors of the fixed-stress split are different from those of the undrained split and are identical to the fully coupled method. Unconditional stability of the fixed-stress split is also obtained from the energy method for poroelastoplasticity. We show that the fixed-stress split is contractive and B-stable when α ⩾ 0.5. We also estimate the convergence behaviors for the two sequential methods by the matrix based and spectral analyses for the backward Euler method in time. From the estimates, the fixed-strain split may not be convergent with a fixed number of iterations particularly around the stability limit even though it is stable. The fixed-stress split, however, is convergent for a fixed number of iterations, showing better accuracy than the undrained split. Even when we cannot obtain the exact local bulk modulus (or exact rock compressibility) at the flow step a priori due to complex boundary conditions or the nonlinearity of the materials, the fixed-stress split can still provide stability and convergence by an appropriate estimation of the local bulk modulus, such as the dimension-based estimation, by which the employed local bulk modulus is less stiff than the exact local bulk modulus. We provide numerical examples supporting all the estimates of stability and convergence for the fixed-strain and fixed-stress splits.

A promising way to mitigate global warming is to sequester CO[subscript 2] in deep saline aquifers. In order to determine which aquifers are the best for sequestration, it is helpful to estimate how much CO[subscript 2] they can store. Currently, this is difficult because both the pressure buildup from injection and the volume of available pore space have been identified as constraints, but have not been compared to determine which is more important. In this study, we evaluate their relative importance using simple, but dynamic models of how pressure rises during injection and how CO[subscript 2] becomes trapped in the pore space. Our results show that the more important constraint depends on the properties of the aquifer and how the CO[subscript 2] is injected.

We incorporate CO2 dissolution due to convective mixing into a sharp-interface mathematical model for the post-injection migration of a plume of CO2 in a saline aquifer. The model captures CO2 migration due to groundwater flow and aquifer slope, as well as residual trapping and dissolution. We also account for the tongued shape of the plume at the end of the injection period. We solve the model numerically and identify three regimes in CO2 migration with dissolution, based on how quickly the brine beneath the plume saturates with dissolved CO2. When the brine saturates slowly relative to plume migration, dissolution is controlled by the dimensionless dissolution rate. When the brine saturates “instantaneously” relative to plume migration, dissolution is instead controlled by the solubility of CO2 in brine. We show that dissolution can lead to a several-fold increase in storage efficiency. In a companion paper, we study migration and pressure limitations on storage capacity [Szulczewski et al., GHGT-10, Paper 917 (2010)].

Safely sequestering CO[subscript 2] in a deep saline aquifer requires calculating how much CO[subscript 2] the aquifer can store. Since offsetting nationwide emissions requires sequestering large quantities of CO[subscript 2], this calculation should apply at the large scale of geologic basins. The only method to calculate storage capacity at the basin scale, however, is not derived from multiphase flow dynamics, which play a critical role in CO[subscript 2] storage. In this study, we explain a new model to calculate basin-scale storage capacity that is derived from flow dynamics and captures the dynamic phenomena of gravity override and capillary trapping. Despite the fact that the model is dynamic, it is simple since it is a closed form expression with few terms. We demonstrate how to apply it on the Fox Hills Sandstone in the Powder River Basin.

We propose a novel approach to flash calculation, with particular application to negative flash. The basis of the method is a parameterization of the tie-line field. Rather than solving the Rachford–Rice equation (or any of its variants) we solve directly for the parameters defining the tie line. For an N-component system, our approach leads to a system of N − 2 quadratic equations, which we solve efficiently using a Newton method. The iterative method is very robust: unlike other negative flash procedures, the solution displays continuous dependence on the overall composition, even in the transition to negative concentrations. We illustrate the properties and behavior of the proposed approach on three-component and four-component systems, and we then generalize the method to systems of N components. From the global triangular structure of the system with constant K-values, it follows that the system of N−2 quadratic equations can only have two roots. For the important case of three components, the flash calculation is explicit.

In this paper we extend to three-phase flow the nonequilibrium formalism proposed by Barenblatt and co-workers for two-phase porous media flow. The underlying idea is to include nonequilibrium effects by introducing a pair of effective water and gas saturations, which are linked to the actual saturations by a local evolution equation. We illustrate and analyze how nonequilibrium effects lead to qualitative and quantitative differences in the solution of the three-phase flow equations.

Streamline methods for subsurface-flow simulation have received renewed attention as fast alternatives to traditional finite volume or finite element methods. Key aspects of streamline simulation are the accurate tracing of streamlines and the computation of travel time along individual streamlines. In this paper, we propose a new streamline tracing framework that enables the extension of streamline methods to unstructured grids composed of triangular or quadrilateral elements and populated with heterogeneous full-tensor coefficients. The proposed method is based on the mathematical framework of mixed finite element methods, which provides approximations of the velocity field that are especially suitable for streamline tracing. We identify and implement two classes of velocity spaces: the lowest-order Raviart-Thomas space (low-order tracing) and the Brezzi-Douglas-Marini space of order one (high-order tracing), both on triangles and quadrilaterals. We discuss the implementation of the streamline tracing method in detail, and we investigate the performance of the proposed tracing strategy by means of carefully designed test cases. We conclude that, for the same computational cost, high-order tracing is more accurate (smaller error in the time-of-flight) and robust (less sensitive to grid distortion) than low-order tracing.

Pore-scale physics, laboratory investigations, and field experience, dictate that three-phase relative permeabilities exhibit strong dependence on the saturation path and the saturation history. Such dependence is especially relevant in immiscible water-alternating-gas (WAG) processes, which are characterized by a sequence of three-phase drainage and imbibition cycles. In this paper, we study the influence of relative permeability hysteresis on the field-scale predictions of WAG injection. Because their measurement is difficult and time-consuming, three-phase relative permeabilities are usually interpolated from two-phase data. The errors associated with this procedure have been investigated by Oak (SPE 20183), who reported that interpolated values might differ significantly from experimental ones. The effect of using different interpolation models in field-scale simulations has been illustrated by a number of authors, who found that recovery predictions could be significantly different depending on the three-phase relative permeability model. Here, we study the impact of using history-dependent saturation functions in reservoir simulations. First, we investigate the degree of accuracy with which different hysteretic models reproduce Oak's three-phase relative permeability data. In doing so, we assess the validity of existing models, and we identify the model parameters subject to most uncertainty. Our analysis suggests that current models account for some of the hysteretic behavior observed experimentally, but do not reproduce experimental measurements adequately during cyclic water/gas injection. Second, we illustrate how the use of a hysteretic relative permeability model affects reservoir simulations. We use a synthetic model of a quarter five-spot pattern in a homogenous reservoir, and a more realistic heterogeneous reservoir modified from the PUNQ-S3 model. We find that there is striking disparity in the simulation results depending on whether a hysteretic or a nonhysteretic model is employed. Therefore, we conclude that (1) it is essential to incorporate hysteresis in the relative permeabilities in order to obtain accurate predictions of realistic immiscible WAG processes; and (2) enhancements are needed to improve the predictive capabilities of current relative permeability hysteresis models.

We present a variational multiscale formulation for the numerical solution of one-dimensional systems of conservation laws. The key idea of the proposed formulation, originally presented by Hughes [Comput. Methods Appl. Mech. Engrg., 127 (1995) 387–401], is a multiple-scale decomposition into resolved grid scales and unresolved subgrid scales. Incorporating the effect of the subgrid scales onto the coarse scale problem results in a finite element method with enhanced stability properties, capable of accurately representing the sharp features of the solution. In the formulation developed herein, the multiscale split is invoked prior to any linearization of the equations. Special attention is given to the choice of the matrix of stabilizing coefficients and the discontinuity-capturing diffusion. The methodology is applied to the one-dimensional simulation of three-phase flow in porous media, and the shallow water equations. These numerical simulations clearly show the potential and applicability of the formulation for solving highly nonlinear, nearly hyperbolic systems on very coarse grids. Application of the numerical formulation to multidimensional problems is presented in a forthcoming paper.

We present a stabilized finite element method for the numerical solution of multiphase flow in porous media, based on a multiscale decomposition of pressures and fluid saturations into resolved (or grid) scales and unresolved (or subgrid) scales. The multiscale split is invoked in a variational setting, which leads to a rigorous definition of a grid scale problem and a subgrid scale problem. The subgrid problem is modeled using an algebraic approximation. This model requires the definition of a matrix of intrinsic time scales, which we design based on stability considerations. We illustrate the performance of the method with simulations of a waterflood in a heterogeneous oil reservoir. The proposed method yields stable, highly accurate solutions on very coarse grids, which we compare with those obtained by the classical Galerkin method or the upstream finite difference method. Although this paper is restricted to multiphase flow in porous media, the formulation is quite general and can be applied to other nonlinear systems of conservation laws, like the shallow water equations.

Groundwater inflows into tunnels constructed in fractured bedrocks not only constitute an important factor controlling the rate of advancement in driving the tunnel but may pose potential hazards. Drawdowns caused by tunnel construction may also induce geotechnical and environmental impacts. Here we present a numerical methodology for the dynamic simulation of the hydrogeological transient conditions induced by the tunnel front advance. The methodology is based on the use of a Cauchy boundary condition at the points lying along the tunnel according to which water discharge, Q , is computed as the product of a leakage coefficient, α , and the head difference, ( H − h ), where H is the prescribed head at the tunnel wall and h is the hydraulic head in the fractured rock in the close vicinity of the tunnel. At a given position of the tunnel, α is zero until the tunnel reaches such position when it is assigned a positive value. The use of step-wise time functions for α allows an efficient and accurate simulation of the transient hydrogeological conditions at and around the tunnel during the excavation process. The methodology has been implemented in TRANMEF-3, a finite element computer code for groundwater flow in 3D fractured media developed at the University of A Coruna, Spain, and has been used to simulate the impact of a tunnel on the groundwater system at the Aspo island (Sweden). This tunnel was constructed to access an underground laboratory for research on radioactive waste disposal. The large amount of available data at this site provides a unique opportunity to test the performance of the numerical model and the proposed methodology for tunnel advance. With just minor calibration, the numerical model is able to reproduce accurately the measurements of inflows into the tunnel at several reaches and hydraulic heads at surface-drilled boreholes. These results obtained at the Aspo site lead us to conclude that accurate predictions of the transient hydrogeological responses induced by tunneling works in fractured bedrocks, can be achieved provided that a sound hydrogeological characterization of large-scale fracture zones is available.