Elsevier

Science of The Total Environment

Volume 408, Issue 24, 15 November 2010, Pages 6108-6116
Science of The Total Environment

Quantitative assessment of groundwater vulnerability using index system and transport simulation, Huangshuihe catchment, China

https://doi.org/10.1016/j.scitotenv.2010.09.002Get rights and content

Abstract

Groundwater vulnerability assessment has been an increasingly important environment management tool. The existing vulnerability assessment approaches are mostly index systems which have significant disadvantages. There need to be some quantitative studies on vulnerability indicators based on objective physical process study. In this study, we tried to do vulnerability assessment in Huangshuihe catchment in Shandong province of China using both contaminant transport simulations and index system approach. Transit time of 75% of hypothetical injected contaminant concentration was considered as the vulnerability indicator. First, we collected the field data of the Huangshuihe catchment and the catchment was divided into 34 sub areas that can each be treated as a transport sub model. Next, we constructed a Hydrus1D transport model of Huangshuihe catchment. Different sub areas had different input values. Thirdly, we used Monte-Carlo simulation to improve the collected data and did vulnerability assessment using the statistics of the contaminant transit time as a vulnerability indicator. Finally, to compare with the assessment result by transport simulation, we applied two index systems to Huangshuihe catchment. The first was DRASTIC system, and the other was a system we tentatively constructed examining the relationships between the transit time and the input parameters by simply changing the input values. The result of comparisons between the two index systems and transport simulation approach suggested partial validation to DRASTIC, and the construction of the new tentative index system was an attempt of building up index approaches based on physical process simulation.

Introduction

Environment protection is an important and difficult part of sustainable development, and groundwater remediation is an especially difficult aspect because of the prohibitive costs and time requirements. Groundwater vulnerability assessment has recently become an increasingly important environment management tool for local governments. It allows for better understanding of the vulnerabilities associated with the pollution of local groundwater sub areas, according to local hydrological, geological or meteorological conditions. There have been many index systems for groundwater vulnerability mapping, including DRASTIC (Aller et al., 1987), SINTACS (Civita and De Maio, 1997), GOD (Foster, 1987), AVI (Stempvoort et al., 1993), PI (Goldscheider et al., 2000), and GLA (Hölting et al., 1995). Of these, DRASTIC is the most popular. It has seven parameters including depth to water table (D), recharge (R), aquifer type (A), soil type (S), topography (T), impact of vadose zone (I) and conductivity (C). Vulnerability index is defined as a weighted sum of ratings of these parameters. It has been applied in many cases in central Japan (Babiker et al., 2005, Pathak et al., 2009), Europe (Martino et al., 2004, Crema et al., 1998, Panagopoulos et al., 2006), United States (Plymale and Angle, 2002, Fritch et al., 2000), South Korea (Lee and Choi, 1997, Kim and Hamm, 1999), South Africa (Lynch et al., 1997), Iran (Mohammadi et al., 2009) and China (Wang et al., 2006, Bazimenyera and Tang, 2008, Chen and Fu, 2003, Zhang et al., 2009). Valuable experiences have been gained but experiences showed DRASTIC ineffective and hard to validate since it was built on the experiences of hydrogeologists rather than quantitative physical process study of groundwater system.

One of the main reasons that index assessment systems have become popular is due to their little requirement on field data inputs. They only require ratings of each parameter range for each subarea. Parameters' values are separated into several ranges or types, each range or type is assigned a rating. Then the study area is divided into sub areas accordingly so that parameter values in each subarea are in the same range. All of the above steps could be carried out artificially according to hydrogeologists' experiences. Additionally, the calculation is quite simple; no physical details of the groundwater are required. Such features of index systems have made them a convenient approach dating back to the time when use of GIS and computers was not widespread. However, there are several disadvantages of index system that often questioned by researchers. First, index system is based on an assumption that the vulnerability is a linear superposition of the artificially specified parameters. But some studies have shown a nonlinear relationship between vulnerability and specific parameters (Neukum and Azzam, 2009). Second, the ratings and weights are all specified artificially by experts according to subjective judgments. This could introduce a huge subjective effect into the results (Frind et al., 2006, Gogu et al., 2003, Merchant, 1994). Third, values of the parameters are discretized into ratings, which may introduce additional errors. Therefore, generally speaking, the index system is quite rough. There need to be some quantitative studies on vulnerability indicators based on objective physical process study. However, the concept of vulnerability still needs to be more specific. The definition of vulnerability is not quite clear now even though it has been discussed for many years. There has not been any general accepted quantitative technical indicator to measure groundwater system vulnerability (Frind et al., 2006, Gogu et al., 2003, Merchant, 1994).

Examining the actual physical or chemical processes when the pollution occurred might be a better way of assessing the vulnerability of contaminants to the groundwater environment. Numerical simulations such as MT3DMS (Zheng and Wang, 1998) or Hydrus3D (Simunek et al., 2008) models might be better choice than index systems. Analytical solution or numerical simulation of transport model has been applied in vulnerability assessment in some studies (Jeannin et al., 2001, Stewart and Loague, 2004, Witkowski and Kowalczk, 2004, Voigt et al., 2004, Connell and Daele, 2003). Mostly one dimensional transport models were constructed containing different parameters, such as saturated hydraulic conductivity, thickness of the relevant geological layers and the dilution ratio. Neukum and Azzam (2009) constructed a Hydrus1D transport model considering soil type, thickness of soil layer and time varying recharge conditions. Compared with index systems, such simulation based vulnerability assessment approaches are still in their infancy. Mean transit time of the contaminant from the source to the target was used as an indicator (Witkowski and Kowalczk, 2004, Voigt et al., 2004, Connell and Daele, 2003), as well as variance of the contaminant's break through curve at the target and the ratio between the maximum concentrations (Brouyere et al., 2001). Neukum and Azzam (2009) defined four indicators to estimate vulnerability: transit time of 50% solute mass breakthrough at groundwater table; maximum solute relative concentration Cmax/C0 in the percolation water when entering the groundwater table to the injected mass or solute concentration; duration of breakthrough defined by the difference in 75% and 25% solute mass breakthrough at groundwater table; temporal shape of the breakthrough curve expressed with the quotient (t25%  t50%)/(t25%  t75%). During the simulation, these four indicators were calculated separately and each of them was provided for vulnerability assessment, but there is not an integrated index of these indicators yet.

Therefore, in this study, we tried to do vulnerability assessment in Huangshuihe catchment in Shandong province of China using both contaminant transport simulations and index system approach. The basic idea and main steps of our approach is following Neukum and Azzam's (2009) way, but we made several modifications according to the specific hydrogeological conditions of our study area. Transit time of 75% of hypothetical injected contaminant concentration was considered as the vulnerability indicator. First, 1), we collected the field data of Huangshuihe catchment. For lack of detailed field survey, the collected parameter distribution was regionalized into sub areas that each had uniform parameter value. Thus the Huangshuihe catchment was divided into 34 sub areas that each can be taken as a transport sub model. Second, 2), we constructed a Hydrus1D transport model according to the hydrogeological conditions of the Huangshuihe catchment. Different sub areas had different input values. Third, 3), we used Monte-Carlo simulation to improve the collected data and did vulnerability assessment using the statistics of the contaminant transit time as a vulnerability indicator. Finally, 4), to compare with the assessment result by transport simulation, we applied two index systems to Huangshuihe catchment. One was DRASTIC system. The other was tentatively built up by us by examining the relationships between the transit time and the input parameters by simply changing the input values. The result comparisons of the two index systems and transport simulation approach suggested partial validation of DRASTIC. The construction of the new tentative index system was an attempt of building up index approaches based on physical process simulation (Fig. 1).

Groundwater vulnerability consists of two parts, intrinsic vulnerability and extrinsic vulnerability (Villa and McLeod, 2002). Intrinsic vulnerability is independent of the nature of contaminants and the contamination scenario. It takes into account the geological, hydrological and hydrogeological characteristics (Zwahlen, 2004). Extrinsic vulnerability contains factors external to the system, incorporates the two components of risk: exposure and hazard acting on it (Villa and McLeod, 2002). According to different human activities, we think there are two types of vulnerabilities, from contamination and from being exhausted by over exploitation. In practice, including in this study, what concerns us most is the vulnerability from contaminants. Our focus was on the solute transport process from the surface through the vadose zone during rains. Once the contaminant reached the water table, it indicated groundwater had been contaminated. This is a reasonable assumption, especially for cases on plains or catchments in which the groundwater table is shallow and there exists a large population density. Usually, such areas exist in areas with highly developed economies, such as Yangtze River Delta, Zhujiang River Delta, and the Bohai Sea Ring Area of China.

We took root water uptake, soil thickness, soil media property, and precipitation flux into consideration in our case study. Root water uptake was affected by vegetation, closely related to land use. It was considered because we wanted to take the human's agricultural activities into consideration to represent the impact of different land use. Soil thickness, soil property and infiltration flux are parts of intrinsic vulnerability factors. Precipitation flux was not considered as time varying, but a hypothetical constant average precipitation rate during rains. We set up a transport conceptual model under a hypothetical steady state to examine the intrinsic potential vulnerability of hydrogeological conditions. The transit time of 75% of the hypothetical contaminant concentration from the surface to the water table was taken as the vulnerability indicator to measure vulnerabilities. It's more vulnerable when transit time is smaller.

The governing equation of Hydrus1D describing the water flow movement in the unsaturated zone is a modified form of the Richard Eq. (1):θt=xKhx+1Swhere h is the water pressure head, θ is the volumetric water content, t is time, x is the spatial coordinate (positive upward), S is the sink term, and K is the unsaturated hydraulic conductivity function. We took each sub area as a 1D numerical model. K, h, and θ has relationship between each other according to van Genuchten Eqs. (2) and (3) (van Genuchten, 1980),θ(h)={θr+θsθr(1+αhn)mh<0θsh0K(h)=KsSel[1(1Sel/m)m]2where m = 1  1/n; (n > 1), θr and θs are the residual and saturated water content respectively. α, n, and m are parameters for soil properties, Se is the effective Saturation and l is the pore-connectivity parameter. Hysteresis in the retention curve and in the hydraulic conductivity is not considered in this study. The sink term S is defined as the volume of water removed from a unit volume of soil per unit time due to plant water uptake(Eq. (4)).S(h)=α(h)b(x)Tpwhere α(h) is a prescribed dimensionless function of the soil water pressure head (0  α  1), b(x) is the root distribution function, and Tp is the potential transpiration rate (further referred as root water uptake rate in this paper) (Šimůnek et al., 2009). The root distribution function b(x) followed a modified right trapezoidal distribution (Eq. (5)).b(x)={43LRx<0.5LR83LR8x3LR2x(0.5LR,LR)0x>LR

The lower boundary condition was a constant water content boundary because it was saturated at the water table. And the upper boundary was specified as constant flux boundary with the value of average precipitation rate during rains. We considered the flow model as a steady state model during the rain. Average flux during the rain was more appropriate than average annual precipitation. The differential equation for solute transport is a general advective–dispersive Eq. (6), ignoring chemical reactions or source-sink item during the transport.t(θc)=xθDcxx(qc)where c is the solute concentration, D is the dispersion coefficient and q is the volumetric fluid flux density calculated from the flow model described by Eq. (1). Absorption and hysteresis were not considered in this hypothetical case in our study. Crank–Nicholson scheme was used for time weighting and Galerkin finite elements for space weighting. The upper boundary condition for contaminant transport model was specified as a hypothetical constant concentration boundary with the value of 1 g/cm³, but the results were analyzed in terms of relative concentration instead of absolute concentration.

The study area is the Huangshuihe catchment in Shandong province of China. This catchment is located in the north of Shandong Peninsula, with an area of 1034.47 km2, 20% of which is plain. Huangshuihe catchment covers an area of 51.8 hectare, total capacity of 5359 × 104 m3 and maximum usable storage of 3929 × 104 m3. The catchment lies to the south of Bohai Sea and is surrounded by uplands with elevations between 1.7 m and 43 m in the east, south and west. The Huangshuihe River lies in the middle of the catchment and flows north straight to Bohai sea (Fig. 2). Rainfall infiltration is the main source of groundwater recharge. Most of the groundwater in the Huangshuihe catchment is Quaternary groundwater. The downstream area of the Huangshuihe catchment was taken as the study area, covering 30.9% of the whole catchment. Shandong, known for its agricultural economy, is now facing a serious shortage of water resources as well as non-point source pollution by agricultural activities. The Huangshuihe catchment, one part of Yantai, which is famous for the “Yantai Apple,” supplies a population of 0.63 million and GDP of 8 billion dollars. Locating on the edge of the North China Plain, Shandong is in an arid and semi-arid region. Groundwater is vitally important to the economic survival of this region. Groundwater vulnerability assessment is now an important part of local groundwater management for local government.

The soil thickness, also the depth from the surface to the groundwater table, varies mainly from less than 2 m to over 10 m in this area. Since the depth of topsoil is about 20–30 cm and subsoil is almost the same vertically, the soil profile is assumed homogeneous according to the subsoil. For the lack of monitoring wells, depths of water table were just recorded at several points, which created the difficulty of producing a depth map by interpolation. Therefore, the values of depth were divided into three ranges (1.5–4.5 m; 4.5–9 m; 9–15 m). Accordingly, the study area was divided into four sub areas for depth (Fig. 3).

Similarly, the soil type of this area is mainly sandy loam, with a sand section in the north adjacent to the Bohai Sea, and another section of loamy sand along the Huangshiuhe in the center of the area. There are five typical types of land use in this area, including orchard, vegetable farm, wheat country, grass land and residential area (Fig. 3). The grass land is distributed mainly along the coast, and the residential area is in the southwest. This area contains plenty of agricultural activities, with the apple as the most well-known product. Feddes parameters for oranges were applied for our numerical model instead of those of apples because of the lack of a previous study on apple trees. Many kinds of vegetables were planted in the south section of the area near the residential section, but we chose celery to typify all the vegetables. Different land use could affect root water uptake rate and root depth affect the depth of water uptake, and root water uptake would affect a lot on the available percolation water. Since solute transit time is heavily depending on the percolation rate, land use may have a significant impact on transit time.

Therefore, four typical plants were chosen to be simulated and represent the effects of different land uses on transport in the vadose zone. The root water uptake rates for these plants were specified according to previous study of Zhang et al. (1997) for vegetables, Gong and Kang (2004) for apple trees, Zhang et al. (2008) for wheat and Wang et al. (2003) for grass. Because these four plant types were chosen to typify the four land uses, the corresponding root water uptake rates were used for the four land uses.

According to the local average precipitation rate during rains, the specified upper boundary condition for each subarea varied from 1 cm/day to 3.3 cm/day, and a hypothetical injected contaminant concentration of 1 g/cm³ was specified as the upper boundary concentration condition. Elements of our model, soil types and Feddes parameters for the plants are listed in Table 1. The study area was divided into 34 subareas. Every sub area was taken as a 1D model with the same parameters but different input values.

For lack of sufficient accurate field data, parameter values of depth, infiltration rate, root water uptake rate, and soil media properties were all of significant error. Values were divided into ranges but not continuously distributed. Accordingly, the study area was regionalized into subareas. Therefore to improve the field data and reduce the error caused by such discretization, we used Monte-Carlo simulations to carry out the transport modeling for each sub area. Depth and infiltration rate values were assumed to follow a uniform distribution and Gaussian distribution was assumed for other parameter values (Table 1). Mean, range and standard deviation were estimated according to the collected field data and experiences. Vulnerabilities for sub areas were assessed from the mean transit time of Monte-Carlo simulations.

During the Monte-Carlo simulation for each sub area, there were one hundred samplings by acceptance–rejection sampling method for Gaussian distributed parameters and fifty samplings for uniform distributed parameter. Mean of the transit time for 75% of the injected concentration were calculated as the vulnerability indicator. Relative vulnerability index for sub area i was calculated by normalization (Eq. (7), Table 2) for comparison with the results of other approaches carried out later in our study.VulRi=TmaxTiTmaxTmin×100

Besides the Monte-Carlo simulation carried out to assess the vulnerabilities for sub areas of Huangshuihe catchment, we also applied two index systems. One was DRASTIC system. The other was tentatively built up by us by examining the relationships between the transit time and the input parameters by simply changing the input values. There was an assumption made that the effective factors to vulnerability were independent of each other, which was certainly too strong for representing the true physical process, however, it might be applicable just for index improvement which was only concerned with the low precision proportional relationship between them. We attempted to take as many other factors into consideration when calculating index for a certain parameter.

Since the manner in which depth affects the transit time is already well known (Neukum and Azzam, 2009), the depth was set to 2 m in this step. The ranges and steps of examined parameters shown in Table 3 are for relationship analysis and index system improvement, not for vulnerability assessment. Further vulnerability assessment was carried out by the tentatively improved index system. Soil properties for soil types that did not exist in Huangshuihe catchment were specified as Hydrus1D recommended. Other parameter values were specified as the same in Monte-Carlo simulation. Note that parameter values in this step were not according to actual conditions of Huangshuihe catchment. Transit time of 75% injected concentration was saved as T(S, L, I, Rroot), where S was soil type, L was land use, I was infiltration rate and Rroot was root water uptake rate.

Transit times of 75% injected concentration in each soil media type and land use were calculated. The results in orchard were plotted as an example in Fig. 4. An inverse relationship with transit time and infiltration rate was appropriate in all the soil types. Relative vulnerabilities of different soil media could be estimated from their transit times. Assuming that the average transit time in soil media S1 and S2 were T1 and T2, the relative indexes of S2 to S1 could be T1 to T2. Transit time in silt loam under the minimum infiltration rate in wheat country under the max root water uptake rate T(silt loam,wheet,1.1,Rrmax) was the largest of all the transit times, therefore it's considered as the least vulnerable condition. Taking T(silt loam,wheet,1.1,Rrmax) as the reference, vulnerability index in each soil media under each infiltration rate was calculated as an average of other parameters (Eq. (8)).Indicator(S,I)=RrootLT(siltLoam,wheet,1.1,Rrmax)T(S,L,I,Rr)NLNRrootwhere NL, NRroot were the number of land uses and root water uptake rate sampling points. To determine the vulnerability index for single parameter of soil type S, the evaluated index was calculated as about three times of the average of all indicators in soil S so as to make the evaluated index integers (Table 4).IndexS=Round3×IIndicator(S,I)NSwhere NS was the number of soil media types. The Round() function returned the nearest integer.

Similarly, in order to check the effects of different land uses, transit times in different soil types and all land uses were calculated. Sand, sandy loam and silt loam were plotted to typify soil media types that they could respectively represent the good, medium and the bad measurements of permeability (Fig. 5). The maximum transit time was taken to be a reference again as the least vulnerable condition. Vulnerability index of I was taken as its value (Eq. (10)) and indexes for land use were calculated using Eq. (11).IndexI=IIndexL=Round3×S,I,RrootT(siltLoam,vegetable,Imin,Rrmax)T(S,L,I,Rr)NSNINRr

To check the sensitivity with transit time and root water uptake rate, results in loamy sand was plotted as an example (Fig. 6). Transit time almost increased along with Rroot2. Relative transit times of different land uses almost still kept their proportions. In this study, we simply set the index of root water uptake rate to be Rroot2 (Eq. (12)).IndexRroot=Rroot2

Depth was a direct ratio of transit time. Therefore, it was an inverse ratio of vulnerability. If the soil type's and land use's indexes of one subarea i was IndexiS and IndexiL, the vulnerability and the relative vulnerability index of subarea i was estimated by Eq. (13) and Eq. (14), where Vulmax was the maximum of all Vuli, Di, Riroot and Ii were the value of infiltration rate, root water uptake rate and depth of sub area i. Parameter indexes and vulnerability index calculation of our tentative model were listed (Table 5), Local data of Huangshuihe catchment was applied to our index system (Table 6).Vuli=Ii(Rrooti)2Di×IndexSi×IndexLiVulRi=VuliVulminVulmaxVulmin×100

Finally DRASTIC index system was also applied. Different from average flux during the rain used in transport simulation of this study, infiltration rate was specified as average annual precipitation as recommended in DRASTIC. Weights of DRASTIC were specified as scheme for agricultural areas recommended by Aller et al. (1987). For comparison with other two sets of relative vulnerabilities, relative vulnerabilities of DRASTIC were also recalculated using Eq. (13). Thus, local data of Huangshuihe catchment was applied to Hydrus1D simulation approach (Table 2), our tentative index system and DRASTIC (Table 6).

The results of DRASTIC turned out to be right in the middle of the two set result of other approaches. Correlation plot (Fig. 7) showed a rough accordance between results of DRASTIC and the other approaches. The result of transport simulation approach could provide partial validation to DRASTIC. Though they were much different in values, but quite similar in tendency. From the vulnerability maps by transport simulation and DRASTIC assessment results, sub areas looked similar, while they looked less vulnerable by the new index system (Fig. 8). Although we failed to find a way to help in validating which result was better in this study, authors still thought this was a good attempt to find other approaches besides index systems using quantitative process study.

Section snippets

Conclusion and discussion

In this research we applied the Neukum and Azzam's (2009) vulnerability assessment approach by transport process simulations in Huangshuihe catchment in Shandong province of China. In this study we were concerned with such plains or catchments exhibiting high levels of agriculture and industry, where there is more concern for groundwater vulnerability than in other areas. Our focus was on the solute transport process during the rain. Average precipitation rate during the rain, land use types or

Acknowledgements

This research was supported by international scientific and technical cooperation and exchange program of the Department of Science and Technology of China (2007DFB70200) and the Major Project of Chinese National Programs for Fundamental Research and Development (2006CB403404). We are grateful to Water Affairs Bureau of Longkou in Shandong province of China for their help on field data support.

References (41)

  • E.O. Frind et al.

    Well vulnerability: a quantitative approach for source water protection

    Ground Water

    (2006)
  • T.G. Fritch et al.

    An aquifer vulnerability assessment of the Paluxy Aquifer, Central Texas, USA, using GIS and a modified DRASTIC approach

    Environ Manage

    (2000)
  • R.C. Gogu et al.

    Comparison of aquifer vulnerability assessment techniques. Application to the Ne'blon river basin (Belgium)

    Environ Geol

    (2003)
  • N. Goldscheider et al.

    The PI method: a GIS based approach to mapping groundwater vulnerability with special consideration of karst aquifers

    Z Angew Geol

    (2000)
  • Gong D.Z., Kang S.Z., 2004. Measuring and Estimating Evaportranspiration of an Apple Orchard (in Chinese), Journal of...
  • B. Hölting et al.

    Conception for the evaluation of the protective function of the unsaturated stratum above the groundwater table

    Geol Jahrb

    (1995)
  • P.Y. Jeannin et al.

    VULK: a tool for intrinsic vulnerability assessment and validation

  • Y.J. Kim et al.

    Assessment of the potential for groundwater contamination using DRASTIC/EGIS technique, Cheongju area, South Korea

    Hydrogeol J

    (1999)
  • S. Lee et al.

    Groundwater pollution susceptibility assessment of Younggwang area using GIS technique (in Korean)

    J Korean Soc Groundw Environ

    (1997)
  • S.D. Lynch et al.

    A DRASTIC approach to groundwater vulnerability in South Africa

    South African J Sci

    (1997)
  • Cited by (0)

    View full text