Parallel Simulation of Groundwater Flow in the North China Plain

0. INTRODUCTION

The North China Plain (NCP) is the economic and political center of China; and a high-profile area which is in great shortage of water resources in the world. In the NCP, groundwater is the main source of water supply. However, the overexploitation of groundwater over the past several decades has lead to serious water crisis and environmental problems such as land subsidence, which have received great concerns in past decades. To cope with the crisis, lots of researches have been focusing on the sustainable development and management of groundwater resources in the NCP (Wang S F et al., 2013; Currell et al., 2012; Li et al., 2012; Yang et al., 2012; Zheng et al., 2010; Wang S Q et al., 2009, 2008; Liu et al., 2008; Liu and Xia, 2004). One of the key issues is to establish a numerical model which can properly describe the groundwater flow system. To describe the hydraulic variety of the underlying system, physical domain must be discretized at a fine resolution. Also, large time ranges must be considered in order to determine an asymptotic behavior. However, previous research works often focused on modeling a representative area within the NCP or modeling the entire domain using coarse numerical grids (Shao et al., 2013, 2009; Wang et al., 2009, 2008; Liu et al., 2008). For example, in our previous work, there are only 164 rows and 148 columns in the horizontal direction (4 km×4 km per cell) and two-year (2002–2003) duration in the established NCP flow model (Shao et al., 2013, 2009). The long-term prediction of groundwater flow in the NCP at a fine scale remains a challenge today. One of the reasons is that it takes great effort to obtain large quantities of data for the long-term modeling over such vast area. Secondly, because the observed data gives a rather scarce description of the medium hydraulic properties, the uncertainties associated with the characterization of the entire NCP area are large. Thirdly, the simulation time on traditional single-CPU program can be unacceptable when performing high-resolution simulation over long time span and quantifying the uncertainties coming from lack of data.

In recent years, with the shortage of water resources and related environmental problems aggravated, the State Council of China has enacted the most stringent policy on groundwater resources management (the State Council of China, 2012) and has been putting more emphasis on monitoring and controlling the overexploitation of groundwater in the whole country. With both the increased national investment on groundwater monitoring and the advance in detecting technologies, more data can be obtained at finer scale in the NCP today. Based on the data collected, we have established a detailed NCP flow model with much refined spatial discretization and long time span. Traditionally, the popular program MODFLOW is used to resolve the flow model and provides increased flexibility in the scenarios simulated (Huang et al., 2008). However, to address the issue of uncertainty, the number of numerical simulations should be run several times for sensitivity analysis or parameter variations. Consequently, the computational efforts required are much more intensive than before. To overcome the computational difficulty in solving groundwater flow problems, several research works have been devoted to speeding up computations in MODFLOW for large-scale simulations (Hughes and White, 2013; Ji et al., 2012; Dong and Li, 2009; Huang et al., 2008; Schreuder, 2006).

From the perspective of the underlying architecture, these works can be mainly divided into two types: parallel computing works on the shared memory system and on the distributed memory system. The shared memory system mainly refers to the multi-core workstation and the one with the GPU device. Though the shared memory system is more affordable and available, the parallel scalability was often limited by the memory bandwidth and cache coherence problem. The distributed memory system mainly refers to the Linux clusters and the supercomputers. Parallel computing on distributed memory system can have good scalability. However, in order to obtain optimal performance, it often requires lots of programming efforts and taking good care of the data locality, data communication and load balancing etc. This may not be easy for modeling-oriented scientists and unfortunately only tens of processors can be effectively utilized from the previous works.

To further improve performance and parallel scalability, we have recently developed a scalable parallel flow program, which is done by following the idea of rebuilding the existing software by introducing advanced parallel-computing libraries and putting more emphasis on some critical enabling techniques such as data communication, load balance for modeling large-scale problems by parallel computing (Cheng et al., 2014). This paper describes the application of parallel program in simulating the detailed NCP flow model. The parallel simulation is aimed at speeding up the time-consuming simulation of the NCP model with much refined grids and long time periods. The refined NCP model provides a better understanding of the characterization of the groundwater flow system and constitutes the foundation for the long-term prediction and management of the flow system. This study will focus on (1) establishing a NCP flow model (from Jan. 2001 to Dec. 2010) using more than one million grids, and (2) simulating the basin-scale model using the parallel flow program.

1. METHODOLOGY OF THE PARALLEL IMPLEMENTATION

MODFLOW is a well-known computer program for simulating the movement of groundwater flow through porous media and its related problems. We have developed an efficient parallel program, which preserves the capabilities and functionalities of standard MODFLOW program for modeling typical basin-scale problems by rebuilding MODFLOW 2000 on JASMIN (J Adaptive Structured Meshes applications Infrastructure) (Mo et al., 2010), a parallel framework which is built on top of message passing interface (MPI). Briefly speaking, the rebuilding process is mainly achieved by designing a main program based on JASMIN-specific parallel data structures and implementing JASMIN numerical routines which perform operations on patch (which refers to a group of cells in rectangular region and can be perceived as a subregion of a subdomain). To implement the numerical routines, only slight modifications are required to add to the existing MODFLOW subroutines and several routines are added to invoke the modified subroutines from the main program.

In order to simulate the NCP model using the JASMIN implementation, two additional input files are required in addition to the MODFLOW model files. The first is a coordinates file for describing space-discretized geometric information, which is generated from the MF2K DIS file by a simple utility program. The other is a JASMIN database file, which stores various parameters for the linear solver, time profiling, visualization, and so on. In the detailed NCP model, refinements both on spatial grid discretization and sink/source data are performed based on the parallel Kriging algorithm we have recently developed (Cheng et al., 2014). Since the detailed NCP model involves more than one million grids and a decade of time span, it requires tremendous input data to form the MODFLOW model files. Particularly, large quantities of sink/stress data should be read at the beginning of each stress periods. Consequently, the frequent input of large data impairs the overall performance of the simulation program. In order to decrease both the size of data and the times of input operation, all sink/source data are organized by zones rather than by cells, which can reduce the size of input data. Besides, all sink/source data files are read only once at the beginning of the entire simulation rather than one time of each stress period (what is done in standard MODFLOW program), which can further improve the I/O performance. As shown in Fig. 1, the NCP model data are first read by the master processors and then partitioned to each processor in the initialization stage. In parallel processing, each processor owns one subdomain of the entire NCP domain and is responsible for the finite-differnece discretizations on its own subdomain and the formation of one part of the large linear system. The linear system set up at each time step is then solved by the linear solvers provided in JASMIN. In both the formation and solution of the large linear system, data communications are inevitable to exchanged data among different subdomains (different processors). The data communications involved are facilitated by adding ghost nodes to each patch. The ghost nodes refer to the additional memory space used for data exchange; and this approach has been widely applied in modern parallel-computing programs. To avoid branch calculations and thus the parallel performance loss in handling the constant-head/inactive cells, we present a novel approach in our implementations: the active cells are firstly reordered to form the linear system and then handled by the parallel solver. To balance the workload, an efficient load balancing strategy based on the number of active cells on each processor is presented. For details, please refer to (Cheng et al., 2014) for discussion of the parallel implementation with JASMIN.

Figure  1.  Simple compute flow of the JASMIN implementation.

DownLoad: Full-Size Img PowerPoint

2. BRIEF DESCRIPTION OT THE NORTH CHINA PLAIN

2.1 Basic Information

As shown in Fig. 2, the NCP is located at the eastern part of China, which includes plain areas of Beijing and Tianjin City, Hebei Province and parts of Henan and Shandong provinces. It covers a total area of approximately 139 000 km2 and bordered on the north by the Yanshan Mountain, on the west by the Taihang Mountain, to the southeast by the Yellow River, and to the northeast by the Bohai Gulf (Liu et al., 2008). The climate of the NCP belongs to the continental climate. The average annual temperature is about 10–15 ℃. The average annual precipitation is about 500–600 mm, which accounts for only 335 m³ of renewable water resources per capita per year.

Figure  2.  Sketch map showing the location of the NCP.

DownLoad: Full-Size Img PowerPoint

2.2 Hydrogeologic Structure

The NCP is a typical alluvial plain and its elevation is less than 100 m. From the foot of the mountains to the Bohai Bay, it can be divided into the alluvial flood plain of the western part, the distribution of multiple sediments and hydraulic property. It is traditional to divide the NCP aquifers into four aquifers (Fig. 3), which have a depth of about 10–60 (the first aquifer), 120–270 (the second aquifer), 250–350 (the third aquifer), 550–650 m (the fourth aquifer) below the land surface respectively. Since the first aquifer is well developed and almost depleted in some area, pumping wells have been extracting groundwater from the second aquifer. Consequently, these two aquifers were considered as one single unconfined shallow layer in our NCP model. According to our previous investigation (Shao et al., 2013, 2009), variation range of hydraulic conductivity and specific yield of unconfined shallow layer are 3–300 and 0.03–0.23 m/d, respectively. Variation range of hydraulic conductivity of confined layer is 1–50 m/d and the order of magnitude of specific storage ranges from 10-4 to 10-6. More details on the hydrogeologic structure of the NCP model can be found in (Shao et al., 2013; Zhang and Fei, 2009).

Figure  3.  Baoding-Changzhou-Binzhou hydrogeologic cross section in Fig. 2.

DownLoad: Full-Size Img PowerPoint

2.3 Model Discretization and Boundary Conditions

In this study, the new NCP flow model has a much refined discretizaiton in the horizontal direction, 656 rows and 592 columns with a constant spacing of 1 km, compared to only 164 rows and 148 columns with a constant spacing of 4 km in our previous NCP model. In the vertical direction, the refined model is consistent with our previous model, which was divided into three layers with a variable spacing, according to the sand distribution, hydraulic property and groundwater exploitation. The simulation time is a decade, from Jan. 2001 to Dec. 2010, whereas only two-year duration (from Jan. 2002 to Dec. 2003) was considered in our previous work.

As depicted in Fig. 2, both the northern and western boundaries in the shallow layer are treated as piedmont lateral boundaries, since these are natural boundaries between the plain and the mountain terrains (Yanshan Mountain in the north and Taihang Mountain in the west). The south and southeast edge is the Yellow River, so it is also treated as flow boundary in the shallow layer. The northeast edge is the Bohai Gulf, which is treated as general head boundary in the shallow layer. The boundaries of both the middle layer and the deep layer are treated as no-flow boundaries. For details about establishing the refined NCP model, readers can refer to Shao et al. (2014).

2.4 Recharge and Discharge

The average recharge and discharge of the NCP model in ten years (from Jan. 2001 to Dec. 2010) are shown in Table 1. The main sources of groundwater recharge include the precipitation infiltration, the lateral flow from the mountain front, and the leakage from river, drain and agriculture irrigation. To calculate the precipitation infiltration, the Thiessen Polygon Method (Brassel and Reif, 1979) was used to partition the study area and the precipitation infiltration coefficient method was used to compute each subdomain. The lateral flow from the mountain front was calculated according to Darcy's Law. The river leakage is mainly from the Yellow River. Other minor river leakages are from the Hutuo, Chaobai, Tuhai, Majia and Wei rivers.

Table  1.  Average recharge and discharge budget in the NCP from 2001 to 2010

 | Show Table

DownLoad: CSV

The major groundwater discharge item refers to the groundwater pumping, which includes agricultural, municipal, industrial and ecological usage. For the shallow aquifer, the groundwater pumping is mainly concentrated in the piedmont belt. And for the deep aquifer, it is mainly concentrated in the central plain area of Hebei Province. Another important discharge item is the evapotranspiration, which varies in different areas. The influential factors include the groundwater table, the lithology of unsaturated zone, the ground vegetations and the climate, etc.. According to our previous investigation, the cutoff depth in this area is 4 m. The evapotranspiration rate was calculated by the Averiyanov's Empirical Formula, according to the statistics of meteorological stations of each county (Shao et al., 2013).

3. SIMULATION RESULTS AND PERFORMANCE ANALYSIS

The main objective of this study is to simulate the flow in the NCP, using parallel computing techniques and a multiyear, high-resolution grid model. To estimate the accuracy of the parallel simulation, the solution results were verified against the results from MODFLOW as well as the observed data. To demonstrate the performance enhancements brought by parallel processing, informative comparisons on computational time were made between the parallel program and the MODFLOW program.

3.1 Model Calibration

The parallel code has been verified against the standard MODFLOW 2000 code during the development of the program. We found the parallel code exhibit good accuracy and two codes produce the same or nearly the same results for all the tested steady state/transient flow models (Cheng et al., 2014). The parallel simulation results for the detailed NCP model were examined against the field observed data. To be specific, in each stress period for which the observation data are available, the data obtained from 64 observation wells of shallow layer and 45 observation wells of deep layer within the study area were selected to compare with the simulation results. Here, we only list the comparisons at four representative observation wells of the shallow layer (QBJ-1, QHB-4, QSD-15 and QHN-12) in Fig. 4, since we are mainly concerned with the shallow aquifer. The comparisons indicate that the simulated data are generally in good agreement with the observed data. Through this process, model input parameters, such as hydraulic conductivity, recharge rates and boundary conditions, were adjusted to minimize the discrepancy between the observed and the estimated water levels.

Figure  4.  Observed vs. computed hydrographs at four representative observation wells in the shallow layer (see Fig. 2 for the observation wells' location).

DownLoad: Full-Size Img PowerPoint

3.2 Parallel Performance and Analysis

To investigate the parallel performance enhancement, we run the NCP model with our parallel program and standard MODFLOW program respectively. The AMG (algebraic multigrid) method was used in the parallel program while the MICCG (modified incomplete cholesky preconditioned conjugate gradient) and GMG (geometric multigrid) methods were chosen in the MODFLOW program. The residual tolerance (RCLOSE) and head tolerance (HCLOSE) of these methods were set to 0.01 and 1.0 respectively. The maximum number of Picard iterations and linear iterations were set to 50 and 100. To be fair, all the codes are running on the same Linux workstation which is equipped with 4 AMD OPTERON CPU (2.2 GHz, 12 cores per each) and 64 GB memory. On this workstation, all the numerical codes are compiled by GCC 4.5.1 with -O2 flags.

Figure 5a gives a comparison in the solution time between the single-processor MODFLOW program and the parallel JASMIN implementation. To note, we do not list the time for data input and output here, since they are read and initialized once at the beginning of parameter estimation. From the figure, we found that the single-processor JASMIN program using AMG method is much faster than the GMG-based MODFLOW program but slower than the MICCG-based MODFLOW program. The major differences between these two programs lie in three parts: solving the equations, computing matrix & RHS and filling matrix & RHS. Especially, filling matrix & RHS required once in each Picard iteration is an extra computation introduced in parallelization. Fortunately, this overhead can be successfully parallelized and the negative impact can be neglected when compared with the benefit brought by parallel processing. Computing matrix & RHS refers to the formation of the coefficient matrix and RHS from conductances, sink/source items and boundary conditions. On solving the equations, the JASMIN AMG solver is much faster than MICCG and GMG solvers. This phenomenon is mainly owing to the better convergence rate in AMG: both the Picard iteration number and the total iteration number of AMG are much less than that of MICCG; the Picard iteration number of AMG and GMG are comparable, but the average linear iteration number of AMG is less than that of GMG in each Picard iteration.

Figure  5.  Wall clock time and speedup vs. processor number.

DownLoad: Full-Size Img PowerPoint

Figure 5b demonstrates the relative speedup for the JASMIN program using multiple processors. With the increase in processor number, we found remarkable reductions in total execution time and time for computing matrix & RHS, filling matrix & RHS and solving the equations at first; but the reductions become less obvious when the processor number increases to some extent. Consequently, the deviation between the ideal and observed speedup curve becomes larger. One reason behind this phenomenon is though the computational efforts are shared among processors, the setup and data communication costs remain (e.g., the setup cost for the AMG solver) or increase as the processor number increases. Another reason is the load imbalance becomes serious, since the load balancing strategy based on the active cell number owned by each processor is mainly devoted to balancing the load of solving the equations. Nevertheless, the JASMIN program using 32 cores still can be 6 times faster than the MICCG-based MODFLOW program and 11 times faster than the GMG-based MODFLOW program in the solution time. To deal with the uncertainties and perform parameter variations, the NCP flow model may be resolved tens or hundreds of times; thus using the JASMIN parallel program result in considerable savings in compute time.

4. CONCLUSIONS

This paper presents an application of massive parallel computing techniques to speed up the multiyear simulation of subsurface flow in the NCP. The application utilizes over one-million grids to develop a basin-scale NCP flow model and involves hundreds times of parameter estimation, which requires intensive computational efforts. By taking advantage of the parallel program, the compute time is shared among multiple processors, which enables simulating problems with much refined spatial and temporal resolution and long time periods. Simulation results conclude that the proposed parallel computing scheme has improved computational efficiency dramatically in simulating the refined NCP model. Currently, efforts are underway to establish a more refined NCP model (500 m×500 m per cell in the horizontal direction), optimize I/O performance as well as add the capability of local grid refinement.