The Spatial and Temporal Variability of Groundwater Vulnerability and Human Health Risk in the Limin District , Harbin , China

This study aimed to analyze the variations in groundwater quality, vulnerability and potential health risk from 2006 to 2016 in the Limin District, Harbin, China. Groundwater geochemical characteristics were described using statistical analysis and Piper diagrams. A modified DRASTIC model that combined factors of intrinsic aquifer vulnerability and land use was applied to assess groundwater vulnerability. The weights of parameters were adjusted by using the analytic hierarchy process (AHP) to optimize the model. The non-carcinogenic health risk was estimated by the Unites States Environmental Protection Agency (USEPA) model. Results suggested that concentrations of NH4-N, Fe and Mn in groundwater exceeded the limits both in 2006 and in 2016. The concentration of Fe in the groundwater showed more significant variation between 2006 and 2016 than the other parameters. Very high vulnerability zones increased from 6.3% in 2006 to 16.9% in 2016, and distributed on agricultural land, indicating that agriculture was still a major source of pollution. Mn and NO3-N contributed the most to human health risks in 2006 and 2016, respectively. This study highlights the influence of groundwater quality variation in decadal exploitation on human health.


Introduction
Water resources are one of the primary factors that determine the degree of regional economic development.In turn, the development of the region will also affect local water quality and quantity [1].In many developing countries, where agriculture is a fundamental part of the economy [2], rapid industrialization is resulting in greater pressure on water resources [3,4].Groundwater has been an important source of water supply because of its abundance, stable quality, and the relatively low cost of exploitation [5][6][7].
As water consumption continues to increase, the release of wastewater from factories, agriculture, and domestic sources increases [8,9].As a result of these problems, groundwater resource depletion and quality degradation are becoming more obvious.Groundwater vulnerability and human health risk assessment are useful tools for groundwater exploitation and management.There are many methods of groundwater vulnerability assessment, including: statistical, process-based and index-overlay including DRASTIC, GOD, SI, etc. [6].The DRASTIC model [10] is the most established method worldwide for presenting the intrinsic vulnerability of groundwater due to its minimal data requirements, relatively accurate results and flexible application [11][12][13].Therefore, various studies have modified the original model in consideration of area-specific ratings and weights [14][15][16][17][18][19][20][21][22].For example, using the analytic hierarchy process (AHP) to optimize the rates and weights of the DRASTIC model makes it possible to obtain better results [13,23].In addition, establishing the link between intrinsic groundwater vulnerability and land use can give a clear view of the impact of human activities on the environment [24][25][26].
Groundwater is a valuable freshwater source, and its physical, chemical, and biological parameters are of serious concern when being used for drinking purposes [27].Human health risk assessment can reflect the safety of groundwater.The model for human health risk assessment of USEPA has been widely used in China.China's Ministry of Environmental Protection proposed some new technical guidelines [28,29] to put forward the study of domestic environmental health risk assessment.The Geographical Information System (GIS) provides a powerful tool for mapping and displaying spatial information on contaminants and health risks, and can offer more specific advice to government managers on groundwater prevention and treatment.
This study was carried out in the Limin District, which is in the northern part of Harbin city in Heilongjiang province, northeast China.The Limin District is an economic development zone with fast growth, and is facing water resource management problems.With regard to long-term exploitation of groundwater, issues such as what happens to the groundwater environment and how these changes affect human health have received less attention.The main objectives of our previous studies were to : (1) compare the temporal and spatial variation of Limin District's groundwater geochemical characteristics and vulnerability and discuss the reasons for these; and (2) analyze the impact of groundwater exploitation in the long-term on environment and human health.The modification of the DRASTIC model was done by replacing net recharge and topography with river distance and land use.This study integrated the human health and groundwater vulnerability associated with the modified DRASTIC model to reflect the relationship between groundwater environment and human activities.

Study Area
Limin District is at the confluence of the perennial Songhua and Hulan Rivers (Figure 1a).It covers an area of about 658 km 2 and lies between 45 • 41 50 and 46 • 03 30 N lat and 126 • 20 16 ~126 • 48 30 E long.The landscape of Limin District is flat, mostly consisting of alluvial plain, except for a small part of high plain in the northwest corner.The area lies within the temperate continental monsoon, with an average annual temperature and precipitation of, respectively, 3.6 • C and 547.5 mm, of which 70~80% occurs in summer (July~September).Cretaceous Nenjiang Formation (K 1n ) and Quaternary sediments (Q 1 and Q 4 ) are generally located in the area (Figure 1b).Because of basement fluctuation, Quaternary sediments have varying thicknesses ranging from 24 to 53 m.Silty clays comprise the upper part of the Quaternary sequence, while the lower parts are sands and gravels.The porous aquifer developed in the Quaternary sediments is main aquifer and can be divided into high floodplain semi-confined-phreatic aquifer and low floodplain phreatic aquifer.Both of them are characterized by highly permeable sand and gravel layers interlayered with about 20~40 m thickness totally.The regional groundwater piezometric head decreases mainly in a west to east direction and according to pumping data, the aquifer is sorted into three water abundance divisions: very rich water area (prediction of single well water inflow above 5000 m 3 /day), rich water area (prediction of single well water inflow at 3000~5000 m 3 /day), and middle water area (prediction of single well water inflow at 1000~3000 m 3 /day) [30].Groundwater from Quaternary sediments is more significant in the local public water supply system.In 2006, the comprehensive exploitation of groundwater was 1.3 million cubic meters per year, which was mainly used for domestic, irrigation and industrial purposes.For the year 2016, water demand reached about 6.1 million cubic meters per year.Additionally, although the study area is dominated by agriculture, industry is developing rapidly.

Groundwater Sampling and Analyses
During the wet seasons of July 2006 and 2016, 13 and 12 groundwater samples (Figure 1), respectively, were collected in the shallow (depth generally < 30 m) Quaternary aquifer.Some samples were collected from private hand-dug wells and others from boreholes.Although samples from 2006 and 2016 were not always at the same location due to the discontinuation of some wells, nearby sample points were chosen.The water samples were filtrated through a 0.45 μm filter membrane and collected in clean 500 mL plastic bottles.All samples were stored in airtight ice-cold containers and transported to the laboratory as soon as possible.Total dissolved solids (TDS), pH and total hardness (TH) were measured in the field using multi-function water quality tester HI7609829.Others, including major ions (HCO3 − , SO4 − , Cl − , K + , Na + , Ca 2+ , Mg 2+ ), nitrate (NO3-N), nitrite (NO2-N), ammonia nitrogen (NH4-N), manganese chemical oxygen demand (CODMn), and metals (Fe, Mn), were tested and analyzed in the laboratory according to the technical specifications for the environmental monitoring of groundwater in the national professional standards of China [31].

Statistical and Spatial Analyses
Statistical information of the main physical and chemical components including the minimum, maximum, mean, and standard deviation of mean were calculated with Microsoft Excel 2013.Waterchemistry characteristics were presented by using Aq•QA version 1.1 (RockWare, Colorado, USA).A Kruskal-Wallis test as a non-parameter test (p < 0.05) was evaluated by SPSS 20.0 (IBM Corporation, New York, USA).ArcGIS 10 software (ESRI, Toronto, Canada) was used to display spatial structures variation of concentrations, groundwater vulnerability and human health risks with time.Groundwater from Quaternary sediments is more significant in the local public water supply system.In 2006, the comprehensive exploitation of groundwater was 1.3 million cubic meters per year, which was mainly used for domestic, irrigation and industrial purposes.For the year 2016, water demand reached about 6.1 million cubic meters per year.Additionally, although the study area is dominated by agriculture, industry is developing rapidly.

Groundwater Sampling and Analyses
During the wet seasons of July 2006 and 2016, 13 and 12 groundwater samples (Figure 1), respectively, were collected in the shallow (depth generally < 30 m) Quaternary aquifer.Some samples were collected from private hand-dug wells and others from boreholes.Although samples from 2006 and 2016 were not always at the same location due to the discontinuation of some wells, nearby sample points were chosen.The water samples were filtrated through a 0.45 µm filter membrane and collected in clean 500 mL plastic bottles.All samples were stored in airtight ice-cold containers and transported to the laboratory as soon as possible.Total dissolved solids (TDS), pH and total hardness (TH) were measured in the field using multi-function water quality tester HI7609829.Others, including major ions (HCO 3 − , SO 4 2− , Cl − , K + , Na + , Ca 2+ , Mg 2+ ), nitrate (NO 3 -N), nitrite (NO 2 -N), ammonia nitrogen (NH 4 -N), manganese chemical oxygen demand (COD Mn ), and metals (Fe, Mn), were tested and analyzed in the laboratory according to the technical specifications for the environmental monitoring of groundwater in the national professional standards of China [31].

Statistical and Spatial Analyses
Statistical information of the main physical and chemical components including the minimum, maximum, mean, and standard deviation of mean were calculated with Microsoft Excel 2013.Water-chemistry characteristics were presented by using Aq•QA version 1.1 (RockWare, Colorado, USA).A Kruskal-Wallis test as a non-parameter test (p < 0.05) was evaluated by SPSS 20.0 (IBM Corporation, New York, NY, USA).ArcGIS 10 software (ESRI, Toronto, Canada) was used to display spatial structures variation of concentrations, groundwater vulnerability and human health risks with time.

Groundwater Vulnerability Assessment Using Modified DRASTIC Model
The DRASTIC model is a stacking exponential method that was proposed by USEPA to protect groundwater resources.The model is based on based on the following seven hydrogeological parameters to evaluate intrinsic vulnerability of aquifer system: Depth to groundwater (D), Net Recharge (R), Aquifer media (A), Soil media (S), Topography (T), Impact of vadose zone (I), and Hydraulic Conductivity (C).The DRASTIC Index (DI) is calculated by Equation (1), below: where r is the rating of parameters from 1 to 10 decided to their relative effect on the groundwater vulnerability, and w is the weight of each parameter that reflects its importance to vulnerability and is computed by analytic hierarchy process (AHP).
The DRASTIC model is mature, and has been applied widely.Combined with actual conditions of different study areas, many researchers have modified the model (Table 1).Changes in rainfall and topography at Limin District are negligible according to meteorological and topographic data in Harbin; therefore, in this study, the model was revised as follows: (1) qualitative parameter of aquifer media was replaced with quantitative parameter of aquifer thickness (A); (2) net recharge and topography were omitted; (3) the factors of river distance (R) and type of land use (L) were added.The ratings are illustrated in Table 2 and the modified model-DRASICL-was applied in the present study.Then, by using the spatial analysis technology of ArcGIS Software, vulnerability mapping of groundwater was generated.Land use type should be considered to be one of the groundwater vulnerability indexes in order to achieve more objective vulnerability.

Huan et al. [21]
DRASTICL-Lin model AHP The classes of vulnerability were defined more clearly with the use of weight values determined by the AHP method.
Sener et al. [33] DRASTIL model Referring Aller et al. [10] Land use was found to be more effective parameters in assessing groundwater vulnerability than assumed by original model.

Sinha et al. [34]
DRTILPQ model AHP Land use was closely related to groundwater quality and had been adopted as an index in many vulnerability models.

Wu et al. [23]
Rd-distance to river network; N-average pollution load of NO 3 -N; V-groundwater velocity; Lin-Lineaments; P-impact of pollution sources; Q-water quality.Agricultural land 9

Analytic Hierarchy Process
The analytic hierarchy process (AHP), a multi-criteria decision making (MCDM) technique, was proposed by Saaty [35].AHP has been a good method for solving complicated decision-making problems [13,33,36,37].In this study, AHP was used to modify the weights of seven parameters based on the local hydrogeological conditions.The weights of parameters reflect their relative importance on the result of vulnerability by a series of pairwise comparison judgments.The judgments were determined depending on the weights of parameters in the initial DRASTIC and references [21,23,33].There was a slight difference in the weights between 2006 and 2016, considering that the contribution of some factors changed with time.For example, land use type is increasingly important to groundwater vulnerability assessment, and had a higher value of weight compared to depth of groundwater in 2016.The pairwise comparison matrix is carried out using the following formulas: w i is the weight of parameter i, and n is the number of parameters.a ij is an element of the pairwise comparison matrix and the more important parameter i than j, the bigger the a ij .
In addition, the quality of the result of the AHP is significantly related to the consistency of the pairwise comparison judgments.In this approach, the consistency of the matrix requires the consistency ratio CR ≤ 0.1.The formula for the CR is as follows: where CI is the consistency, RI stands for random index, and the value of RI is 1.36 (n = 7).

Human Health Risk Assessment
The potential non-cancer risk for individual contaminants in the dose-response assessment can be characterized with a hazard quotient (HQ) [38].Hazard index (HI) is used to assess the overall potential non-cancer risks posed by more than one contaminant.Non-carcinogenic risk assessment through drinking water intake is reckoned as follows: where Intake oral-water is the daily average exposure dosage by drinking water intake per unit weight (mg kg −1 day −1 ).i refers to some pollutant in groundwater and C i is concentration of pollutant (mg L −1 ).IR represents the ingestion rate of water, and EF is exposure frequency.ED indicates exposure duration and BW is the average body weight.In this study, AT denotes the average time of non-carcinogenic effects on adults and children.HQ i and RfD express the hazard quotient and the reference dosage, respectively, for non-carcinogenic contaminants through oral intake pathway.Values of the above exposure parameters and RfD are shown in Table 3.In this study, NO 3 -N, NH 4 -N, Mn, and Fe were selected as the risk assessment parameters, because most of them are primary contaminants in the area.Non-carcinogenic risk assessment was applied, as these parameters are non-carcinogenic pollutants according to the International Agency for Research on Cancer (IARC) and USEPA.Since oral intake by drinking water is the most common exposure pathways for contaminated groundwater [40,41], oral intake was taken into account in the present study.

Groundwater Geochemical Characteristics
Statistical information of groundwater geochemical parameters, including the minimum, maximum, mean, and standard deviation, is presented for 2006 and 2016 in Table 4.The mean value of pH was 7.28 in 2006 and 7.71 in 2016.TH, TDS, and COD Mn were relatively stable in 2006 and 2016.All of them were in the acceptable levels of drinking water quality set by the Ministry of Health of the P.R.China (PRC) and the World Health Organization (WHO) (Table 4).Table 4. Statistical summary of geochemical parameters and their comparison with the limits of PRC drinking water standards [42] and WHO guidelines values [43] (units = mg/L).Table 4 shows that the average concentrations of NH 4 -N, TFe and Mn exceeded the acceptable levels both in 2006 and 2016.Since in more than half of the groundwater samples, NO 2 -N was not detected, we decided not to take it into account in this study.Figure 3 shows that the concentration of NO 3 -N had a decreasing trend from 2006 to 2016.In 2006, NO 3 -N concentrations were 2~10 mg/L, while in 2016 it was less than 5 mg/L in most areas.High NO 3 -N concentration in groundwater occurs frequently due to excessive irrigation water and fertilizer [44].Nevertheless, there was only one sample in 2006 with a NO 3 -N concentration exceeding (>10 mg/L) the water quality standard.Reduction of agricultural area (Figure 4) and restrictions on the use of chemical fertilizers [45] contributed to these low concentrations.Denitrification has a key role in nitrate attenuation in shallow aquifers, because it greatly reduces the nitrate concentrations in water [5,46,47].We suspected that the fresh water recharging mainly from precipitation had some dilution effect on the groundwater concentration of NO 3 -N, and this was consistent with the findings of Zhai et al. [48].The concentrations of NH 4 -N in large areas were 0.5~1.0mg/L in 2006 and it increased to 1.0~1.5 mg/L in 2016.Release of wastewater from factories, agricultural runoff, and domestic sewage constituted both point and non-point sources of NH 4 -N [49].The rapid development of the study area contributed to the increasing of NH 4 -N pollution from 2006 to 2016.There are many sources of pollution distributed along the Songhua River and Hulan River, including both point and non-point sources, resulting in large emissions of NH 4 -N.These sources have lowered the water quality to a point where it has reached the surface water quality standards IV grade (>1 mg/L) of the PRC [50][51][52].This may increase the risk of groundwater pollution with the decreasing river water quality [53,54].It has been reported that high contents of Fe and Mn in groundwater are common in northeastern China, with high Fe concentrations generally occurring where Mn concentrations are high [55][56][57].Especially in the center of the Songnen Plain, in which our study area occurs, contents of Fe and Mn are as high as 11.2-44.4mg/g and 0.5-1.0mg/g, respectively, due to the abundance of iron-manganese nodules in local Quaternary aquifer sediments [55,58].Therefore, the background concentrations of Fe and Mn in the study area are very high and usually exceed the acceptable levels.Groundwater Fe concentration showed a considerably decrease from 2006 (Figure 3e) to 2016 (Figure 3f) and this change corresponded with the results of Sharma, et al. [59] and the Kruskal-Wallis test (Table 4).The Mn concentration decreased insignificantly and varied from 0.57 mg/L in 2006 to 0.51   It has been reported that high contents of Fe and Mn in groundwater are common in northeastern China, with high Fe concentrations generally occurring where Mn concentrations are high [55][56][57].Especially in the center of the Songnen Plain, in which our study area occurs, contents of Fe and Mn are as high as 11.2-44.4mg/g and 0.5-1.0mg/g, respectively, due to the abundance of iron-manganese nodules in local Quaternary aquifer sediments [55,58].Therefore, the background concentrations of Fe and Mn in the study area are very high and usually exceed the acceptable levels.Groundwater Fe concentration showed a considerably decrease from 2006 (Figure 3e) to 2016 (Figure 3f) and this change corresponded with the results of Sharma, et al. [59] and the Kruskal-Wallis test (Table 4).The Mn concentration decreased insignificantly and varied from 0.57 mg/L in 2006 to 0.51 It has been reported that high contents of Fe and Mn in groundwater are common in northeastern China, with high Fe concentrations generally occurring where Mn concentrations are high [55][56][57].Especially in the center of the Songnen Plain, in which our study area occurs, contents of Fe and Mn are as high as 11.2-44.4mg/g and 0.5-1.0mg/g, respectively, due to the abundance of iron-manganese nodules in local Quaternary aquifer sediments [55,58].Therefore, the background concentrations of Fe and Mn in the study area are very high and usually exceed the acceptable levels.Groundwater Fe concentration showed a considerably decrease from 2006 (Figure 3e) to 2016 (Figure 3f) and this change corresponded with the results of Sharma, et al. [59] and the Kruskal-Wallis test (Table 4).The Mn concentration decreased insignificantly and varied from 0.57 mg/L in 2006 to 0.51 mg/L in 2016.For the year 2016, areas of higher concentrations of Mn exhibited a decrease compared with 2006.Dissolution processes of minerals, which are often controlled by the redox level of groundwater, affect reductive dissolution of Fe/Mn oxides.Increasing the exploitation of groundwater will cause water levels to drop, leading to changes in the redox conditions as water table pumping delivers more O 2 to the subsurface [4,60].Dissolved Fe from the reducing zone may precipitate as in the oxidation condition [4].

Pairwise Comparison Matrix Used in AHP
The values of judgment and weights calculated are shown in Tables 2 and 5, respectively.The consistency ratio CR 2006 = 0.070 and CR 2016 = 0.052 were both less than 0.1, which indicated that the judgment matrix passed the consistency check.Note: the values below the diagonal are the proportion of weight of the criteria regarding their relative significance in the assessments of groundwater vulnerability issue in 2006, while the values above the diagonal are the proportion of weights of the criteria in 2016.

Groundwater Vulnerability Assessment
The assessment of groundwater vulnerability in this study was conducted based on modified DRASICL model.Spatial maps of seven parameters are shown in Figure 5. Variation of water table depth and land use were compared in our study based on the assumption of other parameters had no change over time.During the study period, water table depth in both 2006 and 2016 were less than 10 m, which indicated the groundwater to be more susceptible to contamination, based on Table 2. Most of the area in 2006 had a groundwater depth range of 2 to 5 m, which was associated with rating score of 8 (Figure 5c).Areas with deeper depths in 2016 had increased, and were mainly were found in the central part.Although consumption of groundwater in 2016 increased significantly, recharge of groundwater during the rainy season reduced the drop of water table [44,61].
Free maps of nationwide land use type are provided every five years at the website of Geospatial Data Cloud, so maps of land use in 2005 and 2015 were evaluated.It can be seen from Figure 4 that agricultural land occupied the largest area within the Limin District.In 2015, the areas of unused land, woodland, and grassland changed slightly compared with 2005.The drastic changes occurred at construction land that ranged from 8.0% in 2005 to 14.3% in 2015.It seems that an increase of construction land led to a decrease in agricultural land according to the comparison with Figure 5a,b.To a certain degree, types of land use have determined the variety and quantity of pollutants [6].Irrational use of land will directly result in groundwater pollution, thereby affecting future land planning.With population growth and industrial development, area of construction land nearly doubled between 2005 and 2015.
The impact of rivers water on groundwater in the study area is obvious as varying with distance from rivers where interactions between groundwater and surface water occur.Anthropogenic sources such as industrial sources, fecal pollution, livestock wastewater, and agricultural pollution increase pollution of the rivers, which intensify groundwater vulnerability [32], especially during the wet season, when groundwater is recharged from the rivers.Song et al. [54] analyzed the distance between the groundwater sampling points and Songhua riverbank and inferred that the K + , NH 4 -N, and Cl − concentrations of groundwater increased owing to the recharge of polluted surface water.Figure 5e shows the distance from the rivers.Areas with closer distance from the rivers were associated with higher rating values that indicated stronger hydraulic connection between groundwater and the rivers.Both the level and quality of the river have impacts on groundwater system.and Cl − concentrations of groundwater increased owing to the recharge of polluted surface water.Figure 5e shows the distance from the rivers.Areas with closer distance from the rivers were associated with higher rating values that indicated stronger hydraulic connection between groundwater and the rivers.Both the level and quality of the river have impacts on groundwater system.In general, aquifer thickness in the study area was relatively homogeneous (Figure 5f), and thickness, which ranged from 30 to 40 m, occupied most of the area.Different soil media can reflect their infiltration capacity.Lower rating values are associated with lower infiltration capacity, which provides lower risk of groundwater pollution.Figure 5h shows that eastern and northeastern parts of the vadose zone are occupied by silt and silt sand, which are assigned lower rating values.There are mainly two ranges of hydraulic conductivity and are assigned rating scores of 7 and 9. Hydraulic conductivity in the eastern part of study area ranged between 30 and 50 m/day, corresponding to a higher rating score of 7 (Figure 5i).
The groundwater vulnerability maps of 2006 and 2016 and areas of different classes of vulnerability are shown, respectively, in Figure 6 and in Table 6.Very low-and low-vulnerability zones in 2006 and 2016 accounted for 21.6% and 23.9% of the total area, respectively.Low vulnerability areas were mainly distributed in the low floodplain.The largest areas in 2006 and 2016 were determined to have high vulnerability and were, respectively, 47.5% and 37.0%.A very highvulnerability area, which was found near the northwest corner, increased by 10% in 2016, compared with 2006.High-and very high-vulnerability zones in 2006 and 2016 corresponded to cultivated land, and this suggests that agriculture contributed more pollutants to shallow groundwater aquifers.In general, aquifer thickness in the study area was relatively homogeneous (Figure 5f), and thickness, which ranged from 30 to 40 m, occupied most of the area.Different soil media can reflect their infiltration capacity.Lower rating values are associated with lower infiltration capacity, which provides lower risk of groundwater pollution.Figure 5h shows that eastern and northeastern parts of the vadose zone are occupied by silt and silt sand, which are assigned lower rating values.There are mainly two ranges of hydraulic conductivity and are assigned rating scores of 7 and 9. Hydraulic conductivity in the eastern part of study area ranged between 30 and 50 m/day, corresponding to a higher rating score of 7 (Figure 5i).In general, aquifer thickness in the study area was relatively homogeneous (Figure 5f), and thickness, which ranged from 30 to 40 m, occupied most of the area.Different soil media can reflect their infiltration capacity.Lower rating values are associated with lower infiltration capacity, which provides lower risk of groundwater pollution.Figure 5h shows that eastern and northeastern parts of the vadose zone are occupied by silt and silt sand, which are assigned lower rating values.There are mainly two ranges of hydraulic conductivity and are assigned rating scores of 7 and 9. Hydraulic conductivity in the eastern part of study area ranged between 30 and 50 m/day, corresponding to a higher rating score of 7 (Figure 5i).
The groundwater vulnerability maps of 2006 and 2016 and areas of different classes of vulnerability are shown, respectively, in Figure 6 and in Table 6.Very low-and low-vulnerability zones in 2006 and 2016 accounted for 21.6% and 23.9% of the total area, respectively.Low vulnerability areas were mainly distributed in the low floodplain.The largest areas in 2006 and 2016 were determined to have high vulnerability and were, respectively, 47.5% and 37.0%.A very highvulnerability area, which was found near the northwest corner, increased by 10% in 2016, compared with 2006.High-and very high-vulnerability zones in 2006 and 2016 corresponded to cultivated land, and this suggests that agriculture contributed more pollutants to shallow groundwater aquifers.

Human Health Risk Assessment
Higher concentrations of TFe, Mn, NH 4 -N, and NO 3 -N in groundwater are more likely to pose high risks to human health.They were used to conduct the non-carcinogenic health risk assessment, and the assessment results for adults and children in both years are shown in Table 7.We found that the average HQ Mn for adults and children was higher than other constituents' HQ in 2006, while HQ NO3-N was the highest in 2016.The non-carcinogenic risk values of TFe and Mn in 2016 for adults and children were markedly lower than in 2006.Recent research has shown that concentrations of Mn exceeding the drinking water standard may do harm to children's intelligence [62,63].Although the concentration of Mn is high in the Limin District, it is not difficult to remove it simply by using treatments such as aeration and sand filtration [64].There were slightly higher health risks of NH 4 -N and NO 3 -N for adults and children in 2016.High values of NO 3 -N may pose higher risk of human illness, such as methemoglobinemia, gastric cancer, goiter, hypertension, etc. [65].Spatial distribution of HI for adults (Figure 7a,b) show that most areas have little to no health risks except for a small part of the northwest in 2016.HI for children was generally higher than for adults (Table 7), whether in 2006 (Figure 7c) or 2016 (Figure 7d).In 2006, 91.8% of the study area had HI between 1.0 and 2.0, while this figure decreased to 34.4% in 2016.Additionally, in 2016, 58.4% had HI < 1.0, and 7.3% located at the northwestern corner had HI > 2.0.
Water 2018, 10, x FOR PEER REVIEW 13 of 18

Human Health Risk Assessment
Higher concentrations of TFe, Mn, NH4-N, and NO3-N in groundwater are more likely to pose high risks to human health.They were used to conduct the non-carcinogenic health risk assessment, and the assessment results for adults and children in both years are shown in Table 7.We found that the average HQMn for adults and children was higher than other constituents' HQ in 2006, while HQNO3-N was the highest in 2016.The non-carcinogenic risk values of TFe and Mn in 2016 for adults and children were markedly lower than in 2006.Recent research has shown that concentrations of Mn exceeding the drinking water standard may do harm to children's intelligence [62,63].Although the concentration of Mn is high in the Limin District, it is not difficult to remove it simply by using treatments such as aeration and sand filtration [64].There were slightly higher health risks of NH4-N and NO3-N for adults and children in 2016.High values of NO3-N may pose higher risk of human illness, such as methemoglobinemia, gastric cancer, goiter, hypertension, etc. [65].Spatial distribution of HI for adults (Figure 7a,b) show that most areas have little to no health risks except for a small part of the northwest in 2016.HI for children was generally higher than for adults (Table 7), whether in 2006 (Figure 7c) or 2016 (Figure 7d).In 2006, 91.8% of the study area had HI between 1.0 and 2.0, while this figure decreased to 34.4% in 2016.Additionally, in 2016, 58.4% had HI < 1.0, and 7.3% located at the northwestern corner had HI > 2.0.

Conclusions
In the present study, temporal and spatial variation of groundwater geochemical characteristics, vulnerability and health risk assessment were investigated in the Limin District.Statistical analysis associated with Piper diagrams was used to display the changes in groundwater geochemical characteristics from 2006 to 2016.A GIS-based modified DRASICL model was used to assess groundwater vulnerability, and AHP was applied to fix the weights of seven parameters.The health risk due to oral ingestion of contaminated groundwater was also conducted.
Groundwater geochemical characteristics from 2006 to 2016 changed slightly, and the major facies were HCO 3 -Ca•Na-type water in 2006 and HCO 3 -Ca-type water in 2016.Shallow groundwater was contaminated by NH 4 -N, Fe, and Mn.The average concentrations of Fe and Mn showed a significant difference between 2006 and 2016, and both exhibited a decreasing trend.Vulnerability maps of groundwater in 2006 and 2016 had similar spatial distribution with both showing high and very high-vulnerability areas.The percentage of high-and very high-vulnerability areas was greater in 2016 than in 2006.The high-and very high-vulnerability areas corresponded to agricultural land.The total risk from multiple contaminants ranged from 0.42 in 2006 to 0.31 in 2016 for adults, and ranged from 1.45 in 2006 to 1.06 in 2016 for children.Health risk of children was higher than adults.Although the concentration of NO 3 -N was less than drinking water standards, the average value of HQ NO3-N increased for both adults and children between 2006 and 2016.
Many factors affect the groundwater environment, including natural factors such as high background concentrations of Fe and Mn, and human activities such as groundwater exploitation.Meanwhile, the influence of natural factors on the groundwater environment is stable and can be very small, but the impact of human interference may be irreversible.Human activities change not only water quantity but also water quality.The assessment results of vulnerability and human health risk suggest that water resources management should pay attention to groundwater exploitation and agricultural non-point source pollution.Furthermore, introducing new irrigation techniques and enhancing water-saving measures is vital for sustainable development in agriculture-based districts.
Unfortunately, there are some limitations that need to be discussed in present study.For example, large amounts and accuracy of required data are one of the disadvantages of the assessment method applied.Furthermore, there is a need for strong evidence to explain the relationship between groundwater pollution and surface water pollution.Therefore, additional future research, like combining other analysis methods such as the isotope method, will be needed.Nevertheless, our work is an attempt to explore some useful information for local water resources management.

Figure 1 .
Figure 1.Study area location and distribution of sampling points (a) and geological section map (b).

Figure 1 .
Figure 1.Study area location and distribution of sampling points (a) and geological section map (b).

For 3 −
cations, groundwater had K + concentrations of 1.03~2.57mg/L in 2006 and 1.53~3.74mg/L in 2016.There was the same trend to the concentration of Na + , with an average of 46.5 mg/L in 2006 and an average of 47.3 mg/L in 2016.The average concentration of Ca 2+ decreased from 50.8 mg/L in 2006 to 37.5 mg/L in 2016.The concentration of Mg 2+ in 2006 was in the range of 7.12~22.82mg/L, whereas it was lower, ranging from 6.4 to 14.8 mg/L, in 2016.The concentrations of anions were lower than the maximum allowable values of drinking water quality set by the PRC and WHO.Cl − had a higher concentration in 2016, with an average of 19.4 mg/L.The concentration of SO 4 2− had significantly decreased from 2006 to 2016 with the average of 24.3 mg/L and 15.5 mg/L, respectively.HCO had the concentration in the range of 173~583 mg/L in 2006 and 98~423 mg/L in 2016.In general, the analyses above implied that changes in groundwater geochemical characteristics from 2006 to 2016 were not significant.Na + , Ca 2+ and HCO 3 − were major ions in the groundwater.Similarly, it can be found from Figure 2 that the aquifer was mainly dominated by HCO 3 -Ca•Na-type water and HCO 3 − Ca-type water in 2006 and 2016, respectively.

Water 2018 ,
10, x FOR PEER REVIEW 8 of 18The concentrations of NH4-N in large areas were 0.5~1.0mg/L in 2006 and it increased to 1.0~1.5 mg/L in 2016.Release of wastewater from factories, agricultural runoff, and domestic sewage constituted both point and non-point sources of NH4-N[49].The rapid development of the study area contributed to the increasing of NH4-N pollution from 2006 to 2016.There are many sources of pollution distributed along the Songhua River and Hulan River, including both point and non-point sources, resulting in large emissions of NH4-N.These sources have lowered the water quality to a point where it has reached the surface water quality standards IV grade (>1 mg/L) of the PRC[50][51][52].This may increase the risk of groundwater pollution with the decreasing river water quality[53,54].

Figure 2 .
Figure 2. Piper diagrams for the groundwater major ions in 2006 (a) and 2016 (b).

Figure 2 .
Figure 2. Piper diagrams for the groundwater major ions in 2006 (a) and 2016 (b).

Figure 2 .
Figure 2. Piper diagrams for the groundwater major ions in 2006 (a) and 2016 (b).

Figure 4 .
Figure 4. Pie charts of areas of land use types in 2005 (a) and 2015 (b).

Figure 4 .
Figure 4. Pie charts of areas of land use types in 2005 (a) and 2015 (b).

Figure 4 .
Figure 4. Pie charts of areas of land use types in 2005 (a) and 2015 (b).

Figure 5 .
Figure 5. Spatial maps of land use types (a,b); depth to groundwater (c,d); river distance (e); aquifer thickness (f); soil media (g); impact of the vadose zone (h); and hydraulic conductivity (i).

Figure 5 .
Figure 5. Spatial maps of land use types (a,b); depth to groundwater (c,d); river distance (e); aquifer thickness (f); soil media (g); impact of the vadose zone (h); and hydraulic conductivity (i).

Figure 7 .
Figure 7. Spatial distribution of HI for adults (a,b) and children (c,d) in 2006 and 2016.

Figure 7 .
Figure 7. Spatial distribution of HI for adults (a,b) and children (c,d) in 2006 and 2016.

Table 1 .
Groundwater vulnerability assessments conducted by modified DRASTIC model.

Table 2 .
Rating and weight scheme of DRASICL parameter maps.

Table 3 .
The intake values of exposure parameters and RfD for the risk assessment due to exposure.

Table 5 .
Pairwise comparison matrix for calculating weights.

Table 6 .
Proportion of groundwater vulnerability areas in 2006 and 2016.

Table 6 .
Proportion of groundwater vulnerability areas in 2006 and 2016.

Table 6 .
Proportion of groundwater vulnerability areas in 2006 and 2016.

Table 7 .
Assessment results of health risks for adults and children through drinking water intake.

Table 7 .
Assessment results of health risks for adults and children through drinking water intake.