Numerical Simulation of Groundwater Flow in a River Valley Basin in Jilin Urban Area, China

In order to evaluate the groundwater resources and aquifer system of the Jilin urban area (JUA), a groundwater numerical flow model was established by using the groundwater modeling system based on data from 190 boreholes. River stages were interpolated to control the groundwater flow field. The input parameters such as hydraulic conductivity and specific yield were based on data from 260 pumping test data. The model was calibrated by trial and error, simulated results were compared to the observed head and contour maps, which were generally in good agreement, and the root mean squared error was 0.66 m. Sensitivity analysis was carried out and recharge proved to be the most sensitive factor in this model. The water budget showed that the input was 2.07× 10 8 m 3 /a, which was smaller than the output of 2.21 × 10 8 m 3 /a. A groundwater level decline and cone of depression exist in the Songhua and Aolong river valley. The JUA aquifer systems can be well and efficiently modeled by constructing a numerical model. Based on the supply and demand analysis of water resources, the established model would finally provide a scientific basis to use the groundwater resources sustainably in JUA.


Introduction
The Jilin urban area (JUA) is located on the banks of the Songhua River (Figure 1).Surface water is abundant in the JUA and is the main source of water, while groundwater is an auxiliary water supply.Excessive reliance on surface water may result in a domestic water supply crisis in the city if the river becomes polluted.For example, an explosion occurred in 2005 at a petrochemical plant owned by the Jilin Petrochemical Corporation, which caused a large spill of toxic substances consisting of a mixture of benzene, aniline, and nitrobenzene (NB) to contaminate the second Songhua River [1,2].In addition, barrels of trimethyl chlorosilane material were washed into the second Songhua river by a flood in 2010 [3].Groundwater is a relatively stable and essential source of drinking water, especially when compared to surface water.Before using the groundwater, a proper understanding of the groundwater system behavior is very important [4][5][6].
Groundwater simulation modeling is recognized as a very useful tool for understanding regional groundwater systems, management of groundwater resources, and planning for future water consumption [7][8][9].Over the years, various models have been developed for modeling subsurface flow.These include the following: Groundwater Modeling System (GMS) [10,11], developed by U.S. Department of Defense, Visual MODFLOW [10] developed by the Waterloo Hydrogeologic Inc., FEFLOW [11] developed by the US Department of Interior, Groundwater Vistas [12] developed by Environmental Simulation, Inc., Interactive Ground Water (IGW) [13] developed by Michigan State University, SWIFT [14] developed by Sandia National laboratories and so on.GMS which supports groundwater modeling computer code MODFLOW 2000 and FEFLOW, is often used for groundwater model set up and can be used for different types of groundwater flow systems, including confined and unconfined aquifers [15,16], simple or multiple aquifers [4], and homogeneous or heterogeneous layered aquifer systems [17,18].Different types of modeling techniques have already been used by GMS: full [4] and quasi-three dimensional (3D) solutions [19], finite element [20] and finite difference solutions [4]; steady state [20] and transient solutions [21].GMS allows for model conceptualization, mesh and grid generation, geostatistics and post-processing of output [22].Auto CAD and Geographic Information System (GIS) [23][24][25] technology can be integrated with the MODFLOW code in GMS for water resources management [26].
The groundwater flow system for the whole JUA district has not been explored and simulated, and there is little guidance on how to deal with the complicated geological conditions, boundaries, the relationship between the groundwater and rivers, and the construction and calibration of the model for the JUA.Zhang only focus on the GMS groundwater quality in the Mangniuhe district, whereas in the east of the JUA [27], Wang analyzed the quantity of surface water using Mikebasin [28].Meanwhile, the Jilin government began to evaluate and investigate the JUA groundwater resources in order to be better prepared for surface water pollution emergencies and to protect the water supply.Therefore, establishment of a regional groundwater model for the JUA has practical implications for the quantitative evaluation of groundwater resources and the long-term prediction of the groundwater regime in this area.The main purpose of this study was to develop a numerical groundwater flow model to characterize the intermountain river valley groundwater flow system and evaluate the JUA water resources using the GMS software.Based on the supply and demand analysis of water resources, the established model could finally provide a scientific basis for the sustainable use of groundwater resources in the JUA.

Materials and Methods
In this paper, GIS was used for analyzing geographic information while conceptualizing hydrogeological systems.All GIS data required for conceptual modeling were converted to the appropriate formats for numerical modeling using ArcGIS 10.0 (ESRI) [6].The database and groundwater modeling include some stages, which are shown in Figure 2.

Study Area
The JUA covers 1267.89km 2 and is located between 4840,000 and 4905,000 north and 260,000 and 320,000 east in the central of Jilin Province, China (Figure 1).JUA has a population of 1.60 million and is an important irrigation and industrial district of China.JUA belongs to the hilly area of Changbai Mountain, and the Songhua River valley plain.JUA stretches from the Xiaobai, Jian, and Mao mountains in the south to Aolongbei, Jiabanling, and the Tiger mountains in the north.It belongs to the temperate zone continental climate.The mean annual rainfall for the last 30 years (from 1981 to 2010) is 580-670 mm and the average evaporation from land is 367 mm (Figure 3).Average yearly temperature is 3.9 °C, with moderate warm summers in July with an average temperature of 21 °C-23 °C and cold winters in January with an average temperature of −20 °C to −18 °C.The second Songhua River, which is the largest river in Jilin province, flows across JUA in a south-north direction.There are four tributaries (Wende River, Mangniu River, Aolong River, and Tuanshanzi River) of the second Songhua River and 10 smaller branches that flow across the JUA (Figure 1).According to the discharge data of the Jilin hydrologic station from 1954 to 2010, the mean annual discharge of the second Songhua River was 414 m 3 /s.JUA can be divided into irrigation and industrial zones that are divided by Jiuzhan.

Geological and Hydrogeological Condition
The lithological map of JUA (Figure 4) was obtained using the published geological maps and borehole data.An aquifer formed in Quaternary unconsolidated sediments consisting mainly of clays, sands, and gravels is broadly distributed within the JUA.The oldest sediment outcropping in the JUA is the Permian (tuff), and the Cretaceous system bedrock (sandstone).The study area is divided into three hydrogeological zones according to the locations and geomorphic units (Figure 5).Zone I is the valley of the upstream of the second Songhua River, the Wende River and the Mangniu River, Zone II is the valley of the Aolong River, and Zone III is the platform area.The study area can be divided into 37 sub-regions according to the boundary conditions, the water yield properties, the difference of the landforms, the precipitation and irrigation conditions, and the exploitation of the groundwater.The irrigation zone is in the North and the industrial zone is in the South (Figure 5).The aquifer formed by the Quaternary unconsolidated deposits is simulated in this work.The thickness of the Quaternary unconsolidated deposits is variable.Three cross-sections are made (location of the section lines in Figure 1) and presented in Figure 6 to analyze the sequence of the layers.The cross-section of the study area shows that the thickness of the Quaternary deposits is variable, from 20 to 70 m of sand and gravel deposits in alluvial fans and less than 5 m in the upland areas.The thickness of upper unconfined Quaternary aquifer (loam and clayey silt) is 2-35 m. Figure 6.Hydrogeological cross-section of Jilin study area (location of the section lines in Figure 1).

Numerical Model Development
Geological insight and interpretation are crucial, because a conceptual model of the system based only on hydraulic data would be incorrect and misleading [29].Development of a numerical groundwater model usually requires two steps: developing an accurate geological model that integrates all relevant data; and building a groundwater flow model based on the geological model to numerically represent the conceptual groundwater system.

Geological Model
High-quality data are very important for constraining physical and hydrogeological conditions [30].All relevant features, such as topography, digital elevation model (DEM), borehole data, geological structure, drainage density, groundwater flow, and boundary conditions of the aquifer system are very important for setting up the model [31,32].
A total of 190 borehole data points and 75 observation data points were collected, and integrated in GMS to build the geological model.The hydrogeological profiles were drawn to control the interpolation of the top and bottom layers in the conceptual model.There are two hydrogeological systems in the JUA aquifer system.The shallow system is made of silt and the deep one is made of sand and gravel.The hydraulic connection is good between the two systems in the valley region because they are all unconfined aquifers.
The surface elevation is very important in groundwater numerical simulation as the groundwater level, evaporation, and infiltration have a close relationship with the surface elevation.A 1:50,000 topographic map is used to integrate the geological field mapping data, but the contour interval of the 1:50,000 topographic map is 50 m; in the hilly regions it is enough to describe the terrain, but it is rarely enough to describe the terrain in valley districts because of subtle terrain changes.The surface elevation cannot be delineated by only using the borehole data because of the JUA's large area.Therefore, the JUA was delineated using the Shuttle Radar Terrain Mapping (SRTM) digital elevation model (DEM) data with 90 m resolution from the US Geological Survey (USGS) [33,34].ArcGIS 10.0 software was used to reset the resolution of the SRTM DEM.Lambert Conformal Conic was adopted as the projection between SRTM and GMS.The Gauss-Krueger projection and Beijing1954 Coordinate System were selected to integrate the different spatial data [35].ArcGIS was used to analyze, pick up, and section the data and set the contour interval at 10 m.The lines were directly read by GMS and used to scatter the 2D scatters.Then, the 2D contour map was obtained by interpolating the 2D scatters, and the elevations of the 190 boreholes were interpolated to revise the error of the SRTM DEM.
JUA is bounded by the mountains on the south, east, and west and the river on the north.The grid cells were designated as "inactive" and "active", outside and inside, respectively, the model domain.
The boundaries between the shallow and deeper systems are different as the sand and gravel pinch out and the rocks outcrop.Figure 4 shows that the gravel strata is thinning up in the boundaries, and the geologic model will be anamorphic in the boundary using Kriging that was based on the borehole data distributed throughout the valley.Interpolating boreholes are needed for boundary constraint.The hydrogeological cross-section and the borehole near the boundaries were especially important as constraint conditions.
The top and bottom layers and the groundwater table were interpolated using data from the 190 boreholes; artificial interpolating is needed in the districts where the thickness of lithology changes greatly.The geometrical model of JUA is shown in Figure 7a.

Model Discretization
A groundwater flow model in the JUA was built using MODFLOW-2000 in GMS.The area of JUA is 1267.89km 2 and it is assumed to be a rectangular shape.The grid cells were designated as "inactive" outside the model domain and as "active" inside the model domain.The model was divided into 120 rows and 124 columns with the size of each cell being 500 × 500 m.The model includes two layers to represent the hydrogeological condition and data on the porous aquifers at the study site.The simulation period was from 1 January 2001 to 31 December 2010, with 120 stress periods defined for the simulation.Two time steps were defined for each stress period.Groundwater levels observed in 2001 were used as the initial conditions of the model.

Boundary Conditions
The boundary conditions are shown in Figure 7b.For the lateral boundary conditions, the Dirichlet boundary condition and the Neumann boundary condition exist in the model.The river boundary is the Dirichlet boundary condition because there is a good connection between the groundwater and the river.In the interface of hilly regions and at the places perpendicular to the groundwater flow direction, the boundary is defined as a "no flow" boundary as they are impermeable boundaries.Along the ground water flow direction, the boundary is defined as a "specific flow" boundary.The flow rate is based on the hydraulic conductivities (K), the thickness (M), and the hydraulic gradient (I).The flow boundary was defined as wells in MODFLOW.For the vertical boundary conditions, the top boundary was defined for recharge and discharge flux, infiltration of precipitation, irrigation, river leakage, pumping wells, and phreatic water evaporation on the surface of the water table.The bottom boundary condition at the base of the deep aquifer was treated as no-flow boundaries.

Parameters
The formation and distribution of groundwater is controlled by the formation lithology, the geologic structure and the geomorphologic shape.The field geological observations show that the porous aquifer consists of quaternary clay, sands and gravel, and is the main aquifer in the study area.The fissured aquifer before Quaternary is undeveloped and forms a lower confining bed.Pumping test data from more than 260 tests were obtained from previous hydrologic studies and used to quantify hydraulic characteristics of the aquifer and identify the boundary conditions [36].The Theis equation and Cooper-Jacob graphical method were used to calculate parameters such as hydraulic conductivity (K) and specific yield (μ).For the sand and gravel aquifer, the calculated K ranges from 10 to 80 m/day, and µ ranges from 0.1 to 0.3.The final results in every zone are shown in Table 1.K is used to identify the boundary conditions.

Recharge and Discharge
Precipitation and withdrawal are the main recharge and discharge, respectively, in JUA.Irrigation is prominent from May to October, the growing season, in every year.The quantity of irrigation is based on the Jilin province Local Standard of Water Quota (DB22/T 389-2010) [37].Irrigation water is derived from river and groundwater pumping in the irrigation district, but wells in the district are not monitored, so area features were used to simulate the withdrawal of irrigation water with the recharge (RCH) package [38].The direct evaporation and evapotranspiration (ET) package are based on the data of the four meteorological stations.The precipitation (positive recharge rate) is simulated with the RCH package as well as the irrigation (positive recharge rate), and irrigation from groundwater pumping (negative discharge rate); however, the recharge package does not allow for recharge and discharge to occur simultaneously in the same vertical column [10].Therefore, precipitation and irrigation were superimposed using Microsoft Excel and assigned to the area as a transient recharge rate to the first layer.In the ET package, the parameters of elevation, ET extinction depth, and maximum ET rate are needed.The evaporation extinction depth was 5 m based on reports and literature; the evapotranspiration occurs at the maximum ET rate until the water table rises above the elevation.If the water table elevation lies between these two extremes, the evapotranspiration rate varies linearly with depth [39].Precipitation, direct evaporation, and irrigation data were collected from Jilin, Dakouqin, Soudenghe, and Huapichang weather stations and the statistical yearbook of the Jilin province.Zone polygons were performed using these four stations by combining the land use, lithology, and the irrigation district, which resulted in 50 recharge zones being defined in the RCH Package and 45 evapotranspiration zones in the ET Package.The zones are shown in Figure 8.
There is a large amount of centralized industrial water in the industrial zone, and the pumping wells are near the second Songhua River; the locations of the industrial pumping well are shown in Figure 5 and the pumping rate is variable with the time.Withdrawal is simulated by the well (WEL) package as point features based on monitoring data.The pumping discharge was based on the 2006 monitoring data.The quantity of irrigation water was large in the agricultural district in the North.

Interaction between Groundwater and Rivers
JUA has an asymmetrical topography.The groundwater level changes substantially because of the varied topography and the crisscross drainage system.Groundwater levels obtained from boreholes were interpolated by Kriging and integrated into the model.There is a good connection between rivers and the aquifer in JUA.Rivers are recharged by groundwater during the dry season, but groundwater is recharged by rivers during the wet season.Rivers are defined as lines in the river package in MODFLOW.Elevation of river water levels and bottoms of each point along the river channel are needed.From Figure 9a we can see that the groundwater level, incorrectly, crosses the river twice in some places.The values of the river stage along the river channel should be added to control the groundwater flow direction when interpolating the groundwater initial head with the existing borehole.There are six steps for interpolating the river stage.
(1) The 1:50,000 topographic map is imported and registered into GMS and calibrated by using the coordinates in GMS; (2) Then, the origin of a river can be defined as (X1, Y1); (3) Each coordinate along the river can be read as (Xi, Yi); (4) And the length (Li) from each point (Xi, Yi) along the river to the original point on the river (X1,Y1) can be calculated by using the formula below: (5) According to the river stage data (Zi) which are monitored by several hydrometric stations within the Jilin province, the scatters of the river stage (Zi) versus the length (Li) from the river stage to the original point on the river can be obtained (Zi, Li); (6) From the trend line of the scatters, we can determine the river stage anywhere along the rivers.
(the Songhua River and Mangniu River stage interpolations are shown in Figure 10.
After adding the river stage, the groundwater flow field is shown in Figure 9b.

Model Calibration and Sensitivity Analysis
It is important to decide which result is reasonable in the calibration process.In this procedure, two types of data were used to calibrate the regional groundwater flow model for the JUA.First, observed groundwater level contour maps of the unconfined aquifer derived from the 190 borehole data points and 75 observation wells were used to match the computed water levels (Figure 11) as a qualitative step before more comprehensive and quantitative calibration was performed [32].Second, the variable heads of the 75 observation wells from 2001 to 2010 were used as calibration targets (well locations as shown in Figure 5.The observation wells are distributed throughout the study area and, thus, are representative of the regional aquifer conditions.Trial and error was used to calibrate the model using the following procedure [31].First, the recharge rate was adjusted to minimize the difference between the observed and calculated groundwater level contour maps and heads.When this was achieved, the hydraulic conductivity and the specific yield were adjusted to minimize the variance of the differences [40].The final hydraulic conductivity and the specific yield are shown in Table 1.This procedure was repeated until the variance was in a reasonable range. It can be seen from Figure 11 that the contour maps of the observed groundwater level and the calculated groundwater level agreed well.Variance existed in the nearby Mangniu River and the upstream of the second Songhua River.The following reasons may lead to this result: (1) there are few monitoring wells for constructing an accurate water level contour map; (2) there are many pumping wells (Figure 5) that influence the water level, but the pumping wells are located close to the rivers, and while there is a good connection between the groundwater and the river, these problems require further study.
The variable observed hydraulic head and the calculated head were compared with the scatter shown in Figure 12.It can be seen that there was a good match between the observed and calculated head values.The mean error was 0.547 m and the root mean squared error (RMSE) was 0.66 m.The difference of the highest and the lowest water level was over 60 m so the discrepancy was reasonable.Sensitivity analysis was performed on parameters such as the horizontal and vertical hydraulic conductivity, specific yield, recharge, and river conductance.In a sensitivity analysis, the input parameter values are varied across the range of likely values [20].The parameter is increased or decreased 10% or 20% while other parameters remain unchanged.The model runs once for each change in parameters, and the effect upon the computed head is noted.The model is responsive to changes in recharge and horizontal hydraulic conductivity (Kh), river conductance (Cond.),vertical hydraulic conductivity (Kv), and specific yield (sy) are relatively small.The results are shown in Table 2.

Water Budget
It is a difficult task to calculate the water balance in an intermountain river valley model that has complex boundary conditions, and complex interaction between surface water and groundwater.
The water balance equation can be defined as: where Qp is the infiltration of precipitation, Qrl is the river leakage, Qirr is the infiltration of the irrigation in irrigation districts, Qrb is the recharge from the lateral boundaries, QE is the ET, Qr is the discharge from groundwater to a river, Qdb is the discharge from the lateral boundaries, Qw is the artificial pumping (large quantity wells and pumping water for irrigation), and ΔQ is the storage change of the groundwater system.
The "water budget" in the tools menu of GMS provides an effective tool for analyzing the water budget for the JUA.The repeated recharge from irrigation and precipitation should be subtracted to avoid adding the recharge rate as the irrigation does not occur during a rain spell.Precipitation is based on monitoring data from the four meteorological stations from 2001 to 2010 and is given the same as the step of the model.The average annual precipitation of the 10 years is 1.40 × 10 8 m 3 /a, making the total inflow 67.68%.The inflow from the lateral boundaries, the infiltration of the irrigation and the river leakage are the annual average recharge to the groundwater system.The boundary inflow from the mountains is 0.17× 10 8 m 3 /a, which makes the total inflow 8.27%.The irrigation is 0.05 × 10 8 m 3 /a and the river leakage to groundwater is 0.45 × 10 8 m 3 /a, which makes the total inflow 21.73%.
The withdrawal is based on the monitoring data of 2006, with a quantity of 1.57 × 10 8 m 3 /a and is the main discharge of the groundwater system, which makes the total outflow 71%.ET from the natural groundwater flow system is 0.10 × 10 8 m 3 /a.ET should be further assessed because irrigation may increase the value of ET.Groundwater discharge to the river occurs when the groundwater level is higher than the river stage and was as much as 0.41 × 10 8 m 3 /a in the dry season, while outflow from the lateral boundaries was 0.13 × 10 8 m 3 /a.The water balance of the model calculated with the JUA flow budget module shows that recharge from precipitation and discharge from withdrawal were the main terms of the water balance and had major influences on the groundwater levels.
In the water budget, the input was 2.07 × 10 8 m 3 /a, which was smaller than the output of 2.21 × 10 8 m 3 /a.From Figure 11, we see that, in the irrigation district and the industrial district near the big pumping wells, the cone of the depression and the groundwater level decline.

The Parameters of Each Sub-Region
The study area was divided into three hydrogeological zones with 37 sub-regions.For the top aquifer, the lithology of the Songhua River valley is silt and the K is defined as 1-2.7 m/day.The μ value is 0.07-0.15,which is larger than the other locations where the lithology consists of clay and loam and the K is 0.3-0.96m/day and μ is 0.05-0.08.For the lower sand and gravel aquifer, the calculated K ranges from 10 to 80 m/day, and µ ranges from 0.1 to 0.3.The K for the Songhua River valley is larger than that of the Aolong River valley and larger than the platform areas.The parameters are based on the average pumping test results in each sub-region, and the parameter results are within the expected empirical values and are similar with Qi's result [41] in I1-6 where the K is 33.6 m/day, μ is 0.24, and Shi's [42] result in the Songhua river valley where the K is 25-80 m/day.

Model Limitations
Numerical modeling can be used to simulate the real groundwater flow system and is based on the literature, borehole data and the experience of the modeler, but limitations of the technology and the data exist: (1) Limited borehole and groundwater level monitoring data, especially in the mountains and the outcrops of the underlying basement and where the thickness of the aquifer changes greatly, affect the reliability of the model.(2) Aquifer and irrigation recharge rates and spatial distribution of aquifer recharge areas.
The irrigation recharge rates are estimated from the literature because there are no monitoring wells in this area, which results in some uncertainty.(3) Withdrawal was based on the 2006 data, which may overestimate reality prior to 2006 and underestimate that following 2006, as the plant area increased and there was growth in industry and increased population.Further investigation is required to make a more accurate model for the JUA.(4) There are 190 borehole data points and 75 observation wells in this model, but the initial head is mainly based on the observation wells because many of the boreholes were developed in the 1980s and the water level has already changed greatly, which introduces another source of uncertainty.(5) Infiltration of irrigation is based on the average value of water-use quota of the Jilin province.This rate is not accurate because different crops need different quantities of water; therefore, future studies should investigate land types and crop species.

Conclusions
A groundwater flow model was developed using GMS to help understand a complex groundwater system in in Jilin urban area of Jilin province, China.Based on the geological and hydrogeological investigations, the aquifer system here was divided into two aquifers that mainly composed of silt, loam, sand and gravel, respectively.Both aquifers are commonly unconfined, but are confined at those places overlain by the loam.
The elevation of land surface was inferred from the DEM data and the 1:50,000 topographic map.The initial water levels were interpolated by Kriging method based on the water level measurements.
The water stage was interpolated by using the observation data at several hydrometric stations and the coordinate system in GMS, and then it was combined with the observation wells to get groundwater level accuracy where it crosses the river.
The trial-and-error method was used in calibration of the model.In this method, expert judgment plays a more important role than automatic procedures where calibration is restricted to the calibration zones with a limited amount of variables [19].Sensitivity analysis was used to explain the role of calibration parameters in the groundwater response, and the result showed that recharge was the most sensitive factor.
Based on the supply and demand analysis of water resources, the established model would finally provide a scientific basis to use the groundwater resources sustainably in JUA.The groundwater flow model can be used as the basis for a water quality model and the water flow and quality prediction model of JUA.This is the first step to studying the JUA; however, further detailed investigations are needed to accurately model JUA and provide a basis for groundwater management in Jilin city.

Figure 2 .
Figure 2. Procedure for data processing and groundwater modeling studies.

Figure 3 .
Figure 3. Annual precipitation and evaporation of the study area for the last 30 years.

Figure 5 .
Figure 5. Location of the boreholes, observation and pumping wells and zones.

Figure 7 .
Figure 7. Geometrical model of JUA.(a) Geometrical model solid and elevation; (b) Grid, boundary conditions and pumping wells locations of the model.

Figure 8 .
Figure 8. Zones for Recharge (RCH) and evapotranspiration (ET) packages, where different colors represent different zones.(a) Zones are superimposed by precipitation and irrigation for the RCH Package; (b) zones of evaporation for the ET Package.

Figure 9 .
Figure 9. Groundwater level in JUA (a) interpolated just using the hydrogeological wells and (b) interpolated using the hydrogeological wells and the river stage.

Figure 11 .
Figure 11.Comparison of observed and computed groundwater contour level.

Figure 12 .
Figure 12.Scatter plot of observed versus computed values of head for model calibration.

Table 1 .
Final parameters results for each zone.

Table 2 .
Groundwater head changes versus the changes of different parameters.