Large-scale numerical simulation of groundwater flow and solute transport in discretely-fractured crystalline bedrock

A large-scale fluid flow and solute transport model was developed for the crystalline bedrock at Olkiluoto Island, Finland, which is considered as potential deep geological repository for spent nuclear fuel. Site characterization showed that the main flow pathways in the low-permeability crystalline bedrock on the island are 13 subhorizontal fracture zones. Compared to other sites investigated in the context of deep disposal of spent nuclear fuel, most deep boreholes drilled at Olkiluoto are not packed-off but are instead left open. These open boreholes intersect the main fracture zones and create hydraulic connections between them, thus modifying groundwater flow. The combined impact of fracture zones and open boreholes on groundwater flow is simulated at the scale of the island. The modeling approach couples a geomodel that represents the fracture zones and boreholes with a numerical model that simulates fluid flow and solute transport. The geometry of the fracture zones that are intersected by boreholes is complex, and the 3D geomodel was therefore constructed with a tetrahedral mesh. The geomodel was imported into the numerical model to simulate a pumping test conducted on Olkiluoto Island. The pumping test simulation demonstrates that fracture-borehole intersections must be accurately discretized, because they strongly control groundwater flow. The tetrahedral mesh provides an accurate representation of these intersections. The calibrated flow model was then used for illustrative scenarios of radionuclide migration to show the impact of fracture zones on solute transport once the boreholes were backfilled. These mass transport simulations constitute base cases for future predictive analyses and sensitivity studies, since they represent key processes to take into consideration for repository performance assessment.     

Upscaling of a dual-permeability Monte Carlo simulation model for contaminant transport in fractured networks by genetic algorithm parameter identification

The transport of radionuclides in fractured media plays a fundamental role in determining the level of risk offered by a radioactive waste repository in terms of expected doses. Discrete fracture networks methods can provide detailed solutions to the problem of modeling the contaminant transport in fractured media. However, within the framework of the performance assessment (PA) of radioactive waste repositories, the computational efforts required are not compatible with the repeated calculations that need to be performed for the probabilistic uncertainty and sensitivity analyses of PA. In this paper, we present a novel upscaling approach, which consists in computing the detailed numerical fractured flow and transport solutions on a small scale and use the results to derive the equivalent continuum parameters of a lean, one-dimensional dual-permeability, Monte Carlo simulation model by means of a genetic algorithm search. The proposed upscaling procedure is illustrated with reference to a realistic case study of $ {}^{239}{\text{Pu}} $ migration taken from literature.     

Upscaling of flow in heterogeneous porous formations: Critical examination and issues of principle

We consider heterogeneous media whose properties vary in space and particularly aquifers whose hydraulic conductivity K may change by orders of magnitude in the same formation. Upscaling of conductivity in models of aquifer flow is needed in order to reduce the numerical burden, especially when modeling flow in heterogeneous aquifers of 3D random structure. Also, in many applications the interest is in average values of the dependent variables over scales larger or comparable to the conductivity length scales. Assigning values of the conductivity K_b to averaging domains, or computational blocks, is the topic of a large body of literature, the problem being of wide interest in various branches of physics and engineering. It is clear that upscaling causes loss of information and at best it can render a good approximation of the fine scale solution after averaging it over the blocks.The present article focuses on upscaling approaches dealing with random media. It is not meant to be a review paper, its main scope being to elucidate a few issues of principle and to briefly discuss open questions. We show that upscaling can be usually achieved only approximately, and the result may depend on the particular upscaling scheme adopted. The typically scarce information on the statistical structure of the fine-scale conductivity imposes a strong limitation to the upscaling problem. Also, local upscaling is not possible in nonuniform mean flows, for which the upscaled conductivity tensor is generally nonlocal and it depends on the domain geometry and the boundary conditions. These and other limitations are discussed, as well as other open topics deserving further investigation.     

Time-of-flight (TOF)-based two-phase upscaling for subsurface flow and transport

Subsurface formations are characterized by heterogeneity over multiple length scales, which can have a strong impact on flow and transport. In this paper, we present a new upscaling approach, based on time-of-flight (TOF), to generate upscaled two-phase flow functions. The method focuses on more accurate representations of local saturation boundary conditions, which are found to have a dominant impact (in comparison to the pressure boundary conditions) on the upscaled two-phase flow models. The TOF-based upscaling approach effectively incorporates single-phase flow and transport information into local upscaling calculations, accounting for the global flow effects on saturation, as well as the local variations due to subgrid heterogeneity. The method can be categorized into quasi-global upscaling techniques, as the global single-phase flow and transport information is incorporated in the local boundary conditions. The TOF-based two-phase upscaling can be readily integrated into any existing local two-phase upscaling framework, thus more flexible than local–global two-phase upscaling approaches developed recently. The method was applied to permeability fields with different correlation lengths and various fluid-mobility ratios. It was shown that the new method consistently outperforms existing local two-phase upscaling techniques, including recently developed methods with improved local boundary conditions (such as effective flux boundary conditions), and provides accurate coarse-scale models for both flow and transport.     

Upscaling of field-scale soil moisture measurements using distributed land surface modeling

Accurate coarse-scale soil moisture information is required for robust validation of current- and next-generation soil moisture products derived from spaceborne radiometers. Due to large amounts of land surface and rainfall heterogeneity, such information is difficult to obtain from existing ground-based networks of soil moisture sensors. Using ground-based field data collected during the Soil Moisture Experiment in 2002 (SMEX02), the potential for using distributed modeling predictions of the land surface as an upscaling tool for field-scale soil moisture observations is examined. Results demonstrate that distributed models are capable of accurately capturing a significant level of field-scale soil moisture heterogeneity observed during SMEX02. A simple soil moisture upscaling strategy based on the merger of ground-based observations with modeling predictions is developed and shown to be more robust during SMEX02 than upscaling approaches that utilize either field-scale ground observations or model predictions in isolation.     

A simulation of large‐scale groundwater flow and travel time in a fractured rock environment for waste disposal purposes

Two‐dimensional regional groundwater flow was simulated based on a conceptual model of low‐permeability crystalline rocks of the Whiteshell Research Area (WRA) in south‐eastern Manitoba. The conceptual model consists of fracture zones that strike in different directions and dip at various angles in the background rock mass. The thickness and hydraulic properties of the fracture zones in the conceptual model were varied as were the fluid properties and the boundary conditions of the groundwater flow system. The effects of these variations on the groundwater flow pattern and on the convective travel time along pathways from a hypothetical disposal vault at 500 m depth to discharge locations at the ground surface were evaluated. The vault was located in the regional discharge area of the groundwater system. A homogeneous conceptual model of the WRA, having only freshwater flow, formed a groundwater flow pattern with a regional flow system. Local flow systems developed increasingly with the introduction of fracture zones 20 m and 3 m thick, and depth‐dependent fluid density. This indicates a reduction in groundwater residence time by fracture zones and fluid density. Flow pathways were analysed using both a stream‐function and a particle‐tracking technique. The pathways and their lengths from the location of the vault to the surface varied spatially according to the flow patterns. The minimum travel time along these pathways was less than 150 000 and greater than 4 000 000 years in models with and without fracture zones, respectively, indicating that the presence of fracture zones was the major controlling factor. A precise knowledge and refinement of conceptual model parameters is necessary during site selection for waste disposal purposes. Copyright © 2004 John Wiley & Sons, Ltd.     

A parallel mesh-free contaminant transport model based on the Analytic Element and Streamline Methods

This paper introduces a new method for simulating large-scale subsurface contaminant transport that combines an Analytic Element Method (AEM) groundwater flow solution with a split-operator Streamline Method for modeling reactive transport. The key feature of the method is the manner in which the vertically integrated AEM flow solution is used to construct three-dimensional particle tracks that define the geometry of the Streamline Method. The inherently parallel nature of the algorithm supports the development of reactive transport models for spatial domains much larger than current grid-based methods. The applicability of the new approach is verified for cases with negligible transverse dispersion through comparisons to analytic solutions and existing numerical solutions, and parallel performance is demonstrated through a realistic test problem based on the regional-scale transport of agricultural contaminants from spatially distributed sources.     

The role of hand calculations in ground water flow modeling

Most ground water modeling courses focus on the use of computer models and pay little or no attention to traditional analytic solutions to ground water flow problems. This shift in education seems logical. Why waste time to learn about the method of images, or why study analytic solutions to one-dimensional or radial flow problems? Computer models solve much more realistic problems and offer sophisticated graphical output, such as contour plots of potentiometric levels and ground water path lines. However, analytic solutions to elementary ground water flow problems do have something to offer over computer models: insight. For instance, an analytic one-dimensional or radial flow solution, in terms of a mathematical expression, may reveal which parameters affect the success of calibrating a computer model and what to expect when changing parameter values. Similarly, solutions for periodic forcing of one-dimensional or radial flow systems have resulted in a simple decision criterion to assess whether or not transient flow modeling is needed. Basic water balance calculations may offer a useful check on computer-generated capture zones for wellhead protection or aquifer remediation. An easily calculated "characteristic leakage length" provides critical insight into surface water and ground water interactions and flow in multi-aquifer systems. The list goes on. Familiarity with elementary analytic solutions and the capability of performing some simple hand calculations can promote appropriate (computer) modeling techniques, avoids unnecessary complexity, improves reliability, and is likely to save time and money. Training in basic hand calculations should be an important part of the curriculum of ground water modeling courses.     

Reactive transport upscaling at the Darcy scale: A new flow rate based approach raises the unsolved issue of porosity upscaling

This article demonstrates that permeability upscaling, which can require complex techniques, is not necessary to significantly decrease the CPU time in reactive transport modeling. CPU time depends more on the geochemistry than the flow calculation. Flow rate upscaling is proposed as an alternate method to permeability upscaling, which is more suited to time-consuming flow resolution. To apply this method, a finite volume approach is most convenient.Considering the equality of flow as the equivalence criterion, when the coarse grid overlays the fine grid, flow rate upscaling leads, by construction, to the exact results, whereas the accuracy of permeability upscaling methods often depends on specific conditions. Some focus is put on the limitations of a common permeability upscaling technique, the simplified renormalization. In stationary flow, the gain in CPU time is the same for both flow rate upscaling and permeability upscaling. In transient flow, flow rate upscaling is slightly less time-efficient but the ratio between both CPU times decreases when the geochemistry is more complex.Working with an accurate flow rate field in the upscaled case reveals that porosity upscaling is a surprisingly tricky issue. Solution mixing is induced and residence times can be significantly affected. These changes have potentially important consequences on reactive transport modeling. They are not specific to the flow rate upscaling method; they are a general issue. Some simplified cases, assuming a homogeneous mineralogy, are examined. At this stage, a simple heuristic method is proposed, which yields reliable results under particular conditions (high heterogeneity). Porosity upscaling remains an open research field.     

A test case for the simulation of three-dimensional variable-density flow and solute transport in discretely-fractured porous media

A test case has been developed for three-dimensional simulations of variable-density flow and solute transport in discretely-fractured porous media. The simulation domain is a low-permeability porous matrix cube containing a single non-planar fracture. The initial solute concentration is zero everywhere. A constant solute concentration is assigned to the top of the domain, which increases near-top fluid density and induces downward density-driven flow. The test case is therefore comparable to downwelling of a dense brine below a saline disposal basin or a waste repository. Numerous fingers and distinct convection cells develop early in the fracture but the fingers later coalesce and convection becomes less apparent. To help test other variable-density flow and transport models, results of the test case are presented both qualitatively (concentration contours and velocity fields) and quantitatively (penetration depth, mass flux, total mass stored, maximum fracture and matrix velocity).     

Fractured reservoirs with fracture corridors

Outcrop studies reveal a common occurrence of tabular zones of significantly‐increased fracture intensity affecting otherwise well‐lithified rocks. These zones, called fracture corridors, can have a profound effect on multi‐phase fluid flow in the subsurface. Using standard geo‐modelling tools, it is possible to generate 3D realizations of reservoirs that contain distributions of such fracture corridors that are consistent with observations, including the vertical frequency in pseudo‐wells inserted into the model at random locations. These models can generate the inputs to flow simulation. The approach adopted here is to run the flow simulations in a single‐porosity representation where the flow effects of fractures are upscaled into equivalent cell‐based properties, preserving a clear spatial relationship between the input geology and the resulting cellular model. The simulated reservoir performance outcomes are very similar to those seen in real oilfields: extreme variability between wells, early water breakthrough, disappointing recoveries and patchy saturation distributions. Thus, a model based on fracture corridors can provide an explanation for the observed flow performance of a suitable field. However, the use of seismics to identify fracture corridors is not an easy task. New work is needed to predict the seismic responses of fracture corridor systems to be able to judge whether it is likely that we can robustly detect and characterize these flow‐significant features adequately.     

Predicting characteristics of dune bedforms using PSO-LSSVM

Dunes have a large influence on hydraulic roughness, and, thereby, on water levels which could affect the navigability of rivers and performance of hydraulic structures. The present study investigated the variation of geometric and topographic characteristics of dune bedforms and flow features as measured in laboratory studies(data sets from laboratory experiments) to estimate the roughness coefficient and characteristics of dune height. The Least Squares Support Vector Machine(LSSVM), which was optimized using Particle Swarm Optimization(PSO), was used as the Meta model approach to predict the values of interest. Developed models were separated into three categories: modeling using flow characteristics,modeling of flow and bedform characteristics, and modeling by using flow and sediment characteristics.It was found that for estimation of the roughness coefficient in open channels with dune bedforms,models developed based on flow and sediment characteristics performed more successfully. The model with input parameters of flow and grain Reynolds numbers(Re and R_b, respectively) and the ratio of the hydraulic radius(R) to the median grain diameter(D_(50)) yields a squared correlation coefficient(R2) of0.8609, a coefficient of determination(DC) of 0.7361,and a root mean square error(RMSE) of 0.0034 for a test series of Manning's roughness coefficient which was the most accurate model. Results proved the key role of flow Reynolds number(Re) values as an input feature for all models predicting the roughness coefficient. Accordingly, classic approaches led to poor results in comparison. On the other hand, results obtained for estimated values of relative dune height led to moderate prediction quality, which albeit,outperformed classic approaches.     

An Open,Object-Based Framework for Generating Anisotropy in Sedimentary Subsurface Models

The spatial distribution of hydraulic properties in the subsurface controls groundwater flow and solute transport. However, many approaches to modeling these distributions do not produce geologically realistic results and/or do not model the anisotropy of hydraulic conductivity caused by bedding structures in sedimentary deposits. We have developed a flexible object-based package for simulating hydraulic properties in the subsurface—the Hydrogeological Virtual Realities (HyVR) simulation package. This implements a hierarchical modeling framework that takes into account geological rules about stratigraphic bounding surfaces and the geometry of specific sedimentary structures to generate realistic aquifer models, including full hydraulic-conductivity tensors. The HyVR simulation package can create outputs suitable for standard groundwater modeling tools (e.g., MODFLOW), is written in Python, an open-source programming language, and is openly available at an online repository. This paper presents an overview of the underlying modeling principles and computational methods, as well as an example simulation based on the Macrodispersion Experiment site in Columbus, Mississippi. Our simulation package can currently simulate porous media that mimic geological conceptual models in fluvial depositional environments, and that include fine-scale heterogeneity in distributed hydraulic parameter fields. The simulation results allow qualitative geological conceptual models to be converted into digital subsurface models that can be used in quantitative numerical flow-and-transport simulations, with the aim of improving our understanding of the influence of geological realism on groundwater flow and solute transport.     

Hazard area and probability of volcanic disruption of the proposed high-level radioactive waste repository at Yucca Mountain, Nevada, USA

Models that calculate the probability that a new volcano or a dike from a nearby eruption will intersect the footprint of the proposed high-level nuclear waste repository are generalized based on a conceptual model developed for the space transportation industry. The proposed hazard area, defined such that every new eruption that occurs there will disrupt the repository, plays a fundamental role in developing probability models. This hazard area is used not only to hedge the uncertainties in predicting patterns of future volcanic activity, but also to account for the characteristics of a new eruption during the post-closure performance period of an underground geologic repository. The paper discusses the advantages of probability comparisons, capabilities of conservativeness measurements and expert-elicitation on model parameters, and the implications to the proposed repository.Paper funded by a contract from the Agency for Nuclear Projects, State of Nevada, USA.     

Incorporation of conceptual and parametric uncertainty into radionuclide flux estimates from a fractured granite rock mass

Detailed numerical flow and radionuclide simulations are used to predict the flux of radionuclides from three underground nuclear tests located in the Climax granite stock on the Nevada Test Site. The numerical modeling approach consists of both a regional-scale and local-scale flow model. The regional-scale model incorporates conceptual model uncertainty through the inclusion of five models of hydrostratigraphy and five models describing recharge processes for a total of 25 hydrostratigraphic–recharge combinations. Uncertainty from each of the 25 models is propagated to the local-scale model through constant head boundary conditions that transfer hydraulic gradients and flow patterns from each of the model alternatives in the vicinity of the Climax stock, a fluid flux calibration target, and model weights that describe the plausibility of each conceptual model. The local-scale model utilizes an upscaled discrete fracture network methodology where fluid flow and radionuclides are restricted to an interconnected network of fracture zones mapped onto a continuum grid. Standard Monte Carlo techniques are used to generate 200 random fracture zone networks for each of the 25 conceptual models for a total of 5,000 local-scale flow and transport realizations. Parameters of the fracture zone networks are based on statistical analysis of site-specific fracture data, with the exclusion of fracture density, which was calibrated to match the amount of fluid flux simulated through the Climax stock by the regional-scale models. Radionuclide transport is simulated according to a random walk particle method that tracks particle trajectories through the fracture continuum flow fields according to advection, dispersion and diffusional mass exchange between fractures and matrix. The breakthrough of a conservative radionuclide with a long half-life is used to evaluate the influence of conceptual and parametric uncertainty on radionuclide mass flux estimates. The fluid flux calibration target was found to correlate with fracture density, and particle breakthroughs were generally found to increase with increases in fracture density. Boundary conditions extrapolated from the regional-scale model exerted a secondary influence on radionuclide breakthrough for models with equal fracture density. The incorporation of weights into radionuclide flux estimates resulted in both noise about the original (unweighted) mass flux curves and decreases in the variance and expected value of radionuclide mass flux.     

Use of upscaled elevation and surface roughness data in two-dimensional surface water models

In this paper, we present an approach that uses a combination of cell-block- and cell-face-averaging of high-resolution cell elevation and roughness data to upscale hydraulic parameters and accurately simulate surface water flow in relatively low-resolution numerical models. The method developed allows channelized features that preferentially connect large-scale grid cells at cell interfaces to be represented in models where these features are significantly smaller than the selected grid size. The developed upscaling approach has been implemented in a two-dimensional finite difference model that solves a diffusive wave approximation of the depth-integrated shallow surface water equations using preconditioned Newton-Krylov methods. Computational results are presented to show the effectiveness of the mixed cell-block and cell-face averaging upscaling approach in maintaining model accuracy, reducing model run-times, and how decreased grid resolution affects errors. Application examples demonstrate that sub-grid roughness coefficient variations have a larger effect on simulated error than sub-grid elevation variations.     

A comparison of performance of several artificial intelligence methods for forecasting monthly discharge time series

Developing a hydrological forecasting model based on past records is crucial to effective hydropower reservoir management and scheduling. Traditionally, time series analysis and modeling is used for building mathematical models to generate hydrologic records in hydrology and water resources. Artificial intelligence (AI), as a branch of computer science, is capable of analyzing long-series and large-scale hydrological data. In recent years, it is one of front issues to apply AI technology to the hydrological forecasting modeling. In this paper, autoregressive moving-average (ARMA) models, artificial neural networks (ANNs) approaches, adaptive neural-based fuzzy inference system (ANFIS) techniques, genetic programming (GP) models and support vector machine (SVM) method are examined using the long-term observations of monthly river flow discharges. The four quantitative standard statistical performance evaluation measures, the coefficient of correlation (R), Nash–Sutcliffe efficiency coefficient (E), root mean squared error (RMSE), mean absolute percentage error (MAPE), are employed to evaluate the performances of various models developed. Two case study river sites are also provided to illustrate their respective performances. The results indicate that the best performance can be obtained by ANFIS, GP and SVM, in terms of different evaluation criteria during the training and validation phases.     

A new upscaling method for fractured porous media

We present a method to determine equivalent permeability of fractured porous media. Inspired by the previous flow-based upscaling methods, we use a multi-boundary integration approach to compute flow rates within fractures. We apply a recently developed multi-point flux approximation Finite Volume method for discrete fracture model simulation. The method is verified by upscaling an arbitrarily oriented fracture which is crossing a Cartesian grid. We demonstrate the method by applying it to a long fracture, a fracture network and the fracture network with different matrix permeabilities. The equivalent permeability tensors of a long fracture crossing Cartesian grids are symmetric, and have identical values. The application to the fracture network case with increasing matrix permeabilities shows that the matrix permeability influences more the diagonal terms of the equivalent permeability tensor than the off-diagonal terms, but the off-diagonal terms remain important to correctly assess the flow field.