A Parallel PCG Solver for Large-Scale Groundwater Flow Simulation Based on OpenMP

This article presents a parallel simulation solver for groundwater flow on CUDA. Preconditioned conjugate gradient (PCG) algorithm is used to solve the large linear systems arising from the finite-difference discretization of three-dimensional groundwater flow problems. CUDA implementing methods for the two most time-consuming operations in PCG, sparse matrix–vector multiplication and vector inner-product, are given. The experimental results show that CUDA can speed up the solving process of the groundwater simulation significantly. 1.8–3.7 speedup can be achieved with different GPUs for a transient groundwater flow problem.

For a regional land subsidence problem, there are difficulties for the traditional finite element method to solve the regional groundwater flow equation involved. A new finite element method - multiscale finite element method (MsFEM) applied to resolve three dimensional groundwater flow problem- is introduced. Three advantages make the MsFEM very useful for groundwater flow simulation in a large region with highly heterogeneous porous media. One advantage is that this method can significantly save computing efforts by coarse grid mesh. The variation of parameters in an element is brought to the multiscale base functions. Second one is that this method can to a large extent avoid large aspect ratios of elements by combining aquifers and aquitards in one element. The third advantage is that the hydraulic heads in different depths of aquitards or aquifers could be rather accurately interpolated using the multiscale base functions and the nodal heads of the element by MsFEM. The third advantage is useful for the subsidence modeling. A numerical example is done to verify the advantages of MsFEM. Finally, the MsFEM is applied to solve the regional land subsidence model of Shanghai, which includes a three-dimensional groundwater flow mathematical model and a one-dimensional subsidence model. The traditional two-step approach is followed in this study, with the hydrodynamics of the pumped multiaquifer system first simulated by a 3-D groundwater flow model and the subsidence then computed with a 1-D subsidence model with the pore pressure field specified as an external distributed source of strength within the porous medium. The simulation period is from March, 1986 to December, 1998. Using MsFEM, the aquifer system in Shanghai is divided into six layers, which means one aquifer and one aquitard are in one layer with an average thickness of about 50 meters. Every element in MsFEM is subdivided into 7 sub-layers in the vertical direction and every horizontal layer in an element is divided into 4 hexahedral elements. So each element is subdivided into 28 hexahedral sub-cells, which allows the heterogeneity of aquifers and aquitards to be considered by multiscale base functions. The aquifers system in total is discretized into 19,996 elements with a total of 25,000 nodes. The hydraulic heads in any depth of aquitards or aquifers could be interpolated using the multiscale base functions and the nodal heads of the element. It is very useful for the settlement calculation for the aquitards or aquifers after hydraulic heads at nodes of elements are obtained by MsFEM. The application of MsFEM significantly saves computing efforts, effectively decreases the aspect ratios of elements and obtains good results, which proves the MsFEM can be applied to solve the real regional 3D groundwater flow problems, and it is an effective numerical method for regional land subsidence modeling.

This report extends verification and validation testing of the NUFT package of codes (Nitao 2000a,b) to the US1P module of code designed for simulation of single-phase mixed or “coupled” saturated and unsaturated (or “variably saturated”) groundwater flow. Importantly, the verification and validation testing of US1P in this report includes performance evaluation of a modified numerical algorithm that was not included in prior quality assurance (QA) of the NUFT package by Carle et al. (2014) for Underground Test Area (UGTA) project activities, which are directed at assessment of radionuclide contamination in groundwater sourced from underground nuclear test locations at the Nevada National Security Site (NNSS). An immediate purpose of this report is to provide QA for ongoing large-scale, three-dimensional (3-D) groundwater flow modeling of transient water levels associated with longterm water supply pumping and underground nuclear testing at the NNSS (Jackson and Fenelon, 2018; Jackson et al., 2021). For UGTA activities, the QA refers to the standards of ASTM (1996) for verification and validation testing of groundwater modeling codes.

A Three-Dimensional Finite Element - Finite Difference Model for Simulating Confined and Unconfined Groundwater Flow

The preconditioned conjugate gradient (PCG) algorithm is a well-known iterative method for solving sparse linear systems in scientific computations. GPU-accelerated PCG algorithms for large-sized problems have attracted considerable attention recently. However, on a specific multi-GPU platform, producing a highly parallel PCG implementation for any large-sized problem requires significant time because several manual steps are involved in adjusting the related parameters and selecting an appropriate storage format for the matrix block that is assigned to each GPU. This motivates us to propose adaptive optimization modeling of PCG on multi-GPUs, which mainly involves the following parts: (1) an optimization multi-GPU parallel framework of PCG and (2) the profile-based optimization modeling for each one of the main components of the PCG algorithm, including vector operation, inner product, and sparse matrix-vector multiplication (SpMV). Our model does not construct a new storage format or kernel but automatically and rapidly generates an optimal parallel PCG algorithm for any problem on a specific multi-GPU platform by integrating existing storage formats and kernels. We take a vector operation kernel, an inner-product kernel, and five popular SpMV kernels for an example to present the idea of constructing the model. Given that our model is general, independent of the problems, and dependent on the resources of devices, this model is constructed only once for each type of GPU. The experiments validate the high efficiency of our proposed model.

Fault zones intensively disturb local hydrogeologic structures, and, consequently, can play a critical role in governing small to large-scale groundwater flow. Extensive studies have focused on the permeability variation along faults in the light of the conduit or barrier function for the deep groundwater flow. However, little attempt has been made to characterize the hydrological functions of near-surface fault zone.Exposed to atmospheric conditions, fault zones are further disturbed by stress relief and chemical weathering, modifying their structure and generally increasing their permeability. Consequently, the fault zone, which functions as a recharge or discharge zone at the near surface, exerts a non-negligible influence on groundwater flow. However, identifying the hydrological function of near-surface fault zone remains challenging when relying solely on conventional, often non-integrated, geophysical or hydrological investigation approaches.This study proposes a multiphysics coupled strategy to understand the groundwater flow regime around near-surface fault zone. The proposed approach is applied to an active reverse fault zone in Kamikita Plain, NE Japan, which extends for 30 km within the recharge zone of the catchment.The proposed multiphysics approach consists of 5 successive steps:2-D Electrical Resistivity Tomography (ERT) Survey: A 2.3 km-long profile crossing the fault zone, consisting of 7 roll-along surveys with a 6-m electrode spacing. Self-Potential Survey: Conducted along the 2.3 km ERT profile. Rock Property Characterization: A 80 m deep borehole was drilled in the fault zone and physical properties were measured. 3-D Groundwater Flow Simulation of the Fault Zone: Utilizing areal hydrogeological data, measured rock properties, and geophysical imaging. Model Validation Process: Using the results from the groundwater flow simulation, electrical conductivity and self-potential responses were calculated, and compared with observed field data. Preliminary results successfully reproduce the overall resistivity signature and the self-potential anomaly (+35 mV) in the fault zone, attributed to local groundwater upwelling. This newly proposed multiphysics approach could be an essential tool to evaluate the groundwater flow in a region including large-size fault zone, which is important for radioactive waste disposal. Furthermore, this approach could also be effective in capturing the local fluid flow circulation for a variety of applications.Acknowledgements: Main part of this research project has been conducted as the regulatory supporting research funded by the Secretariat of the Nuclear Regulation Authority, Japan.

Modeling the movement of hazardous waste in groundwater was identified by the US Department of Energy (DOE) as one of the grand challenges in scientific computation. In recognition of this need, DOE has provided support for a group of scientists from several national laboratories and universities to conduct research and development in groundwater flow and contaminant transport modeling. This group is part of a larger consortium of researchers, collectively referred to as the Partnership in Computational Science (PICS), that has been charged with the task of applying high-performance computational tools and techniques to grand challenge areas identified by DOE. One of the goals of the PICS Groundwater Group is to develop a new three-dimensional groundwater flow and transport code that is optimized for massively parallel computers. An existing groundwater flow code, 3DFEMWATER, was parallelized in order to serve as a benchmark for these new models. The application of P-FEM, the parallelized version of 3DFEMWATER, to a real field site is the subject of this report.

Efficient parallel multigrid based solvers for large scale groundwater flow simulations

Excessive groundwater exploitation in Sirjan plain, located in south of Iran, has caused severe land subsidence. The main purpose of this paper is to determine the hydrogeological parameters of the aquifer system in Sirjan plain. To this end, the groundwater flow was simulated in both steady and transient states using the United States Geological Survey Modular Three-Dimensional Finite-Difference Groundwater Flow Model, MODFLOW. The hydraulic conductivity parameter of the aquifer system was first estimated and calibrated in the steady state by simulating the hydraulic head at 63 piezometric wells. The specific yield and recharge parameters were then calibrated in the transient state by simulating the drawdown within a period of interest, i.e. 23 September 2004 to 22 September 2005. The main innovation of this paper, however, is to use the potential of Synthetic Aperture Radar (SAR) interferometry time-series analysis in the last step of the proposed algorithm. Small Baseline Subset (SBAS) time-series analysis was applied using Envisat Advanced SAR (ASAR) images acquired from 1 June 2004 to 28 September 2010. The merit of SBAS time series is that it enables us to calculate the amount of subsidence in each period of interest while minimizing the temporal decorrelation effect. In the last step, the SBAS time-series analysis results were obtained to extract the accumulated subsidence from 19 October 2004 to 30 August 2005. The aquifer compaction due to this period was calculated in the Interbed Storage Package (IBS1) for MODFLOW in order to calibrate storage coefficient in a trial-and-error manner so as the calculated subsidence and the observed one by SAR interferometry be in a good agreement. The Root Mean Square Error (RMSE) between the simulated subsidence and the observed one is calculated as 0.019 m which is an indication of the high performance of modelling.

In order to accelerate the Multi-Relaxation-Time Lattice Boltzmann Method (MRT-LBM), it was parallelized on multi-GPUs located in one computer. Lattices of the MRT-LBM were distributed to the different GPUs and they were operated on the GPU threads concurrently. In the step streaming, in order to get the neighbour lattice nformation, all lattices were combined and operated on the master GPU. Three-dimensional (3D) groundwater flow was simulated by using the parallelized MRT-LBM and the experimental results showed that almost 95 times acceleration was attained on a six GPUs computer. The parallelized MRT-LBM based on the multi-GPUs located in one computer is efficient in terms of both time and energy.

Modeling Three-Dimensional Groundwater Flow and Advective Contaminant Transport at a Heterogeneous Mountainous Site in Support of Remediation Strategy Quanlin Zhou, Jens T. Birkholzer, Iraj Javandel, and Preston D. Jordan Ernest Orlando Lawrence Berkeley National Laboratory Earth Sciences Division 1 Cyclotron Road, MS 90-1116, Berkeley CA 94720 USA Abstract A calibrated groundwater flow model for a contaminated site can provide substantial information for assessing and improving hydraulic measures implemented for remediation. A three-dimensional transient groundwater flow model was developed for a contaminated mountainous site, at which interim corrective measures were initiated to limit further spreading of contaminants. This flow model accounts for complex geologic units that vary considerably in thickness, slope, and hydrogeologic properties, as well as large seasonal fluctuations of the groundwater table and flow rates. Other significant factors are local recharge from leaking underground storm drains and recharge from steep uphill areas. The zonation method was employed to account for the clustering of high and low hydraulic conductivities measured in a geologic unit. A composite model was used to represent the bulk effect of thin layers of relatively high hydraulic conductivity found within bedrock of otherwise low conductivity. The inverse simulator ITOUGH2 was used to calibrate the model for the distribution of rock properties. The model was initially calibrated using data collected between 1994 and 1996. To check the validity of the model, it was subsequently applied to predicting groundwater level fluctuation and groundwater flux between 1996 and 1998. Comparison of simulated and measured data demonstrated that the model is capable of predicting the complex flow reasonably well. Advective transport was approximated using pathways of particles originating from source areas of the plumes. The advective transport approximation was in good agreement with the trend of contaminant plumes observed over the years. The validated model was then refined to focus on a subsection of the large system. The refined model was subsequently used to assess the efficiency of hydraulic measures implemented for remediation.

The purpose of this paper is to build a three-dimensional groundwater flow model of the affected zone along the Yellow River, and to conduct equilibrium analysis and resource evaluation to the groundwater of the research area according to the simulation results. In the beginning, the groundwater flow numerical simulation model (1999.1-2009.12) is established and verified through the GMS software on the basis of the establishment of hydrogeological conceptual and mathematical models. Results of the simulation model show that the perennial average of the resources of shallow groundwater recharge is 29.32×108m3/a, the average recharge modulus is 22.35×104m3/km2•a and the safe yield of groundwater resource is 27.04×108m3/a.

The role of karst on large-scale groundwater flow is defined for the Chalk aquifer of Upper Normandy (western Paris Basin), France. In the regional context, chalk plateaus occupy the greater part of watersheds and are the main sites of groundwater recharge. Previous studies focused on karstic output systems in the valleys and less on water-level variations in the recharge zones upstream. This study assesses the relevant hydrogeological processes using time-series data (boreholes and springs) recorded along a down-gradient hydrologeological cross-section in two selected watersheds. These hydrological data are interpreted in the framework of previous descriptions of the morphological organization of the study area’s karst network. The results highlight the hydrological role of (1) the input karst (vertical conduits) which drains recharging water, (2) the output karst (sub-horizontal conduits widely developed in the vicinity of valleys in the surface watersheds) which drains the output flows, and (3) the connections between these two (input and output) networks, which control the upstream water levels and allow quick transfer to springs, particularly after strong rainfall events. A conceptual model of the hydrological functioning of this covered karst aquifer is established, which should serve for the structuring and parameterization of a numerical model.

A new class of normalized explicit approximate inverse matrix techniques, based on normalized approximate factorization procedures, for solving sparse linear systems resulting from the finite difference discretization of partial differential equations in three space variables are introduced. A new parallel normalized explicit preconditioned conjugate gradient square method in conjunction with normalized approximate inverse matrix techniques for solving efficiently sparse linear systems on distributed memory systems, using Message Passing Interface (MPI) communication library, is also presented along with theoretical estimates on speedups and efficiency. The implementation and performance on a distributed memory MIMD machine, using Message Passing Interface (MPI) is also investigated. Applications on characteristic initial/boundary value problems in three dimensions are discussed and numerical results are given.

A multi-GPU parallel optimization model for the preconditioned conjugate gradient algorithm