Subsurface Flow Model Research Articles

Evaluation of the objective functions to improve production history matching performance based on fluid flow behaviour in reservoirs

This work deals with the problem of model performance evaluation, which is a challenge in many science branches, including subsurface fluid flow modelling (e.g. reservoir characterisation using history matching algorithms). The issue is posed as a problem of estimating the model performance by comparing the simulation results with observed data. This is accomplished by defining different statistical indicator objective functions ( OF) to calculate the model efficiency, i.e. how model simulation fits observed data. Several deviation-based statistics used in literature as OF are analysed and applied to two synthetic case studies related to fluid flow in reservoir. Inaccuracies that arise from using the deviation-based statistics are discussed in a comparative way to propose an appropriate approach to boost the model accuracy when selecting the best realisation from multiple equally realisations generated during the history matching algorithm. The evaluations reveal that the suitable approach requires the adoption of an OF, that combines lag time with deviation-based statistic ( SSR), to enhance the history matching process. The proposed approach has proven robust in the sense that it is able to provide an accurate and faster history matching algorithm for reservoir characterisation.

Wavelet Reconstruction of Geologic Facies From Nonlinear Dynamic Flow Measurements

The discrete wavelet transform (DWT) that is widely used in compressing natural images is considered for an effective representation of the geological facies in subsurface flow and transport inverse modeling problems. The inference of the heterogeneous hydraulic rock properties from the scattered dynamic measurements of the flow rates and pressures is a frequently encountered ill-posed inverse problem in subsurface characterization. To better pose this inverse problem, the original grid-based description of the spatial facies maps is replaced with a small number of DWT coefficients that are estimated from indirect nonlinear dynamic measurements. The compressed description of the facies in the wavelet domain after removing the unresolvable high-frequency components leads to an inverse problem with fewer parameters to resolve and improved geologic facies continuity. The main difficulty in the application of the DWT to inverse problems is the lack of sufficient data to resolve higher frequency detail coefficients. Prior information and sensitivity of the flow response to variation in the DWT coefficients are used to infer the location and value of the significant DWT coefficients. The results suggest that the large-scale geologic facies description that control the global flow pattern can be successfully inferred from the dynamic measurements in a reduced wavelet domain. While the flow data may contain information about significant DWT coefficients, a limited observability in ill-posed inverse problems may not allow the identification of these coefficients and the corresponding local spatial features. Therefore, an effective exploitation of the space-frequency localization advantage of the wavelets over the Fourier bases may not be available in solving ill-posed inverse problems.

Simulating unsteady flow, anabranching, and hyporheic dynamics in a glacial meltwater stream using a coupled surface water routing and groundwater flow model

Flooding affects ecosystems by transporting water and solutes across aquatic‐terrestrial interfaces, removing nutrient and organic substrate limitations, and spurring biogeochemical activity. Few studies have considered the influence of flooding on surface water–groundwater interactions. This research examines the temporally variable water storage and exchange in a stream in the McMurdo Dry Valleys (MDV) of Antarctica, where diel flood pulses occur due to glacial melt. Several MDV streams display truncated discharge peaks, suggesting water storage between the source glacier and the gauging station. We tested the hypothesis that stream braids and subsurface water storage contribute to the difference between glacial melt and stream outflow hydrographs by constructing a coupled surface water routing and subsurface water flow model. This model routes water into stream braids at high flows and allows this water to infiltrate and return to the stream via subsurface flow paths as flows recede. Our simulation demonstrates the importance of surface–subsurface water interactions in controlling the hydrograph shape. Maximum simulated discharge was sensitive to storage parameters including aquifer depth and the flooding threshold, while minimum discharge was sensitive to hydraulic conductivity. Subsurface storage volume varied by 38% over a diel cycle and stream‐subsurface exchange rates varied from 0 to 0.19 m3 h−1 m−1, with exchange from the stream to the subsurface during high flows, and vice versa at low flows. These results underscore how unsteady flow can increase hyporheic interactions and ecosystem productivity, and provide support for maintaining natural stream morphology and flow regimes.

An unsplit, higher-order Godunov method using quadratic reconstruction for advection in two dimensions

Linear advection of a scalar quantity by a specified velocity field arises in a number of different applications. Important examples include the transport of species and energy in low Mach number models for combustion, atmospheric flows and astrophysics, and contaminant transport in Darcy models of saturated subsurface flow. In this paper, we present a customized finite volume advection scheme for this class of problems that provides accurate resolution for smooth problems while avoiding undershoot and overshoot for nonsmooth profiles. The method is an extension of an algorithm by Bell, Dawson and Shubin (BDS), which was developed for a class of scalar conservation laws arising in porous media flows in two dimensions. The original BDS algorithm is a variant of unsplit, higher-order Godunov methods based on construction of a limited bilinear profile within each computational cell. The new method incorporates quadratic terms in the polynomial reconstruction, thereby reducing the L1 error and better preserving the shape of advected profiles while continuing to satisfy a maximum principle for constant coefficient linear advection. We compare this new method to several other approaches, including the bilinear BDS method and unsplit piecewise parabolic (PPM) methods.

Identification of nonlinearity in rainfall‐flow response using data‐based mechanistic modeling

Data‐based mechanistic (DBM) modeling is an established approach to time series model identification and estimation, which seeks model structures and parameters that are both statistically optimal and consistent with plausible mechanistic interpretations of the system. This paper describes the application of the DBM method to 10 min, relatively high precision, rainfall‐flow data, including observations of both surface flow and subsurface flow. For a generally wet winter period, the preferred surface flow model is nonlinear in flow generation and linear in routing, while the preferred subsurface flow model is linear in flow generation and nonlinear in routing. These models have mechanistic interpretations in terms of mass balance, hydrodynamics, and conceptual flow pathways. The four‐parameter surface and subsurface flow models explain 91% and 96% of the variance of the corresponding observations. Other plausible models were identified but were less parsimonious or were more reliant on prior perceptions. For a wet summer validation period, the models performed as well as in the calibration period; however, when a long dry spell was included, the performance deteriorated. It is speculated that this is because of complex wetting‐drying dynamics and potential nonstationarity of the soil properties that are not sufficiently revealed in the available data. Conceptual models informed by the DBM results matched the DBM model performance for subsurface flow but gave poorer performance for the more complex surface flow responses. It is concluded that the DBM method can identify nonlinearity in both flow generation and routing and provide conceptual insights that can go beyond prior expectations.

Multi-objective optimisation of cost–benefit of urban flood management using a 1D2D coupled model

This paper presents a multi-objective optimisation (MOO) tool for urban drainage management that is based on a 1D2D coupled model of SWMM5 (1D sub-surface flow model) and BreZo (2D surface flow model). This coupled model is linked with NSGA-II, which is an Evolutionary Algorithm-based optimiser. Previously the combination of a surface/sub-surface flow model and evolutionary optimisation has been considered to be infeasible due to the computational demands. The 1D2D coupled model used here shows a computational efficiency that is acceptable for optimisation. This technological advance is the result of the application of a triangular irregular discretisation process and an explicit finite volume solver in the 2D surface flow model. Besides that, OpenMP based parallelisation was employed at optimiser level to further improve the computational speed of the MOO tool. The MOO tool has been applied to an existing sewer network in West Garforth, UK. This application demonstrates the advantages of using multi-objective optimisation by providing an easy-to-comprehend Pareto-optimal front (relating investment cost to expected flood damage) that could be used for decision making processes, without repeatedly going through the modelling-optimisation stage.

Computational issues and applications of line-elements to model subsurface flow governed by the modified Helmholtz equation

Two new approaches are presented for the accurate computation of the potential due to line elements that satisfy the modified Helmholtz equation with complex parameters. The first approach is based on fundamental solutions in elliptical coordinates and results in products of Mathieu functions. The second approach is based on the integration of modified Bessel functions. Both approaches allow evaluation of the potential at any distance from the element. The computational approaches are applied to model transient flow with the Laplace transform analytic element method. The Laplace domain solution is computed using a combination of point elements and the presented line elements. The time domain solution is obtained through a numerical inversion. Two applications are presented to transient flow fields, which could not be modeled with the Laplace transform analytic element method prior to this work. The first application concerns transient single-aquifer flow to wells near impermeable walls modeled with line-doublets. The second application concerns transient two-aquifer flow to a well near a stream modeled with line-sinks.

Theoretical analysis and field evidence of reciprocity gaps during interference pumping tests

We address the question of the reciprocity of drawdowns observed during a sequence of interference pumping tests performed within a porous medium, in the context of different conceptual models of subsurface flow. We provide a generalization of the work by Bruggeman (1972) [5], extending his results obtained in 1972 for Darcian flows in an unbounded, heterogeneous porous medium. We then analyze reciprocity within a dual-continuum conceptualization, where the medium is viewed as being composed of two overlapping continua, representing the porous matrix and embedded fractures, respectively. We show theoretically, and demonstrate numerically, that only drawdowns associated with the fracture continuum display reciprocity under transient flow conditions. Conversely, non-flowing matrix and fracture continua display reciprocity behavior under steady-state conditions. We then provide field evidence of the insurgence of reciprocity gaps by analyzing interference test data from the karstic limestone aquifer of the Hydrogeological Experimental Site (HES) in Poitiers, France. On the basis of our theoretical results and experimental observations, we discuss different interpretations of the observed reciprocity gaps. These include (a) non-linear dependencies of local hydraulic parameters, (b) occurrence of internal boundaries within the domain, (c) inertial effects that develop through open conduits within the rock matrix, (d) modifications of the aquifer properties between subsequent pumping tests, and (e) significant contribution of the matrix pressure in monitored wells when the behavior of the aquifer is conceptualized by a dual-continuum approach.

A Three-Dimensional, Unsplit Godunov Method for Scalar Conservation Laws

Linear advection of a scalar quantity by a specified velocity field arises in a number of different applications. Of particular interest here is the transport of species and energy in low Mach number models for combustion, atmospheric flows, and astrophysics, as well as contaminant transport in Darcy models of saturated subsurface flow. An important characteristic of these problems is that the velocity field is not known analytically. Instead, an auxiliary equation is solved to compute averages of the velocities over faces in a finite volume discretization. In this paper, we present a customized three-dimensional finite volume advection scheme for this class of problems that provides accurate resolution for smooth problems while avoiding undershoot and overshoot for nonsmooth profiles. The method is an extension of an algorithm by Bell, Dawson, and Shubin (BDS), which was developed for a class of scalar conservation laws arising in porous media flows in two dimensions. The original BDS algorithm is a variant of unsplit, higher-order Godunov methods based on construction of a limited bilinear profile within each computational cell. Here we present a three-dimensional extension of the original BDS algorithm that is based on a limited trilinear profile within each cell. We compare this new method to several other unsplit approaches, including piecewise linear methods, piecewise parabolic methods, and wave propagation schemes.

Development and Validation of a Numerical Model of Subsurface Flow Constructed Wetlands

Development and Validation of a Numerical Model of Subsurface Flow Constructed Wetlands

Impact of grid resolution on the integrated and distributed response of a coupled surface–subsurface hydrological model for the des Anglais catchment, Quebec

AbstractDigital elevation models (DEMs) at different resolutions (180, 360, and 720 m) are used to examine the impact of different levels of landscape representation on the hydrological response of a 690‐km2 catchment in southern Quebec. Frequency distributions of local slope, plan curvature, and drainage area are calculated for each grid size resolution. This landscape analysis reveals that DEM grid size significantly affects computed topographic attributes, which in turn explains some of the differences in the hydrological simulations. The simulations that are then carried out, using a coupled, process‐based model of surface and subsurface flow, examine the effects of grid size on both the integrated response of the catchment (discharge at the main outlet and at two internal points) and the distributed response (water table depth, surface saturation, and soil water storage). The results indicate that discharge volumes increase as the DEM is coarsened, and that coarser DEMs are also wetter overall in terms of water table depth and soil water storage. The reasons for these trends include an increase in the total drainage area of the catchment for larger DEM cell sizes, due to aggregation effects at the boundary cells of the catchment, and to a decrease in local slope and plan curvature variations, which in turn limits the capacity of the watershed to transmit water downslope and laterally. The results obtained also show that grid resolution effects are less pronounced during dry periods when soil moisture dynamics are mostly controlled by vertical fluxes of evaporation and percolation. Copyright © 2010 John Wiley & Sons, Ltd.

Infiltration in Arid Environments: Spatial Patterns between Subsurface Heterogeneity and Water‐Energy Balances

Transient simulations of an arid mountain system were conducted using a fully integrated model of subsurface and overland water flow and the land surface energy balance (ParFlow). An hourly atmospheric time series was constructed and used to force a two‐dimensional model containing detailed, stochastic descriptions of subsurface heterogeneity for alluvial and fractured units. Two cases were simulated for multiple years of simulation: a dry case based on a year with below‐average precipitation, and a wet case with above‐average precipitation. For each case, four realizations of hydraulic conductivity were simulated, along with a homogeneous domain and a domain with lower fracture density. Detailed analysis of water balances from each simulation indicated that recharge (change in subsurface water storage) over multiple years was significant even for the dry case. The fractured model of subsurface heterogeneity, along with episodic infiltration, combined to create localized regions of fully saturated (perched) water. This was shown to be even more substantial in the wet case. By contrast, homogeneous simulations did not exhibit this behavior and estimated lower recharge than the heterogeneous simulations. Geospatial statistics were subsequently used to evaluate correlations in land‐energy fluxes. The land‐energy fluxes exhibited clear spatial correlation and were influenced by both the shallow and deeper soil structures. This indicates that land–atmosphere fluxes may provide an integrated measure of subsurface heterogeneity. Finally, the vertical spatial structure of soil saturation in the fractured system was shown to exhibit multiscale behavior.

Real-Time Management of an Urban Groundwater Well Field Threatened by Pollution

We present an optimal real-time control approach for the management of drinking water well fields. The methodology is applied to the Hardhof field in the city of Zurich, Switzerland, which is threatened by diffuse pollution. The risk of attracting pollutants is higher if the pumping rate is increased and can be reduced by increasing artificial recharge (AR) or by adaptive allocation of the AR. The method was first tested in offline simulations with a three-dimensional finite element variably saturated subsurface flow model for the period January 2004-August 2005. The simulations revealed that (1) optimal control results were more effective than the historical control results and (2) the spatial distribution of AR should be different from the historical one. Next, the methodology was extended to a real-time control method based on the Ensemble Kalman Filter method, using 87 online groundwater head measurements, and tested at the site. The real-time control of the well field resulted in a decrease of the electrical conductivity of the water at critical measurement points which indicates a reduced inflow of water originating from contaminated sites. It can be concluded that the simulation and the application confirm the feasibility of the real-time control concept.

A sparse Bayesian framework for conditioning uncertain geologic models to nonlinear flow measurements

We present a Bayesian framework for reconstructing hydraulic properties of rock formations from nonlinear dynamic flow data by imposing sparsity on the distribution of the parameters in a sparse transform basis through Laplace prior distribution. Sparse representation of the subsurface flow properties in a compression transform basis (where a compact representation is often possible) lends itself to a natural regularization approach, i.e. sparsity regularization, which has recently been exploited in solving ill-posed subsurface flow inverse problems. The Bayesian estimation approach presented here allows for a probabilistic treatment of the sparse reconstruction problem and has its roots in machine learning and the recently introduced relevance vector machine algorithm for linear inverse problems. We formulate the Bayesian sparse reconstruction algorithm and apply it to nonlinear subsurface inverse problems where solution sparsity in a discrete cosine transform is assumed. The probabilistic description of solution sparsity, as opposed to deterministic regularization, allows for quantification of the estimation uncertainty and avoids the need for specifying a regularization parameter. Several numerical experiments from multiphase subsurface flow application are presented to illustrate the performance of the proposed method and compare it with the regular Bayesian estimation approach that does not impose solution sparsity. While the examples are derived from subsurface flow modeling, the proposed framework can be applied to nonlinear inverse problems in other imaging applications including geophysical and medical imaging and electromagnetic inverse problem.

Faults and fault properties in hydrocarbon flow models

Geofluids (2010) 10, 94–113AbstractThe petroleum industry uses subsurface flow models for two principal purposes: to model the flow of hydrocarbons into traps over geological time, and to simulate the production of hydrocarbon from reservoirs over periods of decades or less. Faults, which are three‐dimensional volumes, are approximated in both modelling applications as planar membranes onto which predictions of the most important fault‐related flow properties are mapped. Faults in porous clastic reservoirs are generally baffles or barriers to flow and the relevant flow properties are therefore very different to those which are important in conductive fracture flow systems. A critical review and discussion is offered on the work‐flows used to predict and model capillary threshold pressure for exploration fault seal analysis and fault transmissibility multipliers for production simulation, and of the data from which the predictions derive. New flow simulation models confirm that failure of intra‐reservoir sealing faults can occur during a reservoir depressurization via a water‐drive mechanism, but contrary to anecdotal reports, published examples of production‐induced seal failure are elusive. Ignoring the three‐dimensional structure of fault zones can sometimes have a significant influence on production‐related flow, and a series of models illustrating flow associated with relay zones are discussed.

External Iterative Coupling Strategy for Surface-Subsurface Flow Calculations in Surface Irrigation

Coupling the unsteady open-channel flow equations of surface irrigation with the equation of variably saturated porous media flow is a computationally complex problem because of the dependence of infiltration on surface-flow depths. Several models of this coupled process have been developed, all of which solve the surface and subsurface flow iteratively at each time step of the surface-flow solution. This study presents an alternative strategy, in which stand-alone surface and subsurface-flow models are used to conduct the calculations sequentially and iteratively at the time level of the irrigation event. At each iteration, the subsurface-flow results are computed using the surface-flow depths generated by the surface-flow model at the current iteration. Infiltration results computed at selected computational nodes are fitted to an empirical infiltration function, and then fed back to the surface-flow model. The proposed strategy, labeled external iterative coupling, was prototyped for border and basin irrigation systems using the WinSRFR and HYDRUS-1D models. The proposed procedure produced irrigation performance results comparable to those generated with an internally coupled model, even when using a single representative location to calibrate the empirical infiltration equation used by the surface-flow model. In comparison with models that iterate at each time step, the proposed coupling strategy reduces the computational effort and improves convergence. The approach provides a practical alternative for coupling existing and future surface and subsurface flow models.

Comparison of Data Assimilation Techniques for a Coupled Model of Surface and Subsurface Flow

Data assimilation in the geophysical sciences refers to methodologies to optimally merge model predictions and observations. The ensemble Kalman filter (EnKF) is a statistical sequential data assimilation technique explicitly developed for nonlinear filtering problems. It is based on a Monte Carlo approach that approximates the conditional probability densities of the variables of interest by a finite number of randomly generated model trajectories. In Newtonian relaxation or nudging (NN), which can be viewed as a special case of the classic Kalman filter, model variables are driven toward observations by adding to the model equations a forcing term, or relaxation component, that is proportional to the difference between simulation and observation. The forcing term contains four‐dimensional weighting functions that can, ideally, incorporate prior knowledge about the characteristic scales of spatial and temporal variability of the state variable(s) being assimilated. In this study, we examined the EnKF and NN algorithms as implemented for a complex hydrologic model that simulates catchment dynamics, coupling a three‐dimensional finite element Richards' equation solver for variably saturated porous media and a finite difference diffusion wave approximation for surface water flow. We report on the retrieval performance of the two assimilation schemes for a small catchment in Belgium. The results of the comparison show that nudging, while more straightforward and less expensive computationally, is not as effective as the ensemble Kalman filter in retrieving the true system state. We discuss some of the strengths and weaknesses, both physical and numerical, of the NN and EnKF schemes.

Simulating Dynamic Water Supply Systems in a Fully Integrated Surface–Subsurface Flow and Transport Model

Groundwater pumping is a boundary condition in groundwater models, and surface application, reuse, diversions, and water reclamation plant discharges are boundary conditions in hydrologic or hydraulic (i.e., surface) models. These processes can be simulated, however, as internal flow and solute transport components within a fully integrated groundwater–surface water model. A methodology has been developed and incorporated into the MODHMS fully integrated surface–subsurface modeling framework to enable evaluation of water distribution and application and impacts to water quality on surface and subsurface domains. This methodology allows the supply and delivery of water and associated water quality to be a part of the model solution rather than being predetermined boundary conditions to the system. Furthermore, this methodology handles temporal and spatial changes to land use and water quality and associated effects on outdoor water distribution, infiltration, runoff, and evapotranspiration properties that are related to surface cover or land use classifications. This new methodology was applied to a portion of the Upper Santa Clara River system and underlying groundwater basins of southern California to aid in the evaluation of a Cl− total maximum daily load (TMDL) study. The outcome of this modeling evaluation has led to an enhanced understanding of how Cl− is transported throughout the surface–subsurface model domain and improved collaboration among stakeholders. The model serves as a powerful tool with which to develop and evaluate various approaches for the Cl− TMDL implementation for the Upper Santa Clara River system.

Ensemble Kalman filter data assimilation for a process‐based catchment scale model of surface and subsurface flow

A sequential data assimilation procedure based on the ensemble Kalman filter (EnKF) is introduced and tested for a process‐based numerical model of coupled surface and subsurface flow. The model is based on the three‐dimensional Richards equation for variably saturated porous media and a diffusion wave approximation for overland and channel flow. A one‐dimensional soil column experiment and a three‐dimensional tilted v‐catchment test case are presented. A preliminary analysis of the assimilation scheme is undertaken for the one‐dimensional test case in order to validate the implementation by comparison with published results and to assess the influence of various factors on the filter's performance. The numerical results suggest robustness with respect to the ensemble size and provide useful information for the more complex tilted v‐catchment test case. The assimilation frequency and the effects induced by data assimilation on the surface and/or subsurface system states are then evaluated for the v‐catchment experiment using synthetic observations of pressure head and streamflow. The results suggest that streamflow prediction can be improved by assimilation of pressure head and streamflow, either individually or in tandem, whereas assimilation of streamflow data alone does not improve the subsurface system state. In terms of the global system state, i.e., surface and subsurface variables, frequent updates are especially beneficial when assimilating both pressure head and streamflow. Furthermore, it is shown that better evaluation of the subsurface volume resulting from assimilation of head data is crucial for improving subsequent surface response.

FITOVERT: A dynamic numerical model of subsurface vertical flow constructed wetlands

This paper introduces a mathematical model (FITOVERT) specifically developed to simulate the behaviour of vertical subsurface flow constructed wetlands (VSSF-CWs). One of the main goals of the development of FITOVERT was to keep the complexity of the model to an acceptable level, so as to provide a practical tool for design and operation optimization. The dynamic formulation of the model allows to simulate the typical non stationary feeding-emptying operation of VSSF-CWs. FITOVERT is able to describe the water flow through porous media in unsaturated conditions, combined with evapotranspiration; its biochemical module describes the degradation of both organic matter and nitrogen; the transport in the liquid phase is implemented for both dissolved and particulate components; the oxygen transport in the gaseous phase of the soil and its exchange with the liquid phase are also considered. As a main advantage, compared to the few currently available dedicated numerical models, FITOVERT is able to handle the porosity reduction due to bacteria growth and accumulation of particulate components, so that the clogging process is also simulated as an effect of the pore size reduction on the hydraulic conductivity of the simulated system. The performance of the model was firstly analyzed by comparison with hydrodynamic tests recorded in an experimental VSSF-CW pilot plant: tracer test were carried out in three different saturation conditions (fully saturated, partially saturated, and completely drained). FITOVERT proved to accurately simulate the hydraulic behaviour of VSSF-CWs in both saturated and unsaturated conditions. The needs for model improvements and further calibration are finally discussed.