Showing papers in "Computational Geosciences in 2015"

Reactive transport codes for subsurface environmental simulation

Abstract: A general description of the mathematical and numerical formulations used in modern numerical reactive transport codes relevant for subsurface environmental simulations is presented. The formulations are followed by short descriptions of commonly used and available subsurface simulators that consider continuum representations of flow, transport, and reactions in porous media. These formulations are applicable to most of the subsurface environmental benchmark problems included in this special issue. The list of codes described briefly here includes PHREEQC, HPx, PHT3D, OpenGeoSys (OGS), HYTEC, ORCHESTRA, TOUGHREACT, eSTOMP, HYDROGEOCHEM, CrunchFlow, MIN3P, and PFLOTRAN. The descriptions include a high-level list of capabilities for each of the codes, along with a selective list of applications that highlight their capabilities and historical development.

Phase-field modeling of a fluid-driven fracture in a poroelastic medium

Abstract: In this paper, we present a phase field model for a fluid-driven fracture in a poroelastic medium. In our previous work, the pressure was assumed given. Here, we consider a fully coupled system where the pressure field is determined simultaneously with the displacement and the phase field. To the best of our knowledge, such a model is new in the literature. The mathematical model consists of a linear elasticity system with fading elastic moduli as the crack grows, which is coupled with an elliptic variational inequality for the phase field variable and with the pressure equation containing the phase field variable in its coefficients. The convex constraint of the variational inequality assures the irreversibility and entropy compatibility of the crack formation. The phase field variational inequality contains quadratic pressure and strain terms, with coefficients depending on the phase field unknown. We establish existence of a solution to the incremental problem through convergence of a finite dimensional approximation. Furthermore, we construct the corresponding Lyapunov functional that is linked to the free energy. Computational results are provided that demonstrate the effectiveness of this approach in treating fluid-driven fracture propagation.

Dimensionally reduced flow models in fractured porous media: crossings and boundaries

Abstract: For the simulation of fractured porous media, a common approach is the use of co-dimension one models for the fracture description. In order to simulate correctly the behavior at fracture crossings, standard models are not sufficient because they either cannot capture all important flow processes or are computationally inefficient. We propose a new concept to simulate co-dimension one fracture crossings and show its necessity and accuracy by means of an example and a comparison to a literature benchmark. From the application point of view, often the pressure is known only at a limited number of discrete points and an interpolation is used to define the boundary condition at the remaining parts of the boundary. The quality of the interpolation, especially in fracture models, influences the global solution significantly. We propose a new method to interpolate boundary conditions on boundaries that are intersected by fractures and show the advantages over standard interpolation methods.

Particle tracking approach for transport in three-dimensional discrete fracture networks

Abstract: The discrete fracture network (DFN) model is a method to mimic discrete pathways for fluid flow through a fractured low-permeable rock mass, and may be combined with particle tracking simulations to address solute transport. However, experience has shown that it is challenging to obtain accurate transport results in three-dimensional DFNs because of the high computational burden and difficulty in constructing a high-quality unstructured computational mesh on simulated fractures. We present a new particle tracking capability, which is adapted to control volume (Voronoi polygons) flow solutions on unstructured grids (Delaunay triangulations) on three-dimensional DFNs. The locally mass-conserving finite-volume approach eliminates mass balance-related problems during particle tracking. The scalar fluxes calculated for each control volume face by the flow solver are used to reconstruct a Darcy velocity at each control volume centroid. The groundwater velocities can then be continuously interpolated to any point in the domain of interest. The control volumes at fracture intersections are split into four pieces, and the velocity is reconstructed independently on each piece, which results in multiple groundwater velocities at the intersection, one for each fracture on each side of the intersection line. This technique enables detailed particle transport representation through a complex DFN structure. Verified for small DFNs, the new simulation capability enables numerical experiments on advective transport in large DFNs to be performed. We demonstrate this particle transport approach on a DFN model using parameters similar to those of crystalline rock at a proposed geologic repository for spent nuclear fuel in Forsmark, Sweden.

Iterative regularization for ensemble data assimilation in reservoir models

Abstract: We propose the application of iterative regularization for the development of ensemble methods for solving Bayesian inverse problems. In concrete, we construct (i) a variational iterative regularizing ensemble Levenberg-Marquardt method (IR-enLM) and (ii) a derivative-free iterative ensemble Kalman smoother (IR-ES). The aim of these methods is to provide a robust ensemble approximation of the Bayesian posterior. The proposed methods are based on fundamental ideas from iterative regularization methods that have been widely used for the solution of deterministic inverse problems (Katltenbacher et al. de Gruyter, Berlin 2008). In this work, we are interested in the application of the proposed ensemble methods for the solution of Bayesian inverse problems that arise in reservoir modeling applications. The proposed ensemble methods use key aspects of the regularizing Levenberg-Marquardt scheme developed by Hanke (Inverse Problems 13,79–95 1997) and that we recently applied for history matching in Iglesias (Comput. Geosci. 1–21 2013). Unlike most existing methods where the stopping criteria and regularization parameters are typically selected heuristically, in the proposed ensemble methods, the discrepancy principle is applied for (i) the selection of the regularization parameters and (ii) the early termination of the scheme. The discrepancy principle is key for the theory of iterative regularization, and the purpose of the present work is to apply this principle for the development of ensemble methods defined as iterative updates of solutions to linear ill-posed inverse problems. The regularizing and convergence properties of iterative regularization methods for deterministic inverse problems have long been established. However, the approximation properties of the proposed ensemble methods in the context of Bayesian inverse problems is an open problem. In the case where the forward operator is linear and the prior is Gaussian, we show that the tunable parameters of the proposed IR-enLM and IR-ES can be chosen so that the resulting schemes coincide with the standard randomized maximum likelihood (RML) and the ensemble smoother (ES), respectively. Therefore, the proposed methods sample from the posterior in the linear-Gaussian case. Similar to RML and ES methods, in the nonlinear case, one may not conclude that the proposed methods produce samples from the posterior. The present work provides a numerical investigation of the performance of the proposed ensemble methods at capturing the posterior. In particular, we aim at understanding the role of the tunable parameters that arise from the application of iterative regularization techniques. The numerical framework for our investigations consists of using a state-of-the art Markov chain Monte Carlo (MCMC) method for resolving the Bayesian posterior from synthetic experiments. The resolved posterior via MCMC then provides a gold standard against to which compare the proposed IR-enLM and IR-ES. Our numerical experiments show clear indication that the regularizing properties of the regularization methods applied for the computation of each ensemble have significant impact of the approximation properties of the proposed ensemble methods at capturing the Bayesian posterior. Furthermore, we provide a comparison of the proposed regularizing methods with respect to some unregularized methods that have been typically used in the literature. Our numerical experiments showcase the advantage of using iterative regularization for obtaining more robust and stable approximation of the posterior than unregularized methods.

Data assimilation and uncertainty assessment for complex geological models using a new PCA-based parameterization

Abstract: In this paper, a recently developed parameterization procedure based on principal component analysis (PCA), which is referred to as optimization-based PCA (O-PCA), is generalized for use with a wide range of geological systems. In O-PCA, the mapping between the geological model in the full-order space and the low-dimensional subspace is framed as an optimization problem. The O-PCA optimization involves the use of regularization and bound constraints, which act to extend substantially the ability of PCA to model complex (non-Gaussian) systems. The basis matrix required by O-PCA is formed using a set of prior realizations generated by a geostatistical modeling package. We show that, by varying the form of the O-PCA regularization terms, different types of geological scenarios can be represented. Specific systems considered include binary-facies, three-facies and bimodal channelized models, and bimodal deltaic fan models. The O-PCA parameterization can be applied to generate random realizations, though our focus here is on its use for data assimilation. For this application, O-PCA is combined with the randomized maximum likelihood (RML) method to provide a subspace RML procedure that can be applied to non-Gaussian models. This approach provides multiple history-matched models, which enables an estimate of prediction uncertainty. A gradient procedure based on adjoints is used for the minimization required by the subspace RML method. The gradient of the O-PCA mapping is determined analytically or semi-analytically, depending on the form of the regularization terms. Results for two-dimensional oil-water systems, for several different geological scenarios, demonstrate that the use of O-PCA and RML enables the generation of posterior reservoir models that honor hard data, retain the large-scale connectivity features of the geological system, match historical production data, and provide an estimate of prediction uncertainty. MATLAB code for the O-PCA procedure, along with examples for three-facies and bimodal models, is included as online Supplementary Material.

Non-isothermal, Multi-phase, Multi-component Flows through Deformable Methane Hydrate Reservoirs.

Abstract: We present a hydro-geomechanical model for subsurface methane hydrate systems. Our model considers kinetic hydrate phase change and non-isothermal, multi-phase, multi-component flow in elastically deforming soils. The model accounts for the effects of hydrate phase change and pore pressure changes on the mechanical properties of the soil. It also accounts for the effect of soil deformation on the fluid-solid interaction properties relevant to reaction and transport processes (e.g., permeability, capillary pressure, and reaction surface area). We discuss a ’cause-effect’ based decoupling strategy for the model and present our numerical discretization and solution scheme. We then proceed to identify the important model components and couplings which are most vital for a hydro-geomechanical hydrate simulator, namely, (1) dissociation kinetics, (2) hydrate phase change coupled with non-isothermal two phase two component flow, (3) two phase flow coupled with linear elasticity (poroelasticity coupling), and finally (4) hydrate phase change coupled with poroelasticity (kinetics-poroelasticity coupling). To show the versatility of our hydrate model, we numerically simulate test problems where, for each problem, we methodically isolate one out of the four aforementioned model components or couplings. A special emphasis is laid on the kinetics-poroelasticity coupling for which we present a test problem where an axially loaded hydrate bearing sand sample experiences a spontaneous shift in the hydrate stability curve causing the hydrate to melt. For this problem, we present an analytical solution for pore-pressure, which we subsequently use to test the accuracy of the numerical scheme. Finally, we present a more complex 3D example where all the major model components are put together to give an idea of the model capabilities. The setting is based on a subsurface hydrate reservoir which is destabilized through depressurization using a low pressure gas well. In this example, we simulate the melting of hydrate, methane gas generation, and the resulting ground subsidence and stress build-up in the vicinity of the well.

Implementation and evaluation of permeability-porosity and tortuosity-porosity relationships linked to mineral dissolution-precipitation

Abstract: Changes of porosity, permeability, and tortuosity due to physical and geochemical processes are of vital importance for a variety of hydrogeological systems, including passive treatment facilities for contaminated groundwater, engineered barrier systems (EBS), and host rocks for high-level nuclear waste (HLW) repositories. Due to the nonlinear nature and chemical complexity of the problem, in most cases, it is impossible to verify reactive transport codes analytically, and code intercomparisons are the most suitable method to assess code capabilities and model performance. This paper summarizes model intercomparisons for six hypothetical scenarios with generally increasing geochemical or physical complexity using the reactive transport codes CrunchFlow, HP1, MIN3P, PFlotran, and TOUGHREACT. Benchmark problems include the enhancement of porosity and permeability through mineral dissolution, as well as near complete clogging due to localized mineral precipitation, leading to reduction of permeability and tortuosity. Processes considered in the benchmark simulations are advective-dispersive transport in saturated media, kinetically controlled mineral dissolution-precipitation, and aqueous complexation. Porosity changes are induced by mineral dissolution-precipitation reactions, and the Carman-Kozeny relationship is used to describe changes in permeability as a function of porosity. Archie’s law is used to update the tortuosity and the pore diffusion coefficient as a function of porosity. Results demonstrate that, generally, good agreement is reached amongst the computer models despite significant differences in model formulations. Some differences are observed, in particular for the more complex scenarios involving clogging; however, these differences do not affect the interpretation of system behavior and evolution.

Comprehensive framework for gradient-based optimization in closed-loop reservoir management

Abstract: An efficient, robust, and flexible adjoint-based computational framework for performing closed-loop reservoir management is developed and applied. The methodology includes gradient-based production optimization and data assimilation (history matching). Flexibility is achieved through the use of automatic differentiation (AD) within the reservoir simulation, production optimization, and history matching modules. The use of AD will also facilitate the application of closed-loop reservoir management to physical models of higher complexity. A fast sequential convex programming (SCP) solver based on the method of moving asymptotes (MMA) is applied for the production optimization component of the closed-loop. This technique is shown to outperform the sequential quadratic programming (SQP) method, which is commonly used for production optimization computations. The history matching component of the workflow integrates both production data and proxy seismic measurements into a unified adjoint-based data assimilation framework. The effect of noisy data, and data of different types, on the accuracy of the history matching component is assessed. The overall closed-loop reservoir management methodology is tested using the well-documented Brugge model. Results demonstrate the efficient performance of the individual closed-loop components and the improvement in net present value that is achieved using these procedures.

Multirate time integration for coupled saturated/unsaturated porous medium and free flow systems

Abstract: A model formulation to describe fluid flows in coupled saturated/unsaturated porous medium and adjacent free flow regions is proposed. The Stokes equations are applied in the free flow domain, while the Richards equation is used to model the porous medium system. These two flow problems are coupled at the fluid-porous interface via an appropriate set of interface conditions. A multiple-time-step scheme is developed to solve the coupled problem efficiently. Numerical simulation results are presented for a model problem and a realistic setting that demonstrate the convergence and efficiency of the proposed computational algorithm. Time-splitting multistep methods can be successfully applied for modeling other physical systems where the processes evolve on different time scales, and these potential extensions are discussed.

Benchmark reactive transport simulations of a column experiment in compacted bentonite with multispecies diffusion and explicit treatment of electrostatic effects

Abstract: Bentonite clay is considered as a potential buffer and backfill material in subsurface repositories for high-level nuclear waste. As a result of its low permeability, transport of water and solutes in compacted bentonite is driven primarily by diffusion. Developing models for species transport in bentonite is complicated, because of the interaction of charged species and the negative surface charge of clay mineral surfaces. The effective diffusion coefficient of an ion in bentonite depends on the ion’s polarity and valence, on the ionic strength of the solution, and on the bulk dry density of the bentonite. These dependencies need to be understood and incorporated into models if one wants to predict the effectiveness of bentonite as a barrier to radionuclides in a nuclear repository. In this work, we present a benchmark problem for reactive transport simulators based on a flow-through experiment carried out on a saturated bentonite core. The measured effluent composition shows the complex interplay of species transport in a charged medium in combination with sorption and mineral precipitation/dissolution reactions. The codes compared in this study are PHREEQC, CrunchFlow, FLOTRAN, and MIN3P. The benchmark problem is divided into four component problems of increasing complexity, leading up to the main problem which addresses the effects of advective and diffusive transport of ions through bentonite with explicit treatment of electrostatic effects. All codes show excellent agreement between results provided that the activity model, Debye-Huckel parameters, and thermodynamic data used in the simulations are consistent. A comparison of results using species-specific diffusion and uniform species diffusion reveals that simulated species concentrations in the effluent differ by less than 8 %, and that these differences vanish as the system approaches steady state.

Benchmarks for multicomponent diffusion and electrochemical migration

Abstract: In multicomponent electrolyte solutions, the tendency of ions to diffuse at different rates results in a charge imbalance that is counteracted by the electrostatic coupling between charged species leading to a process called “electrochemical migration” or “electromigration.” Although not commonly considered in solute transport problems, electromigration can strongly affect mass transport processes. The number of reactive transport models that consider electromigration has been growing in recent years, but a direct model intercomparison that specifically focuses on the role of electromigration has not been published to date. This contribution provides a set of three benchmark problems that demonstrate the effect of electric coupling during multicomponent diffusion and electrochemical migration and at the same time facilitate the intercomparison of solutions from existing reactive transport codes. The first benchmark focuses on the 1D transient diffusion of HNO3 (pH = 4) in a NaCl solution into a fixed concentration reservoir, also containing NaCl—but with lower HNO3 concentrations (pH = 6). The second benchmark describes the 1D steady-state migration of the sodium isotope 22Na triggered by sodium chloride diffusion in neutral pH water. The third benchmark presents a flow-through problem in which transverse dispersion is significantly affected by electromigration. The system is described by 1D transient and 2D steady-state models. Very good agreement on all of the benchmarks was obtained with the three reactive transport codes used: CrunchFlow, MIN3P, and PHREEQC.

ParCrunchFlow: an efficient, parallel reactive transport simulation tool for physically and chemically heterogeneous saturated subsurface environments

Abstract: Understanding the interactions between physical, geochemical, and biological processes in the shallow subsurface is integral to the development of effective contamination remediation techniques, or the accurate quantification of nutrient fluxes and biogeochemical cycling. Hydrology is a primary control on the behavior of shallow subsurface environments and must be realistically represented if we hope to accurately model these systems. ParCrunchFlow is a new parallel reactive transport model that was created by coupling a multicomponent geochemical code (CrunchFlow) with a parallel hydrologic model (ParFlow). These models are coupled in an explicit operator-splitting manner. ParCrunchFlow can simulate three-dimensional multicomponent reactive transport in highly resolved, field-scale systems by taking advantage of ParFlow’s efficient parallelism and robust hydrologic abilities, and CrunchFlow’s extensive geochemical abilities. Here, the development of ParCrunchFlow is described and two simple verification simulations are presented. The parallel performance is evaluated and shows that ParCrunchFlow has the ability to simulate very large problems. A series of simulations involving the biologically mediated reduction of nitrate in a floodplain aquifer were conducted. These floodplain simulations show that this code enables us to represent more realistically the variability in chemical concentrations observed in many field-scale systems. The numerical formulation implemented in ParCrunchFlow minimizes numerical dispersion and allows the use of higher-order explicit advection schemes. The effects that numerical dispersion can have on finely resolved, field-scale reactive transport simulations have been evaluated. The smooth gradients produced by a first-order advection scheme create an artificial mixing effect, which decreases the spatial variance in solute concentrations and leads to an increase in overall reaction rates. The work presented here is the first step in a larger effort to couple these models in a transient, variably saturated surface-subsurface framework, with additional geochemical abilities.

Benchmarks for multicomponent reactive transport across a cement/clay interface

Abstract: The use of the subsurface for CO2 storage, geothermal energy generation, and nuclear waste disposal will greatly increase the interaction between clay(stone) and concrete. The development of models describing the mineralogical transformations at this interface is complicated, because contrasting geochemical conditions (Eh, pH, solution composition, etc.) induce steep concentration gradients and a high mineral reactivity. Due to the complexity of the problem, analytical solutions are not available to verify code accuracy, rendering code intercomparisons as the most efficient method for assessing code capabilities and for building confidence in the used model. A benchmark problem was established for tackling this issue. We summarize three scenarios with increasing geochemical complexity in this paper. The processes considered in the simulations are diffusion-controlled transport in saturated media under isothermal conditions, cation exchange reactions, and both local equilibrium and kinetically controlled mineral dissolution-precipitation reactions. No update of the pore diffusion coefficient as a function of porosity changes was considered. Seven international teams participated in this benchmarking exercise. The reactive transport codes used (TOUGHREACT, PHREEQC, with two different ways of handling transport, CRUNCH, HYTEC, ORCHESTRA, MIN3P-THCm) gave very similar patterns in terms of predicted solute concentrations and mineral distributions. Some differences linked to the considered activity models were observed, but they do not bias the general system evolution. The benchmarking exercise thus demonstrates that a reactive transport modelling specification for long-term performance assessment can be consistently addressed by multiple simulators.

Multi-objective optimization for rapid and robust optimal oilfield development under geological uncertainty

Abstract: Optimal well placement is critical to oil and gas field development. Typical workflows involve procedures to place a new well or a group of new wells in a reservoir in order to maximize some pre-defined reservoir performance metric. However, there are two main drawbacks with these traditional optimization approaches: First, the impact of geological uncertainty is often neglected or there may be no framework to include geological uncertainty. Second, traditional optimization techniques normally cannot meet the requirement of optimizing two or more conflicting objectives simultaneously—this may be useful when maximizing oil recovery while also minimizing water production. Consequently, in recent years, multiple objective optimization to obtain robust solutions that minimize the decision risk under geological uncertainty becomes a topic of renewed interest. Therefore, in this work, we develop a new work flow for well placement optimization while considering geological uncertainty in reservoir models. In general, when considering geological uncertainty, the primary goal is to maximize the mean net present value (NPV) over all realizations. However, restricting the search to simply maximizing the mean NPV may be inappropriate or inadequate for decision-making. A more reasonable choice is to maximize the mean NPV while minimizing the spread of the optimal NPV’s obtained for each realization. Therefore, in this work, we apply multi-objective optimization techniques to maximize the mean and minimize the variance of NPV values over all geological realizations to provide robust well placement solutions for decision-makers to select according to their risk attitude towards field development plans.

Reactive transport benchmarks for subsurface environmental simulation

Abstract: Over the last 20 years, we have seen firsthand the evolution of multicomponent reactive transport modeling and the expanding range and increasing complexity of subsurface applications it is being used to address. There is a growing reliance on reactive transport modeling (RTM) to address some of the most compelling issues facing our planet: climate change, nuclear waste management, contaminant remediation, and pollution prevention. While these issues are motivating the development of new and improved capabilities for subsurface environmental modeling using RTM (e.g., biogeochemistry from cell-scale physiology to continental-scale terrestrial ecosystems, nonisothermal multiphase conditions, coupled geomechanics), there remain longstanding challenges in characterizing the natural variability of hydrological, biological, and geochemical properties in subsurface environments and limited success in transferring models between sites and across scales. An equally important trend over the last 20 years is the evolution of modeling from a service sought out after data has been collected to a multifaceted research approach that provides (1) an organizing principle for characterization and monitoring activities; (2) a systematic framework for identifying knowledge gaps, developing and integrating new knowledge; and (3) a mechanistic understanding that represents the collective wisdom of the participating scientists and engineers. There are now large multidisciplinary projects wheremore » the research approach is model-driven, and the principal product is a holistic predictive simulation capability that can be used as a test bed for alternative conceptualizations of processes, properties, and conditions. Much of the future growth and expanded role for RTM will depend on its continued ability to exploit technological advancements in the earth and environmental sciences. Advances in measurement technology, particularly in molecular biology (genomics), isotope fractionation, and high-resolution X-ray spectroscopy, have created new lines of research that can be used to inform the conceptualization of reactions and rate laws and validate mechanistic models. For example, spectroscopy has identified the oxidation states of key components and elemental distributions at increasingly smaller scales and lower concentrations; molecular biology has progressed from identifying the presence of microbes to characterization of which microbial communities are active and what they are doing (i.e., microbial function), which has led in turn to the identification of active processes under conditions beyond what analytical chemistry can discern; isotope ratios in pore water and solid phases that can be used to distinguish between biotic from abiotic processes, sorption from precipitation, and origin and age of groundwater. The other noteworthy development that is expanding the role of RTM in subsurface environmental modeling is he advance in computational technology that is enabling the simulation of more coupled processes with increasing mechanistic detail. In some cases, this involves the inclusion of more reactive species and/or microbial populations in the simulations; in other cases, the impact is through the ability to achieve high resolution of property distributions over longer simulated times. To achieve these ambitious objectives for subsurface reactive transport simulation, the subsurface science and engineering community is being driven to provide accurate assessments of engineering performance and risk for important issues with far-reaching consequences. As a result, the complexity and detail of subsurface processes, properties, and conditions that can be simulated have significantly expanded. This expansion was enabled, in part, by advances in measurement technology, computing technology, and numerical techniques.« less

Investigation of the effective permeability of vuggy or fractured porous media from a Darcy-Brinkman approach

Abstract: In this paper, the macroscopic representation of one-phase incompressible flow in fractured and cavity (or vuggy) porous media is studied from theoretical and numerical points of view. A single-domain (or equivalently a Darcy-Brinkman) type of approach is followed to describe the momentum transport at Darcy scale where the fracture or cavity region and porous matrix region are well identified. The Darcy scale model is upscaled yielding a macroscopic momentum equation operating on the equivalent homogeneous medium. Numerical solution to the associated closure problem is proposed in order to compute the effective permeability. Numerical results on some model fractured and cavity media are discussed and compared to some analytical results.

Iterative Bayesian inversion with Gaussian mixtures: finite sample implementation and large sample asymptotics

Abstract: Approximate solutions for Bayesian estimation in large scale models is a topic under investigation in many scientific communities. We define an iterative method based on the adaptive Gaussian mixture filter with batch updates as a robust alternative to adaptive importance sampling. We prove asymptotic optimality under certain conditions, contrary to other methods discussed where the sample distribution depends on the nonlinearity and scaling of the model. The finite sample implementation of the method is compared to an ensemble smoother with multiple data assimilation and an ensemble-based randomized maximum likelihood approach on a synthetic 1D reservoir model.

A reactive transport benchmark on heavy metal cycling in lake sediments

Abstract: Sediments are active recipients of anthropogenic inputs, including heavy metals, but may be difficult to interpret without the use of numerical models that capture sediment-metal interactions and provide an accurate representation of the intricately coupled sedimentological, geochemical, and biological processes. The focus of this study is to present a benchmark problem on heavy metal cycling in lake sediments and to compare reactive transport models (RTMs) in their treatment of the local-scale physical and biogeochemical processes. This benchmark problem has been developed based on a previously published reactive-diffusive model of metal transport in the sediments of Lake Coeur d’Alene, Idaho. Key processes included in this model are microbial reductive dissolution of iron hydroxides (i.e., ferrihydrite), the release of sorbed metals into pore water, reaction of these metals with biogenic sulfide to form sulfide minerals, and sedimentation driving the burial of ferrihydrite and other minerals. This benchmark thus considers a multicomponent biotic reaction network with multiple terminal electron acceptors (TEAs), Fickian diffusive transport, kinetic and equilibrium mineral precipitation and dissolution, aqueous and surface complexation, as well as (optionally) sedimentation. To test the accuracy of the reactive transport problem solution, four RTMs—TOUGHREACT (TR), CrunchFlow (CF), PHREEQC, and PHT3D—have been used. Without sedimentation, all four models are able to predict similar trends of TEAs and dissolved metal concentrations, as well as mineral abundances. TR and CF are further used to compare sedimentation and compaction test cases. Results with different sedimentation rates are captured by both models, but since the codes do not use the same formulation for compaction, the results differ for this test case. Although, both TR and CF adequately capture the trends of aqueous concentrations and mineral abundances, the difference in results highlights the need to consider further the conceptual and numerical models that link transport, biogeochemical reactions, and sedimentation.

Robust optimization of subsurface flow using polynomial chaos and response surface surrogates

Abstract: This study employs an inclusive framework for surrogate model-based optimization in the presence of parametric and spatial uncertainties. The framework is applied to optimize water injection rate for optimal hydrocarbon recovery from a synthetic subsurface model with uncertainty in the geological and fluid relative permeability properties. In one model of parametric uncertainty, geological properties such as the channel’s absolute permeability and the fault transmissibility multiplier and the fluid relative permeability parameters such as the residual oil saturation to water and the water relative permeability at residual oil are assumed to be non-informative. In another model, the channels positions are assumed uncertain and various realizations of the channelized permeability are parameterized and the spatial uncertainty is accounted for in the optimization. The uncertainty is quantified in each evaluation of the objective function via polynomial chaos expansions. The coefficients of polynomial chaos expansion are solved by probabilistic collocation method. The objective function is assigned with a risk-averse net present value computed from a distribution of values obtained from the probabilistic proxies. The proxies are updated for each round of objective function evaluation. Monte-Carlo simulations are also conducted to verify accuracy and to demonstrate the computational efficiency of the probabilistic collocation approach. The optimization is conducted in various random input cases (depending on the number of uncertain parameters) and for each case net present value is successfully maximized and optimal solutions of the water injection rates are determined.

Benchmarking the simulation of Cr isotope fractionation

Abstract: A benchmark problem set consisting of four problem levels was developed for the simulation of Cr isotope fractionation in 1D and 2D domains The benchmark is based on a recent field study where Cr(VI) reduction and accompanying Cr isotope fractionation occurs abiotically by an aqueous reaction with dissolved Fe 2+ (Wanner et al, 2012, Appl Geochem, 27, 644–662) The problem set includes simulation of the major processes affecting the Cr isotopic composition such as the dissolution of various Cr(VI) bearing minerals, fractionation during abiotic aqueous Cr(VI) reduction, and non-fractionating precipitation of Cr(III) as sparingly soluble Cr-hydroxide Accuracy of the presented solutions was ensured by running the problems with four well-established reactive transport modeling codes: TOUGHREACT, MIN3P, CRUNCHFLOW, and FLOTRAN Results were also compared with an analytical Rayleigh-type fractionation model An additional constraint on the correctness of the results was obtained by comparing output from the problem levels simulating Cr isotope fractionation with the corresponding ones only simulating bulk concentrations For all problem levels, model to model comparisons showed excellent agreement, suggesting that for the tested geochemical processes any code is capable of accurately simulating the fate of individual Cr isotopes

Bridging multipoint statistics and truncated Gaussian fields for improved estimation of channelized reservoirs with ensemble methods

Abstract: In this paper, we present a new parameterization of channelized reservoirs with two facies types, which is coupled with the ensemble Kalman filter (EnKF) and the iterative adaptive Gaussian mixture filter (IAGM) for history matching (HM) of production data. The main objectives are to match the past data within the model and measurement uncertainties and to preserve the geological realism in order to predict the future behavior of the reservoir. The parameterization bridges the method of Gaussian truncation with multipoint statistics from a training image. To generate an ensemble of channelized reservoirs, a multi-point geostatistical tool (SNESIM) is used in combination with a training image. The parameterization is performed in a Gaussian space by drawing from a conditional Gaussian distribution with truncation rules estimated from the ensemble, ensuring that the updates are always facies realizations. The EnKF is a HM method that updates the ensemble based only on the first two statistical moments, which is not enough to characterize the posterior when channelized structures are present. Therefore, we propose using the iterative version of AGM in combination with SNESIM in the resampling step to better handle nonlinearities and preserve the channelized structure of the ensemble members. The results presented show that the IAGM procedure is able to reduce the uncertainty in the updated ensemble with realistic geological structure, in addition to having good predictability and preservation of channels. We performed a comparison with IAGM applied directly to the permeability field.

GPU-accelerated 3D reconstruction of porous media using multiple-point statistics

Abstract: It is very important for the study of predicting fluid transport properties or mechanisms of fluid flow in porous media that the characteristics of porous media can be extracted in relatively smaller scales and then are copied in a larger or even arbitrary region to reconstruct virtual 3D porous media that have similar structures with the real porous media One of multiple-point statistics (MPS) method, the single normal equation simulation algorithm (SNESIM), has been widely used in reconstructing 3D porous media recently However, owing to its large CPU cost and rigid memory demand, the application of SNESIM has been limited in some cases To overcome this disadvantage, parallelization of SNESIM is performed on the compute unified device architecture (CUDA) kernels in the graphic processing unit (GPU) to reconstruct each node on simulation grids, combined with choosing the optimal size of data templates based on the entropy calculation towards the training image (TI) to acquire high-quality reconstruction with a low CPU cost; meanwhile, the integration of hard data and soft data is also included in the processing of CUDA kernels to improve the accuracy Representative elementary volumes (REVs) for porosity, variogram, and entropy are analyzed to guarantee that the scale of observation is large enough and parameters of concern are constant This parallel GPU-version 3D porous media reconstruction only requires relatively small size memory and benefits from the tremendous calculating power given by CUDA kernels to shorten the CPU time, showing its high efficiency for the reconstruction of porous media

High order discontinuous Galerkin method for simulating miscible flooding in porous media

Abstract: We present a high-order method for miscible displacement simulation in porous media. The method is based on discontinuous Galerkin discretization with weighted average stabilization technique and flux reconstruction post processing. The mathematical model is decoupled and solved sequentially. We apply domain decomposition and algebraic multigrid preconditioner for the linear system resulting from the high-order discretization. The accuracy and robustness of the method are demonstrated in the convergence study with analytical solutions and heterogeneous porous media, respectively. We also investigate the effect of grid orientation and anisotropic permeability using high-order discontinuous Galerkin method in contrast with cell-centered finite volume method. The study of the parallel implementation shows the scalability and efficiency of the method on parallel architecture. We also verify the simulation result on highly heterogeneous permeability field from the SPE10 model.

High-resolution finite element methods for 3D simulation of compositionally triggered instabilities in porous media

Abstract: The formation and development of patterns in the unstable interface between an injected fluid and hydrocarbons or saline aqueous phase in a porous medium can be driven by viscous effects and gravity. Numerical simulation of the so-called fingering is a challenge, which requires rigorous representation of the fluid flow and thermodynamics as well as highresolution discretization in order to minimize numerical artifacts. To achieve such a high resolution, we present higherorder 3D finite element methods for the simulation of fully compositional, three-phase and multi-component flow. This is based on a combination of the mixed hybrid finite element (MHFE) method for total fluid velocity and discontinuous Galerkin (DG) method for the species transport. The phase behavior is described by cubic or cubic-plus-association (CPA) equations of state. We present challenging numerical examples of compositionally triggered fingering at both the core and the large scale. Four additional test cases illustrate the robustness and efficiency of the proposed methods, which demonstrate their power for problems of this complexity. Results reveal three orders of magnitude improvement in CPU time in our method compared with the lowest-order finite difference method for some of the examples. Comparison between 3D and 2D results highlights the significance of dimensionality in the flow simulation.

Large deformation failure analysis of the soil slope based on the material point method

Abstract: The failure analysis of the soil slope is a very important topic in the field of geomechanics. Being a fully Lagrangian particle method, the material point method (MPM) has distinct advantages in solving the extremely large deformation problem. For both cohesive and non-cohesive soil slopes, the large deformation failure problems are analyzed using MPM and the Drucker-Prager constitutive model. For verification of the numerical method, the comparison between MPM and analytical solutions of the dam break problem is presented. Moreover, the numerical results by MPM are compared with the experimental results for the collapse of the aluminum-bar assemblage. Simulations reveal the cohesive soil slope under gravity has a shear band failure mode. Computational results show the reposed angle of non-cohesive soil slope is less than the internal friction angle, and the reason for this phenomenon is presented. The purpose of this study is to give a further understanding of the slope failure in different soil types and provide a computational tool for the failure analysis of soil slopes.

High-order finite difference modeling of tsunami generation in a compressible ocean from offshore earthquakes

Abstract: To study the full seismic, ocean acoustic, and tsunami wavefields generated by subduction zone earthquakes, we have developed a provably stable and accurate finite difference method that couples an elastic solid to a compressible fluid subject to gravitational restoring forces. We introduce a linearized dynamic traction-free boundary condition for the moving sea surface that is valid for small amplitude perturbations about an ocean initially in hydrostatic balance. We derive an energy balance for the continuous problem and then use high-order summation-by-parts finite difference operators and weak enforcement of boundary conditions to derive an equivalent discrete energy balance. The discrete energy balance is used to prove stability of the numerical scheme, and stability and accuracy are verified through convergence tests. The method is applied to study tsunami generation by time-dependent rupture on a thrust fault in an elastic solid beneath a compressible ocean. We compare the sea surface evolution in our model to that predicted by the standard tsunami modeling procedure, which assumes seafloor uplift occurs instantaneously and neglects compressibility of the ocean. We find that the leading shoreward-traveling tsunami wave in our model has a noticeably smaller amplitude than that predicted by the standard approach.

A multidomain discretization of the Richards equation in layered soil

Abstract: We consider the Richards equation on a domain that is decomposed into nonoverlapping layers, i.e., the decomposition has no cross points. We assume that the saturation and permeability functions are space-independent on each subdomain. Kirchhoff transformation of each subdomain problem separately then leads to a set of semilinear equations, which can each be solved efficiently using monotone multigrid. The transformed subdomain problems are coupled by nonlinear continuity and flux conditions. This nonlinear coupled problem can be solved using substructuring methods like the Dirichlet–Neumann or Robin iteration. We give several numerical examples showing the discretization error, the solver robustness under variations of the soil parameters, and a hydrological example with four soil layers and surface water.

Benchmark problems for reactive transport modeling of the generation and attenuation of acid rock drainage

Abstract: Acid rock drainage (ARD) is a problem of international relevance with substantial environmental and economic implications. Reactive transport modeling has proven a powerful tool for the process-based assessment of metal release and attenuation at ARD sites. Although a variety of models has been used to investigate ARD, a systematic model intercomparison has not been conducted to date. This contribution presents such a model intercomparison involving three synthetic benchmark problems designed to evaluate model results for the most relevant processes at ARD sites. The first benchmark (ARD-B1) focuses on the oxidation of sulfide minerals in an unsaturated tailing impoundment, affected by the ingress of atmospheric oxygen. ARD-B2 extends the first problem to include pH buffering by primary mineral dissolution and secondary mineral precipitation. The third problem (ARD-B3) in addition considers the kinetic and pH-dependent dissolution of silicate minerals under low pH conditions. The set of benchmarks was solved by four reactive transport codes, namely CrunchFlow, Flotran, HP1, and MIN3P. The results comparison focused on spatial profiles of dissolved concentrations, pH and pE, pore gas composition, and mineral assemblages. In addition, results of transient profiles for selected elements and cumulative mass loadings were considered in the intercomparison. Despite substantial differences in model formulations, very good agreement was obtained between the various codes. Residual deviations between the results are analyzed and discussed in terms of their implications for capturing system evolution and long-term mass loading predictions.

A novel method for NMR data compression

Abstract: Nuclear magnetic resonance (NMR) technique has been widely used to reservoir evaluation and core analysis in oil industry. Rapid and stable inversion of NMR data is very important for NMR logging application. A rapid data compression method with high compression ratio can effectively improve the inversion speed of NMR data. This paper compared and analyzed the window averaging (WA) and singular value decomposition (SVD) methods for NMR data compression. The numerical results show that the WA method has the features of low compression ratio, low computational complexity, and less time-consuming; the SVD method has the features of high compression ratio, large amount of calculation, and more time-consuming. Combining the advantages of these two compression methods, this paper proposed a novel compression method which had achieved good application effects in NMR data compression. The novel method can not only ensure the high compression ratio, but also effectively reduces the time-consuming and possesses more prominent advantages in multi-dimensional NMR data compression.