Authors:Véronique Gervais; Mickaële Le Ravalec Pages: 3 - 28 Abstract: Numerical representations of a target reservoir can help to assess the potential of different development plans. To be as predictive as possible, these representations or models must reproduce the data (static, dynamic) collected on the field. However, constraining reservoir models to dynamic data – the history-matching process – can be very time consuming. Many uncertain parameters need to be taken into account, such as the spatial distribution of petrophysical properties. This distribution is mostly unknown and usually represented by millions of values populating the reservoir grid. Dedicated parameterization techniques make it possible to investigate many spatial distributions from a small number of parameters. The efficiency of the matching process can be improved from the perturbation of specific regions of the reservoir. Distinct approaches can be considered to define such regions. For instance, one can refer to streamlines. The leading idea is to identify areas that influence the production behavior where the data are poorly reproduced. Here, we propose alternative methods based on connectivity analysis to easily provide approximate influence areas for any fluid-flow simulation. The reservoir is viewed as a set of nodes connected by weighted links that characterize the distance between two nodes. The path between nodes (or grid blocks) with the lowest cumulative weight yields an approximate flow path used to define influence areas. The potential of the approach is demonstrated on the basis of 2D synthetic cases for the joint integration of production and 4D saturation data, considering several formulations for the weights attributed to the links. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9663-y Issue No:Vol. 22, No. 1 (2018)

Authors:Andreas S. Stordal; Geir Nævdal Pages: 29 - 41 Abstract: Randomized maximum likelihood is known in the petroleum reservoir community as a Bayesian history matching technique by means of minimizing a stochastic quadratic objective function. The algorithm is well established and has shown promising results in several applications. For linear models with linear observation operator, the algorithm samples the posterior density accurately. To improve the sampling for nonlinear models, we introduce a generalized version in its simplest form by re-weighting the prior. The weight term is motivated by a sufficiency condition on the expected gradient of the objective function. Recently, an ensemble version of the algorithm was proposed which can be implemented with any simulator. Unfortunately, the method has some practical implementation issues due to computation of low rank pseudo inverse matrices and in practice only the data mismatch part of the objective function is maintained. Here, we take advantage of the fact that the measurement space is often much smaller than the parameter space and project the prior uncertainty from the parameter space to the measurement space to avoid over fitting of data. The proposed algorithms show good performance on synthetic test cases including a 2D reservoir model. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9664-x Issue No:Vol. 22, No. 1 (2018)

Authors:M. Baumann; R. Astudillo; Y. Qiu; E. Y. M. Ang; M. B. van Gijzen; R.-É. Plessix Pages: 43 - 61 Abstract: In this work, we present a new numerical framework for the efficient solution of the time-harmonic elastic wave equation at multiple frequencies. We show that multiple frequencies (and multiple right-hand sides) can be incorporated when the discretized problem is written as a matrix equation. This matrix equation can be solved efficiently using the preconditioned IDR(s) method. We present an efficient and robust way to apply a single preconditioner using MSSS matrix computations. For 3D problems, we present a memory-efficient implementation that exploits the solution of a sequence of 2D problems. Realistic examples in two and three spatial dimensions demonstrate the performance of the new algorithm. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9667-7 Issue No:Vol. 22, No. 1 (2018)

Authors:Abiola D. Obembe; M. Enamul Hossain; Ben-Mansour Rached Pages: 63 - 80 Abstract: The Oberbeck-Boussinesq (OB) approximation is widely employed as a simplifying assumption for density-dependent flow problems. It reduces the governing differential equations to simpler forms, which can be handled analytically or numerically. In this study, a modified OB model is formulated to account for the variation of rock permeability and porosity with temperature during the hot fluid injection process in an oil-saturated porous medium under the assumption of local thermal equilibrium (LTE). The mathematical model is solved numerically using a fully implicit control volume finite difference discretization with the successive over relaxation (SOR) method to handle the non-linearity. Subsequently, the numerical model is validated with the analytical solution of the simplified problem successfully. Through detailed sensitivity analyses, the simulation results reveal the hot fluid injection rate as the most important operational parameter to be optimized for a successful thermal flood. The numerical runs show that that for single-phase core-flood simulation, the effect of temperature on the rock absolute permeability and porosity can be neglected without introducing any significant errors in the estimated recovery and temperature profile. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9670-z Issue No:Vol. 22, No. 1 (2018)

Authors:Sujit K Bose Pages: 81 - 86 Abstract: Free surface flow of an incompressible fluid over a shallow plane/undulating horizontal bed is characteristically turbulent due to disturbances generated by the bed resistance and other causes. The governing equations of such flows in one dimension, for finite amplitude of surface elevation over the bed, are the Continuity Equation and a highly nonlinear Momentum Equation of order three. The method developed in this paper introduces the “discharge” variable q = η U, where η = elevation of the free surface above the bed level, and U = average stream-wise forward velocity. By this substitution, the continuity equation becomes a linear first-order PDE and the momentum equation is transformed after introduction of a small approximation in the fifth term. Next, it is shown by an invertibility argument that q can be a function of η: q = F(η), rendering the momentum equation as a first order, second degree ODE for F(η), that can be be integrated by the Runge-Kutta method. The continuity equation then takes the form of a first order evolutionary PDE that can be integrated by a Lax-Wendroff type of scheme for the temporal evolution of the surface elevation η. The method is implemented for two particular cases: when the initial elevation is triangular with vertical angle of 120 ∘ and when it has a sinusoidal form. The computations exhibit the physically interesting feature that the frontal portion of the propagating wave undergoes a sharp jump followed by tumbling over as a breaker. Compared to other discretization methods, the application of the Runge-Kutta and an extended version of the Lax-Wendroff scheme is much easier. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9671-y Issue No:Vol. 22, No. 1 (2018)

Authors:P. Sochala; F. De Martin Pages: 125 - 144 Abstract: A polynomial chaos (PC) surrogate is proposed to reconstruct seismic time series in one-dimensional (1D) uncertain media. Our approach overcomes the deterioration of the PC convergence rate during long time integration. It is based on a double decomposition of the signal: a damped harmonic decomposition combined with a polynomial chaos expansion of the four coefficients of each harmonic term (amplitude, decay constant, pulsation, and phase). These PC expansions are obtained through the least squares method which requires the solution of nonlinear least squares problems for each sample point of the stochastic domain. The use of the surrogate is illustrated on vertically incident plane waves traveling in 1D layered, vertically stratified, isotropic, viscoelastic soil structure with uncertainties in the geological data (geometry, wave velocities, quality factors). Computational tests show that the stochastic coefficients can be efficiently represented with a low-order PC expansion involving few evaluations of the direct model. For the test cases, a global sensitivity analysis is performed in time and frequency domains to investigate the relative impact of the random parameters. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9677-5 Issue No:Vol. 22, No. 1 (2018)

Authors:Ahmad Jan; Ethan T. Coon; Scott L. Painter; Rao Garimella; J. David Moulton Pages: 163 - 177 Abstract: Integrated surface/subsurface models for simulating the thermal hydrology of permafrost-affected regions in a warming climate have recently become available, but computational demands of those new process-rich simu- lation tools have thus far limited their applications to one-dimensional or small two-dimensional simulations. We present a mixed-dimensional model structure for efficiently simulating surface/subsurface thermal hydrology in low-relief permafrost regions at watershed scales. The approach replaces a full three-dimensional system with a two-dimensional overland thermal hydrology system and a family of one-dimensional vertical columns, where each column represents a fully coupled surface/subsurface thermal hydrology system without lateral flow. The system is then operator split, sequentially updating the overland flow system without sources and the one-dimensional columns without lateral flows. We show that the app- roach is highly scalable, supports subcycling of different processes, and compares well with the corresponding fully three-dimensional representation at significantly less computational cost. Those advances enable recently developed representations of freezing soil physics to be coupled with thermal overland flow and surface energy balance at scales of 100s of meters. Although developed and demonstrated for permafrost thermal hydrology, the mixed-dimensional model structure is applicable to integrated surface/subsurface thermal hydrology in general. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9679-3 Issue No:Vol. 22, No. 1 (2018)

Authors:Vadym Aizinger; Andreas Rupp; Jochen Schütz; Peter Knabner Pages: 179 - 194 Abstract: We present an a priori stability and convergence analysis of a new mixed discontinuous Galerkin scheme applied to the instationary Darcy problem. The analysis accounts for a spatially and temporally varying permeability tensor in all estimates. The proposed method is stabilized using penalty terms in the primary and the flux unknowns. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9682-8 Issue No:Vol. 22, No. 1 (2018)

Authors:Shahid Manzoor; Michael G. Edwards; Ali H. Dogru; Tareq M. Al-Shaalan Pages: 195 - 230 Abstract: Grid generation for reservoir simulation must honor classical key constraints and be boundary aligned such that control-volume boundaries are aligned with geological features such as layers, shale barriers, fractures, faults, pinch-outs, and multilateral wells. An unstructured grid generation procedure is proposed that automates control-volume and/or control point boundary alignment and yields a PEBI-mesh both with respect to primal and dual (essentially PEBI) cells. In order to honor geological features in the primal configuration, we introduce the idea of protection circles, and to generate a dual-cell feature based grid, we construct halos around key geological features. The grids generated are employed to study comparative performance of cell-centred versus cell-vertex control-volume distributed multi-point flux approximation (CVD-MPFA) finite-volume formulations using equivalent degrees of freedom. The formulation of CVD-MPFA schemes in cell-centred and cell-vertex modes is analogous and requires switching control volume from primal to dual or vice versa together with appropriate data structures and boundary conditions. The relative benefits of both types of approximation, i.e., cell-centred versus vertex-centred, are made clear in terms of flow resolution and degrees of freedom required. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9686-4 Issue No:Vol. 22, No. 1 (2018)

Authors:Shahid Manzoor; Michael G. Edwards; Ali H. Dogru; Tareq M. Al-Shaalan Pages: 231 - 231 Abstract: Due to an oversight, some author’s corrections were not carried out during Performing proof corrections stage. The Publisher apologizes for these mistakes. The original article was corrected. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9699-z Issue No:Vol. 22, No. 1 (2018)

Authors:Halvor M. Nilsen; Jan Nordbotten; Xavier Raynaud Pages: 233 - 260 Abstract: In this paper, we study newly developed methods for linear elasticity on polyhedral meshes. Our emphasis is on applications of the methods to geological models. Models of subsurface, and in particular sedimentary rocks, naturally lead to general polyhedral meshes. Numerical methods which can directly handle such representation are highly desirable. Many of the numerical challenges in simulation of subsurface applications come from the lack of robustness and accuracy of numerical methods in the case of highly distorted grids. In this paper, we investigate and compare the Multi-Point Stress Approximation (MPSA) and the Virtual Element Method (VEM) with regard to grid features that are frequently seen in geological models and likely to lead to a lack of accuracy of the methods. In particular, we look at how the methods perform near the incompressible limit. This work shows that both methods are promising for flexible modeling of subsurface mechanics. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9687-3 Issue No:Vol. 22, No. 1 (2018)

Authors:Niloofar Misaghian; Mehdi Assareh; MohammadTaqi Sadeghi Pages: 261 - 282 Abstract: The upscaling process of a high-resolution geostatistical reservoir model to a dynamic simulation grid model plays an important role in a reservoir study. Several upscaling methods have been proposed in order to create balance between the result accuracy and computation speed. Usually, a high-resolution grid model is upscaled according to the heterogeneities assuming single phase flow. However, during injection processes, the relative permeability adjustment is required. The so-called pseudo-relative permeability curves are accepted, if their corresponding coarse model is a good representation of the fine-grid model. In this study, an upscaling method based on discrete wavelet transform (WT) is developed for single-phase upscaling based on the multi-resolution analysis (MRA) concepts. Afterwards, an automated optimization method is used in which evolutionary genetic algorithm is applied to estimate the pseudo-relative permeability curves described with B-spline formulation. In this regard, the formulation of B-spline is modified in order to describe the relative permeability curves. The proposed procedure is evaluated in the gas injection case study from the SPE 10th comparative solution project’s data set which provides a benchmark for upscaling problems [1]. The comparisons of the wavelet-based upscaled model to the high-resolution model and uniformly coarsened model show considerable speedup relative to the fine-grid model and better accuracy relative to the uniformly coarsened model. In addition, the run time of the wavelet-based coarsened model is comparable with the run time of the uniformly upscaled model. The optimized coarse models increase the speed of simulation up to 90% while presenting similar results as fine-grid models. Besides, using two different production/injection scenarios, the superiority of WT upscaling plus relative permeability adjustment over uniform upscaling and relative permeability adjustment is presented. This study demonstrates the proposed upscaling workflow as an effective tool for a reservoir simulation study to reduce the required computational time. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9688-2 Issue No:Vol. 22, No. 1 (2018)

Authors:Zhe Liu; Fahim Forouzanfar Pages: 283 - 296 Abstract: In the development of naturally fractured reservoirs (NFRs), the existence of natural fractures induces severe fingering and breakthrough. To manage the flooding process and improve the ultimate recovery, we propose a numerical workflow to generate optimal production schedules for smart wells, in which the inflow control valve (ICV) settings can be controlled individually. To properly consider the uncertainty introduced by randomly distributed natural fractures, the robust optimization would require a large ensemble size and it would be computationally demanding. In this work, a hierarchical clustering method is proposed to select representative models for the robust optimization in order to avoid redundant simulation runs and improve the efficiency of the robust optimization. By reducing the full ensemble of models into a small subset ensemble, the efficiency of the robust optimization algorithm is significantly improved. The robust optimization is performed using the StoSAG scheme to find the optimal well controls that maximize the net-present-value (NPV) of the NFR’s development. Due to the discrete property of a natural fracture field, traditional feature extraction methods such as model-parameter-based clustering may not be directly applicable. Therefore, two different kinds of clustering-based optimization methods, a state-based (e.g., s w profiles) clustering and a response-based (e.g., production rates) clustering, are proposed and compared. The computational results show that the robust clustering optimization could increase the computational efficiency significantly without sacrificing much expected NPV of the robust optimization. Moreover, the performance of different clustering algorithms varies widely in correspondence to different selections of clustering features. By properly extracting model features, the clustered subset could adequately represent the uncertainty of the full ensemble. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9689-1 Issue No:Vol. 22, No. 1 (2018)

Authors:Calogero B. Rizzo; Felipe P. J. de Barros; Simona Perotto; Luca Oldani; Alberto Guadagnini Pages: 297 - 308 Abstract: We study the applicability of a model order reduction technique to the solution of transport of passive scalars in homogeneous and heterogeneous porous media. Transport dynamics are modeled through the advection-dispersion equation (ADE) and we employ Proper Orthogonal Decomposition (POD) as a strategy to reduce the computational burden associated with the numerical solution of the ADE. Our application of POD relies on solving the governing ADE for selected times, termed snapshots. The latter are then employed to achieve the desired model order reduction. We introduce a new technique, termed Snapshot Splitting Technique (SST), which allows enriching the dimension of the POD subspace and damping the temporal increase of the modeling error. Coupling SST with a modeling strategy based on alternating over diverse time scales the solution of the full numerical transport model to its reduced counterpart allows extending the benefit of POD over a prolonged temporal window so that the salient features of the process can be captured at a reduced computational cost. The selection of the time scales across which the solution of the full and reduced model are alternated is linked to the Péclet number (P e), representing the interplay between advective and dispersive processes taking place in the system. Thus, the method is adaptive in space and time across the heterogenous structure of the domain through the combined use of POD and SST and by way of alternating the solution of the full and reduced models. We find that the width of the time scale within which the POD-based reduced model solution provides accurate results tends to increase with decreasing P e. This suggests that the effects of local-scale dispersive processes facilitate the POD method to capture the salient features of the system dynamics embedded in the selected snapshots. Since the dimension of the reduced model is much lower than that of the full numerical model, the methodology we propose enables one to accurately simulate transport at a markedly reduced computational cost. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9693-5 Issue No:Vol. 22, No. 1 (2018)

Authors:Flávio Dickstein; Paulo Goldfeld; Gustavo T. Pfeiffer; Renan V. Pinto Pages: 309 - 327 Abstract: We propose a new algorithm for solving the history matching problem in reservoir simulation, truncated conjugate gradient (TCG), which involves a model reparameterization based on the factorization of the prior covariance matrix, C M = L L T . We also revisit the LBFGS algorithm, framing it into the same reparametrization, introducing M-LBFGS. We present numerical evidence that this reparameterization has an important regularizing impact on the solution process. We show how TCG and M-LBFGS, as well as TSVD, can be implemented without the need of actually computing the factor L. Our numerical experiments, including the PUNQ-S3 and the Brugge cases, indicate that TCG and M-LBFGS are effective schemes for history matching. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9694-4 Issue No:Vol. 22, No. 1 (2018)

Authors:Maziar Veyskarami; Amir Hossein Hassani; Mohammad Hossein Ghazanfari Pages: 329 - 346 Abstract: The network modeling approach is applied to provide a new insight into the onset of non-Darcy flow through porous media. The analytical solutions of one-dimensional Navier-Stokes equation in sinusoidal and conical converging/diverging throats are used to calculate the pressure drop/flow rate responses in the capillaries of the network. The analysis of flow in a single pore revealed that there are two different regions for the flow coefficient ratio as a function of the aspect ratio. It is found that the critical Reynolds number strongly depends on the pore geometrical properties including throat length, average aspect ratio, and average coordination number of the porous media, and an estimation of such properties is required to achieve more reliable predictions. New criteria for the onset of non-Darcy flow are also proposed to overcome the lack of geometrical data. Although the average aspect ratio is the main parameter which controls the inertia effects, the effect of tortuosity on the onset of non-Darcy flow increases when the coordination number of media decreases. In addition, the higher non-Darcy coefficient does not essentially accelerate the onset of inertial flow. The results of this work can help to better understand how the onset of inertial flow may be controlled/changed by the pore architecture of porous media. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9695-3 Issue No:Vol. 22, No. 1 (2018)

Authors:Carsten Burstedde; Jose A. Fonseca; Stefan Kollet Pages: 347 - 361 Abstract: Regional hydrology studies are often supported by high-resolution simulations of subsurface flow that require expensive and extensive computations. Efficient usage of the latest high performance parallel computing systems becomes a necessity. The simulation software ParFlow has been demonstrated to meet this requirement and shown to have excellent solver scalability for up to 16,384 processes. In the present work, we show that the code requires further enhancements in order to fully take advantage of current petascale machines. We identify ParFlow’s way of parallelization of the computational mesh as a central bottleneck. We propose to reorganize this subsystem using fast mesh partition algorithms provided by the parallel adaptive mesh refinement library p4est. We realize this in a minimally invasive manner by modifying selected parts of the code to reinterpret the existing mesh data structures. We evaluate the scaling performance of the modified version of ParFlow, demonstrating good weak and strong scaling up to 458k cores of the Juqueen supercomputer, and test an example application at large scale. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9696-2 Issue No:Vol. 22, No. 1 (2018)

Authors:Fayadhoi Ibrahima; Hamdi A. Tchelepi; Daniel W. Meyer Pages: 389 - 412 Abstract: In the context of stochastic two-phase flow in porous media, we introduce a novel and efficient method to estimate the probability distribution of the wetting saturation field under uncertain rock properties in highly heterogeneous porous systems, where streamline patterns are dominated by permeability heterogeneity, and for slow displacement processes (viscosity ratio close to unity). Our method, referred to as the frozen streamline distribution method (FROST), is based on a physical understanding of the stochastic problem. Indeed, we identify key random fields that guide the wetting saturation variability, namely fluid particle times of flight and injection times. By comparing saturation statistics against full-physics Monte Carlo simulations, we illustrate how this simple, yet accurate FROST method performs under the preliminary approximation of frozen streamlines. Further, we inspect the performance of an accelerated FROST variant that relies on a simplification about injection time statistics. Finally, we introduce how quantiles of saturation can be efficiently computed within the FROST framework, hence leading to robust uncertainty assessment. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9698-0 Issue No:Vol. 22, No. 1 (2018)

Authors:Mohamed Hayek; Anis Younes; Jabran Zouali; Noura Fajraoui; Marwan Fahs Pages: 413 - 421 Abstract: A new analytical solution is developed for interference hydraulic pumping tests in fractal fractured porous media using the dual-porosity concept. Heterogeneous fractured reservoirs are considered with hydrodynamic parameters assumed to follow power-law functions in radial distance. The developed analytical solution is verified by comparison against a finite volume numerical solution. The comparison shows that the numerical solution converges toward the analytical one when the size of the time step decreases. The applicability of the fractal dual-porosity model is then assessed by investigating the identifiability of the parameters from a synthetic interference pumping test with a set of noisy data using Bayesian parameter inference. The results show that if the storage coefficient in the matrix is fixed, the rest of the parameters can be appropriately inferred; otherwise, the identification of the parameters is faced with convergence problems because of equifinality issues. PubDate: 2018-02-01 DOI: 10.1007/s10596-017-9701-9 Issue No:Vol. 22, No. 1 (2018)