Application of grey correlation method to evaluate potential groundwater recharge sites

(1)

O R I G I N A L P A P E R

H. S. Gau Æ C. Y. Hsieh Æ C. W. Liu

Application of grey correlation method to evaluate potential groundwater

recharge sites

Published online: 9 February 2006 Springer-Verlag 2006

Abstract Artificial recharge is a practical tool available for increasing the groundwater storage capacity. The efficiency of artificial recharge is related to various hydrogeological factors of the target area. In this study, a variable saturated groundwater flow model, FEM-WATER, was used to evaluate the arrival times of recharged water that infiltrates from an artificial recharge pond to the groundwater table under various hydrogeological conditions. Forty-five arrival times were generated by FEMWATER. The relationships between the arrival times and hydrogeological factors used in the simulation of FEMWATER were analyzed by the grey correlation method. The results show the order of importance of the factors as they influence the arrival time. In order from high to low importance, they are a, Dg, he, Dp, KSand b. Dgand Dpare interpreted as the potential for movement of the recharge water; heis the water storage capacity of soil, and KSrepresents the ability of soil to transport water. a and b describe the characteristic curve of the unsaturated soil. The method was applied to evaluate a suitable site for artificial recharge in the Yun-Lin area. Grey correlation analysis was performed to obtain the grey correlation grade using the minimum arrival time as a reference sequence. An index is proposed herein to determine the recharge efficiency of 20 sampling sites. A contour mapping of index values at the 20 sampling sites identified three areas for artificial aquifer recharge in Yun-Lin. Area A in the upper plain is considered more appropriate for groundwater recharge than areas B and C in the coast.

Keywords Artiﬁcial recharge Æ Grey correlation analysis Æ Site selection Æ Numerical simulation

1 Introduction

The use of surface water during the wet season and groundwater during the dry season is one means of efficiently utilizing limited water resources. Artificial recharge methods can be applied to increase the groundwater storage capacity in the wet season. Todd (1980) and Bouwer (1996, 2002) introduced several artificial recharge techniques, such as well recharge, and the use of infiltration basin and seepage trenches. An effective artificial recharge method must suit site-specific conditions, including soil permeability, the properties of the unsaturated soil, geological structure, hydrological conditions, topography and land use. Recently, the photographic mapping has been used to determine the ambit of land use, the spread of soil types and the geo-logical structure, based on satellite images, which have thus been used to identify potential recharge areas (Salama et al. 1994; Saraf and Choudhury 1998). Rastogi and Pandey (1998) compared the efficacy of artificial recharge basins with different shapes by finite element modeling. They pointed out that a rectangular artificial recharge basin can store more water in an aquifer than basins of other shapes. Moreover, Gau et al. (2001) investigated the effects of relevant hydro-geological factors, including the depth of the artificial recharge pond water, the hydraulic conductivity and the depth of the groundwater table, on the arrival time of the recharge water to the groundwater table. They concluded that an appropriate site for artificial recharge is not determined by only the hydraulic conductivity.

A suitable artiﬁcial recharge site must have highly permeable soil, with no impermeable layer, an uncon-ﬁned aquifer and an absence of any contamination source (Bouwer1996). Methods for estimating the rate or amount of groundwater recharge have been devel-oped. The simple water balance model has been widely

H. S. Gau Æ C. Y. Hsieh

Department of Environmental Engineering and Science, Tajen Institute of Technology, Ping-Tung, 907 Taiwan, R.O.C C. W. Liu (&)

Department of Bioenvironmental Systems Engineering, National Taiwan University, Taipei, 106 Taiwan, R.O.C

E-mail: lcw@gwater.agec.ntu.edu.tw DOI 10.1007/s00477-006-0034-9

(2)

applied to estimating the amount of groundwater recharge, but any estimate made using this model is sensitive to the forest canopy parameters, and parame-ters in the soil moisture model, particularly the rooting depth, the fractional available water content, and the fractional field of drainable water (Finch 1998). Other simple soil models can also be used to estimate water recharge, by calculating how much water is potentially available for groundwater recharge from the differences among the measured rainfall, the estimated actual evapotranspiration and the soil moisture (Ragab et al. 1997; Chapman and Malone 2002). As well as the aforementioned physical models, the measured distri-bution of environmental tracers, 3H,36Cl and chloride, in a soil water profile, have also been used to estimate the rate of movement of water in soil, to locate a groundwater aquifer (Eriksson and Khunakasem 1969; Sukhija et al. 1997; Cook and Robinson2002).

All of the above methods are appropriate for esti-mating the recharge rate at a site where sufficient in situ data have been measured in the field. However, if no field data are available to estimate the recharge rate, then the use of results simulated using a groundwater model is an effective means of evaluating the feasibility of a site for artificial recharge. Many groundwater sim-ulation models have been developed over the last 20 years for simulating the response of a groundwater head to surface infiltration (van Dam and Feddes2000; Warrick1991; Giao et al.1999; Parlange et al.1999). A variable saturation model is required to simulate the movement of the wetting front of recharge water from unsaturated soil to the saturated ground water aquifer. The recharge water contributes to the groundwater only if the wetting front reaches the groundwater table. Accordingly, the efficacy of an artificial recharge site can be measured by comparing the required arrival time of the wetting front at groundwater table.

Hydrogeological factors, including hydraulic con-ductivity, the depth of water in an artificial recharge pond, the depth of the groundwater table and the unsaturated soil characteristic curve, are required to simulate variable saturated groundwater flow. Some factors may be easily measured in the field, such as hydraulic conductivity, depth of water in the artificial recharge pond and depth of the groundwater table. Other factors, such as parameters that govern the characteristic curve of unsaturated soil, can only be determined with extra time and labor, and expensive instruments. Therefore, the hydrogeological factors that can be measured in the field depend on the budget of the project and the time available. The order of importance of these factors to the movement of recharge water must be established to determine the priority of these factors to be measured in the field.

In this work, the variable saturated groundwater model, FEMWATER, was used to simulate the arrival time of recharge water at the groundwater table and to determine the relationship between the arrival time and each hydrogeological factor. The order of importance of

the hydrogeological factors to the recharge of ground water was evaluated by the grey correlation method (Deng1989). The determined order of importance of the hydrogeological factors coupled with the grey correla-tion analysis was then applied to locate the optimal recharge site in Yun-Lin plain in Taiwan. The results help engineers not only to identify the important factors that need to be measured in the field, but also offer an effective means of evaluation in a preliminary screening of a groundwater recharge site.

2 Materials and methods

2.1 Numerical simulation model

The arrival time of the recharged water front at the groundwater table is an important factor that governs groundwater recharge. Various methods can be used to determine the arrival time, including laboratory, field and numerical experiments. Laboratory and field experiments are expensive and time-consuming, so numerical experiments are conducted to simulate the movement of a wetting front and estimate the arrival time of artificial recharge water at the groundwater table.

The Richards’ equation is commonly used to describe the movement of inﬁltration through an unsaturated zone of soils. It is given as follows (El-Hames and Richards1995)

CðhÞ@h

@t ¼ r KðhÞ rðh þ zÞ½ þ q; ð1Þ where C(h) represents the diﬀerential water capacity (¶h/¶h) (L1); h is the water content (L3L3); h is the water pressure head; t is the time (T); K is the unsatu-rated hydraulic conductivity; z is the potential head (L), and q is the source and/or sinks (L3L3 T1).

Many numerical groundwater models have been developed based on Richards’ equation to simulate the distribution of pressure head and water content, or the saturation of unsaturated soil. In this work, a three-dimensional ﬁnite element computer model, FEM-WATER, was used to simulate density-dependent ﬂow and transport in variably saturated media. (Lin et al. 1997) This model is incorporated into a visualized groundwater modeling system, GMS (EMRL2003), which is adopted herein to perform a numerical experiment.

A homogenous and isotropic soil column of size 14 m·14 m·30 m and a pond of size 12 m·12 m·3 m on top of it, are considered in the numerical experiment (Fig.1). The depth of water in the pond and the groundwater table is maintained at a constant level throughout the numerical simulation. Impermeable walls are set around the soil column and recharge water ﬂow inﬁltrates vertically. The soil column is discretized into 2,790 prism elements with 1,936 nodes. The size of each prism element is 2 m·2 m·1 m (Fig.1).

(3)

The characteristic curve of unsaturated soil is deter-mined by the Mualem-van Genuchten model equations (van Genuchten1980):

Se¼ 1 þ ðawÞh bi1=b1; ð2Þ

KðSeÞ ¼ KSpffiffiffiffiffiSeð1 ð1 S_eb=b1Þ11=bÞ2; ð3Þ where K(Se) (L T1) is the unsaturated hydraulic conductivity; KS (L T1) is the saturated hydraulic conductivity; w is the absolute value of the matrices potential (M T2L1); Se¼ hhr

hshr is the degree of

satu-ration (dimensionless); h is the water content(L3 L3); hS is the saturated water content (L3L3); hr is the residual water content (L3L3),and a (L1) and b (dimensionless) are parameters obtained by curve fitting (Luckner et al. 1989). hr is defined as the value of h at which fluid flow in response to hydraulic gradients ceased, whereas hsis measured experimentally.

The fitting parameters a and b are used to describe the shape of the water-release curve. a is proportional to the reciprocal of the value of w halfway between hSand hr. Small a values represent little change in water content as w becomes more negative, and are generally observed in fine-grained and unstructured soils. Large values of a indicate a rapid change in water content as some pores empty out under very small negative heads. They are typically found in sands or well-structured soils. The parameter b specifies the steepness of the water-release curve. If b is large, then the curve is steep, and the water content decreases quickly as w becomes more negative. If b is small, then the change in water content is slow (Hodnett and Tomasella2002). Five soils, silty clay loam (SCL), yolo light clay (YLC), sand, loamy sand and loam, are selected as the materials in the simulated soil column. Table 1presents the parameters of the five soil characteristics curves. Kao et al. (2001) experimentally measured the parameters of the first three soils. The

14m×14m×30m soil column 10m×10m×3m artificial recharge pond

1m

Fig. 1 Soil column for modeling recharge test by FEMWATER in the numerical experiment

(4)

other parameters were generated using the distribution functions suggested by Hubert et al. (1996).

Constant depths of pond water (Dp) of 1, 2 and 3 m were maintained throughout the simulations. The depths between the bottom of the pond and the groundwater table (Dg) were also ﬁxed at 10, 20 and 30 m. Other sources and sinks, including rainfall inﬁltration, the pumping well and evapotranspiration, were not consid-ered herein. The initial pressure heads were set to 0.5 m in all grids. Under the simulation conditions, a total of 45 (3·3·5) cases were considered.

The FEMWATER outputs the simulated total head and pressure heads every hour. The results of the pres-sure heads were plotted as contours to trace the position of the zero-pressure. The contour of the zero pressure head intersected the groundwater table, and the arrival time of the recharged water front was determined accordingly.

2.2 Grey correlation analysis

A sensitivity analysis is conducted to evaluate the importance of the hydrological factors for the movement of recharge water. The traditional sensitivity analysis is a ‘‘change-one-factor-at-a-time method’’: one factor is varied while all others are ﬁxed. However, the simplicity of the method may yield unreliable results and an inadequate conclusion because factors are uncertain (Deng 1996; Chang 2005). Deng (1989) presented a method of grey correlation analysis, which addresses some of the shortcomings of the traditional method. The basic concept of grey correlation analysis is to determine whether a relationship among a series of data is close, based on the degree of similarity among the geometric shapes of the data series’ curves. Closer curves indicate a greater the correlation among the relative data series. The grey correlation grades are determined from the grey correlation coeﬃcients, which measure the degree of similarity among sequences. The order of sensitivity of the hydrological factors to the movement of recharge water can then be determined based on the calculated grey correlation grades. The method has been exten-sively applied across many areas, such as auto control design, water resource development, hydroelectric gen-eration scheduling and environmental management.

(Fung 2003; Yeh and Chen 2004; Liang 1999; Chang and Lin1999).

2.2.1 Data preprocessing

Grey data processing must be performed before grey correlation coeﬃcients can be calculated. A series of various units must be transformed to be dimensionless. Usually, each series is normalized by dividing the data in the original series by their average.

Let the original reference sequence and sequence for comparison be represented as xo(k) and xi(k), i=1, 2, ..., m; k=1, 2, ..., n, respectively, where m is the total number of factors to be considered, and n is the total number of observation data. Data preprocessing con-verts the original sequence to a comparable sequence. Several methodologies of preprocessing data can be used in grey relation analysis, depending on the characteris-tics of the original sequence (Fung 2003). If the target value of the original sequence is ‘‘the-larger-the-better’’, then the original sequence is normalized as follows;

x_iðkÞ ¼ xiðkÞ min xiðkÞ

max xiðkÞ min xiðkÞ: ð4Þ

If it is ‘‘the-smaller-the-better’’, then the original sequence is normalized as follows;

x_iðkÞ ¼ max xiðkÞ xiðkÞ

max xiðkÞ min xiðkÞ: ð5Þ

However, if there is a speciﬁc target value, then the original sequence is normalized using,

x_iðkÞ ¼ 1 jxiðkÞ OBj

maxfmax xiðkÞOB; OB min xiðkÞg; ð6Þ where OB is the target value.

2.2.2 Grey correlation coeﬃcient and grey correlation grade

After the data have been preprocessed, a grey correla-tion coefficient is determined using the preprocessed sequence. The grey correlation coefficient is defined as (Deng1989; Fung2003; Yeh and Chen2004),

Table 1 Mualem-van Genuchten parameters of the soil materials

SCLSilty clay loam, YLC yolo light clay, SC sand clay

Soil material type Saturated hydraulic conductivity (KS) (m/days) Residual water content (hr) Saturated water content (hs) a (1/m) b SCL 0.0624 0.1000 0.4100 1.90 1.31 YLC 0.4320 0.2300 0.5500 3.60 1.90 SAND 2.4000 0.0500 0.5000 3.70 5.00 LOAM 0.2520 0.0776 0.4300 3.67 1.56 SC 0.0307 0.0993 0.3800 2.70 1.28

(5)

cðx_oðkÞ; x_iðkÞÞ ¼ Dminþ nDmax

DoiðkÞ þ nDmax; ð7Þ

where n is the distinguishing coeﬃcient, n 2(0, 1). Generally n is taken as 0.5 (Deng 1989). Doi(k) is the sequence of deviation of the reference sequence xo*(k) from the sequence xi*(k) for comparison; DoiðkÞ ¼ x oðkÞ xiðkÞ ; Dmax ¼ max 8j2i max8k x oðkÞ xjðkÞ

and Dmin ¼ min 8j2imin8k x

oðkÞ xjðkÞ

: The grey correla-tion grade is an average of the grey correlacorrela-tion coef-ﬁcients and is deﬁned as (Deng 1989),

Ciðxo; xiÞ ¼ 1 n Xn k¼1 cðx_oðkÞ; x_iðkÞÞ: ð8Þ

Herein, the grey correlation grade represents the degree of correlation between the reference sequence and the sequence for comparison. If a particular sequence for comparison is more important to the reference sequence than other comparable sequences, then the grey corre-lation grade between that comparable sequence and reference sequence exceed the other grey correlation grades. If the comparable sequence is consistent with the reference sequence, then the grey correlation grade is close to unity. Grey relational analysis measures the absolute diﬀerences between sequences, and so can be

used to measure the approximate correlation of sequences (Fung2003).

3 Results and discussion

3.1 Numerical experiment

Figure2 plots the simulated arrival times in the numerical experiment. The ordinate represents the arri-val time of recharge water at the groundwater table. The abscissa represents the depth of the pond water. The solid line with rectangular, triangular and circular symbols represents the groundwater tables 10, 20 and 30 m underground, respectively. Four vertical solid lines divide the figure into five zones based on the type of soil. The results indicate that the arrival time increased with the depth of the groundwater table and fell as the depth of the pond water increased, regardless of the soil. Comparing the arrival times in different soils shows that the arrival time is shortest in the sand material and longest in the SCL material. The difference between the shortest and longest arrival times increases with the depth of the groundwater table.

The arrival time is used to evaluate the efficiency of artificial recharge. A short arrival time represents highly efficient artificial recharge because the pond water easily infiltrates the soil and recharges to groundwater aquifer. Conversely, a long arrival time represents a poor

eﬃ-1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 Dp(m) -200 0 200 400 600 800 1000 1200 1400 1600 A V (hr)

SCL YLC SAND LAOM SC

Dg= 10m

D_g= 20m

Dg=30m

(6)

ciency of artificial recharge, because most of the infil-trated pond water is retained in the unsaturated zone and cannot quickly be replenished to the groundwater table. If the required arrival time exceeds the duration of the wet season, then no surface water is available for recharge in the dry season. Therefore, the arrival time is the key to locating an artificial recharge pond. Clearly, the results in Fig.2demonstrate that a suitable site for an artificial recharge pond has a shallow groundwater table, permeable soil like sand, and a deep pond water. Unfortunately, such conditions are not always available. For instance, a permeable soil may have a deep groundwater table. If a hydraulic engineer evaluates a potential site for an artificial recharge pond based solely on a single hydrogeological factor, such as the hydraulic conductivity of soil, an incorrect decision may be made. Hence, grey correlation analysis is useful for calculating the grey correlation grade of influential factors. With reference to the order of the grey correlation grade, the sensitivity of each factor to the movement of recharge

water is analyzed and an optimal site for constructing an artiﬁcial recharge pond can be determined without bias.

3.2 Most inﬂuential factor

Six influential factors are considered herein. Table2 presents the influential factors and factor levels at which the effective water content, he, is defined as he= hS hr. Figure3displays the arrival times of recharge water (AV) at the groundwater table for a total of 45 different factor levels. The arrival times of the 45 numerical runs are assigned to the original reference sequence xAVðkÞ; k=1, 2,

..., 45. The six inﬂuential factors in 45 runs are set as the comparable sequence xDgðkÞ; xDpðkÞ; xKSðkÞ; xheðkÞ;

xa(k), xb(k), k=1, 2, ..., 45.

Figure4is a scatter plot of the factors he, KS, a and b versus arrival time. The results indicate that the arrival time increases as the factors decrease. Figure2indicates that the arrival time increases with Dg and falls as Dp increases. A shorter arrival time generally corresponds to a higher artiﬁcial recharge. Thus, the arrival time is ‘‘the-smaller-the-better’’. The smaller-the-better method, Eq. 5, was applied to preprocess the data. The factor sequence he, KS, a, b and Dpare ‘‘the-larger-the-better’’. Hence, the ‘‘the-larger-the-better’’ method, Eq. 4, was employed to the factors he, KS, a, b and Dp for data preprocessing.

Figure5plots the results of data preprocessing. x_A_VðkÞ and xi*(k), i = Dg, Dp... a, b denote the reference sequence and comparable sequence, respectively. The deviation sequences DAVDgð1Þ; DAVDpð1Þ; DAVKSð1Þ; DAVheð1Þ; DAVa

ð1Þ; DAVbð1Þ are calculated as follows;

0 200 400 600 800 1000 1200 1400 1600 0 5 10 15 20 25 30 35 40 45 A V (hr)

Numerical experimental run number

Fig. 3 Simulated arrival times of artiﬁcial recharge water at ground water table Table 2 Experimental factor values corresponding to factor levels

Levels of factors Factors Dg(m) Dp(m) Soil type Property of soil KS (m/days) he a (1/m) b 1 10 1 SCL 0.0624 0.3100 1.90 1.31 2 20 2 YLC 0.4320 0.3200 3.60 1.90 3 30 3 SAND 2.400 0.4500 3.70 5.00 4 LOAM 0.2520 0.3524 3.67 1.56 5 SC 0.0307 0.2807 2.70 1.28

(7)

DAVDpð1Þ ¼ x AVð1Þ x Dpð1Þ ¼ 0:7241 DAVDgð1Þ ¼ x AVð1Þ x Dgð1Þ ¼ 0:2579 DAVheð1Þ ¼ x AVð1Þ x heð1Þ ¼ 0:5510 DAVKSð1Þ ¼ x AVð1Þ x KSð1Þ ¼ 0:7107 DAVað1Þ ¼ x AVð1Þ x að1Þ ¼ 0:7241 DAVbð1Þ ¼ xAVð1Þ x bð1Þ ¼ 0:7160: ð9Þ

The same method was used for each of the k values (k=1, ..., 45) and Fig. 6plots the results of all deviation sequences. According to Fig.6, Dmax and Dmin are 0.9815 and 0.0, respectively. These Dmax, Dmin and the distinguishing coeﬃcient n (set to 0.5) were substituted into Eq. 7 to yield the grey relational coeﬃcient. Figure7displays the results. Finally, the grey relational

coeﬃcients were averaged under equal weighting con-ditions using Eq. 8 to obtain the grey correlation grade. Table3 presents the grey correlation grade of the ref-erence sequence x

AVðkÞ and the comparability sequence

xi*(k), i = Dg, Dp, he, KS, a, b. Resorting the grey cor-relation grade in order of magnitude, as shown in the Table3, yields C (AV, a) > C (AV, Dg) > C (AV, he) > C (AV, Dp) > C (AV, KS) > C (AV, b). Accordingly, the order of importance of the factors that aﬀect the arrival time of the recharge water at the groundwater table was a, Dg, he, Dp, KS, b. Dg, he, Dpand KSare interpreted as the displacement of the recharge water, the capacity of the soil to store water, the potential for the movement of recharge water, and the ability of the soil material to transport water, respectively. Therefore, if a site can recharge to groundwater quickly, then the sites’ dis-placement of recharge water, the capacity the soil to

0.28 0.32 0.36 0.4 0.44 0.48 θ e D_g 10 m 20 m 30 m 0 0.5 1 1.5 2 2.5

K

s (m /day) 0.28 0.32 0.36 0.4 0.44 0.48 α 0 400 800 1200 1600

A

V

(Hr)

1 2 3 4 5 6 D_g 10 m 20 m 30 m D_g 10 m 20 m 30 m D_g 10 m 20 m 30 m β

(8)

store water, the potential for recharge water movement and the ability of the soil material to transport water must be evaluated in that order.

Figures2 and 4 show that a small Dg, a large he, a large Dp and a large KS promote the rapid recharge of water. Table3 indicates that a is the factor that most strongly affects the arrival time. A large a is typical of sands or well-structured soils and indicates a large change in water content with a small change of heads. Accordingly, the artificial pond water infiltrates easily through the soil with high a. Moreover, the result reveal that the hydraulic conductivity is not the most

important factor in determining the arrival time. Again, an incorrect decision may be made if an engineer sear-ches for a site for an artiﬁcial recharge pond based on only the magnitude of hydraulic conductivity.

4 Application

4.1 Study area

Yun-Lin is located in the southern part of the alluvial fan of the Cho-Shui River (Fig.8). The area is enclosed 0 0.5 1 X * D p 0 0.5 1 X * θ e 0 0.5 1 X * K S 0 0.5 1 X * α 0 5 10 15 20 25 30 35 40 45

Numerical run number

0 0.5 1 X * β 0 0.5 1 X * D g 0 0.5 1 X * A V

(9)

by the Taiwan Strait to the west, the Central Mountain to the east, the Cho-Shui River to the north, and the Pei-Kong River to the south. The area is approximately 1,000 km2, and extends 48 km from east to west and 24 km from north to south. Large gravels and cobble were deposited on the east side; near the coast, thick layers of clay and fine sand dominate. In the middle zone, the gravel layers are interbedded with sand and clay. WRA (1999) divided the groundwater basin of Yun-Lin area into layers stacked in the vertical direc-tion. In the coastal area, five layers of impermeable aquitard and four layers of aquifer are present. In the middle area, three layers of impermeable aquitard and four layers of aquifer are present. On the eastern side, only a single layer of the unconfined aquifer is present.

Figure9 displays the distribution of layers and their thickness.

The conventional source of revenue in the area is agricultural product. However, the average income of farmers in this region is considerably below that of workers in other sectors, because the industrial sector has recently seen rapid growth. Accordingly, many farmers have converted their farmlands into aquacul-tural operations and demand large amounts of water to boost earnings. Illegal pumping wells have been densely placed along the coastal area. The wells continuously extract groundwater from the aquifer to meet aquacul-tural demand and have caused seriously adverse eﬀects, such as land subsidence, seawater intrusion and envi-ronmental deterioration. Over the past two decades,

0 0.5 1 ∆ AV Dp 0 0.5 1 ∆ AV θ e 0 0.5 1 ∆ AV KS 0 0.5 1 ∆ AV α 0 5 10 15 20 25 30 35 40 45

0 0.5 1 ∆ AV β 0 0.5 1 ∆ AV Dg

(10)

around 610.4 km2 of the Yun-Lin plain has suﬀered from land subsidence. The cumulative land subsidence exceeded 215 cm from 1985 to 2002. Therefore, suitable sites must be selected for artiﬁcial recharge to reduce the rate of land subsidence (WRA 1998).

4.2 Data collection

Six hydrogeological factors, Dg, Dp, he, KS, a and b, are considered herein this work. Dp is assigned to 3 m. Other factors are obtained from the literature. The soil parameters are taken from Chen (2002), who measured the unsaturated hydraulic conductivity of 20 soil sam-ples over the Yun-Lin plain, as shown in Fig.8. They are ﬁtted using the Brooks–Corey model (Brooks and Corey 1964). The deﬁnition of parameters in the Mualem-van Genuchten model (van Genuchten1980) is not fully consistent with the parameters of the Brooks– Corey model, so parameters must be transformed between the two models. A parameter transformation developed by Hubert et al. (1996) is adopted and expressed as follows. 0 0.5 1

γ (

AV ,D p ) 0 0.5 1

γ

(AV ,θ e ) 0 0.5 1

γ

(AV ,K S ) 0 0.5 1

γ

(AV ,α ) 0 5 10 15 20 25 30 35 40 45

0 0.5 1

γ

(AV ,β ) 0 0.5 1

γ

(

AV ,D g )

Fig. 7 Calculated grey correlation coeﬃcients

Table 3 Grey correlation grade determined from the average of grey relational coeﬃcients with equal weightings, based on AVas a reference series Hydrogeological parameters a Dg he Dp KS b Grey correlation grade 0.7324 0.7114 0.6438 0.6143 0.5769 0.5747

(11)

11 b¼ 2 p 1; ð10Þ 1 a¼ hce 2pðp 1Þ pþ 3 55:6þ 7:4p þ p2 147:8þ 8:1p þ 0:092p2 ; ð11Þ

where hceis the entry pressure; k is the ﬁtting parameter in the Brooks–Corey model, and p=2/k + 3. At the same time, the depths of groundwater below the surface must be measured at the sampling sites. No observation well is present at the sampling site, so the values of Dgin the sampling sites were determined by interpolation from the values of Dg in the 42 observation wells that were constructed in Yun-Lin plain by WRA (2002). The Dgat an observation well was calculated by subtracting the daily average level form the level at the top of observation well. Table 4lists the values of Dg, Dp, he, KS, a and b in the 20 sampling sites.

4.3 Evaluation of artiﬁcial recharge site

An index must be defined for comparison among the sampling sites to identify an optimal site for artificial recharge from all of the sampling sites. An index, ID, was defined as follows;

ID¼X

m

i¼1

Wi FCi; ð12Þ

where Wiis the weight of factor i; FCiis the observed value of factor i, and m is the number of inﬂuential factors.

In this work, the depth of pond water was set to 3 m at all sampling sites. Hence, this factor has the same eﬀect on the movement of recharge water at all sampling sites. Accordingly, the number of factors m is reduced to

five. The grey correlation grade can be used to represent the grade of the factor’s contribution to the movement of the recharge water, so the grey correlation grade of the influence factor was taken as the weight in Eq. 12. The observed value yielded the normalized value of the influence factor. Equation 4 was used to determine the factors in the sequence he, KS, a, b and Dpbecause these factors are ‘‘the-larger-the-better’’. Equation 5 was used to evaluate the factor sequence Dg, because this factor is ‘‘the-smaller-the-better’’. Figure10 presents the result-ing indicator values at the samplresult-ing sites.

The index value supported the classiﬁcation of good sites, normal sites and poor sites. Good sites have index values from 2 to 3; normal sites have index values from 1 to 2; poor sites have index values from 0 to 1. Figure10 shows that sites S9, S12 and S13 are good sites; site S20 is a poor site, and the others are normal sites.

If the suitable area for artificial recharge is not the sampling site, then it can be identified by constructing a contour map of index values. Figure11plots a contour map of the distribution of index values, based on values at 20 sampling sites. The shading represents the value of index. As the gray color changes from dark to light, the index increases from small to large. A region with a large value indicates is suitable for artificial recharge. If a critical value of 1.6 is set which is about the mean value of the normal site, then the three areas A, B and C, displayed in Fig.11, are areas in which the index values exceed the critical value. Therefore, those areas are considered to be potentially good for artificial recharge in the Yun-Lin plain. Areas B and C are near the coast, and the recharge water from these two areas can only reach the upper aquifer. Area A is in the upper plain, and can directly recharge water to the lower layer of the aquifer. Area A is regarded as more suitable for groundwater recharge than area B and C.

S1 S2 S3 S4 S5 S6 S7 S9 _S11 S12 S20 S19 S18 S17 S13 S14 _S16 S15 S10 S8 Cho-Shui River Pei-Kong River 0 5Km N Samping site Taiwan Fig. 8 Geographical location of

Yun-Lin plain and sampling sites of the unsaturated soil experiment

(12)

5 Conclusion

Artiﬁcial recharge is a practical approach for increas-ing groundwater storage in an aquifer. The eﬃciency

of artificial recharge is related to site-specific hydro-geological factors. Therefore, the relationship between the hydrogeological factors and the movement of the recharge water front must be investigated to identify a suitable site for artificial recharge. A variable

satu-Elevation(m)

Taiwan Samping site N 5Km 0 Pei-K ong Rive r Cho-S hui Riv_er S8 S10 S15 S16 S14 S13 S17 S18 S19 _S20 S12 S11 S9 S7 S6 S5 S4 S3 S2 S1

Dou Liu Hill

Gravel mixed with sand layer Clay layer -300 -350 -250 -200 -150 -100 -50 200 150 100 50 0 250 300 350 400 450

(13)

rated groundwater model, FEMWATER, was used herein to simulate the wetting front of recharge water inﬁltration from the artiﬁcial recharge pond and to determine the arrival time of wetting front contacts with the groundwater table. The relationship between the arrival time and the hydrogeological factors used in FEMWATER was analyzed by the grey correlation method.

The six factors used in the FEMWATER were Dg, Dp, he, KS, a and b. Dg, Dp, he and KS represent the displacement of the recharge water, the capacity of the soil to store water, the potential for movement of the recharge water and the ability of the soil to transport water, respectively. Forty-five cases were simulated to obtain the corresponding arrival time. The grey corre-lation method was used to analyze the recorre-lationships between arrival time and the above six factors. The results indicate that the order of importance of the fac-tors to the arrival time, is in the order a, Dg, he, Dp, KS, b. The hydraulic conductivity is not the most important factor in determining the arrival time. An incorrect decision may be made if an engineer selects a site for artificial recharge based only on hydraulic conductivity. Additionally, the results provide a rule for selecting an artificial recharge site. The factors a, Dp, he, KS and b can be compared in that order to locate a suitable site for artificial recharge. High factor values indicate suit-ability for artificial recharge, except for factor Dg. Factors a and b are determined by fitting to an unsat-urated soil characteristic curve. However, evaluating a and b is time-consuming and expensive. If the project fund or time is limited, factors should be compared in order Dg, he, Dp, KSto select an optimal site for artificial recharge.

The method is also applied to determine a suitable site for artiﬁcial recharge in the Yun-Lin plain. The results of grey correlation analysis show that sampling sites S9, S12 and S13 are more suitable for artiﬁcial recharge than the others sites. Moreover, contour mapping of the distribution of index values reveals cases

0 0.5 1 1.5 2 2.5 3 S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 S12 S13 S14 S15 S16 S17 S18 S19 S20 Sampling site Index value

Fig. 10 Values of index at sampling site

Table 4 Soil parameters and depth of groundwater table below the surface at the sampling sites

Number of sampling site a b he KS(m/h) Dp(m) Dg(m) S1 0.0915 1.4500 0.3000 0.00201 3 18.75 S2 0.2479 1.1500 0.3500 0.04140 3 4.99 S3 0.6347 1.6200 0.2700 0.00445 3 8.79 S4 0.1228 1.3100 0.3700 0.00173 3 4.84 S5 0.5096 1.4200 0.3600 0.05044 3 6.16 S6 0.6498 1.4600 0.2500 0.02526 3 17.17 S7 0.3017 1.6400 0.3700 0.00401 3 9.39 S8 0.3449 1.3200 0.3000 0.00226 3 16.53 S9 0.6688 1.6300 0.3800 0.12577 3 3.13 S10 0.2194 1.4000 0.3800 0.00123 3 4.75 S11 0.1495 1.2800 0.3300 0.04765 3 12.37 S12 1.0036 1.9500 0.1800 0.02472 3 6.55 S13 0.8312 1.5100 0.4100 0.06488 3 20.01 S14 0.4311 2.0000 0.2500 0.01332 3 6.08 S15 0.3159 1.4200 0.2100 0.05682 3 21.19 S16 0.4909 1.6600 0.2200 0.02198 3 23.7 S17 0.3391 1.3300 0.2800 0.00255 3 16.76 S18 0.4587 1.4800 0.3200 0.01096 3 3.76 S19 0.4682 1.4400 0.4000 0.00725 3 5.93 S20 0.3354 1.5900 0.1300 0.01627 3 41.65

(14)

in which the area of artificial recharge falls outside the sampling sites. The results demonstrate that three areas are the suitable for artificial recharge. Two areas B and C are near the coast, where infiltrated water can only recharge to the upper aquifer, while area A is in the upper plain, where the infiltrated water can recharge directly to the lower aquifer. Area A is considered to be better suited to groundwater recharge than areas B and C in the Yun-Lin plain.

Acknowledgements We would like to thank the National Science Council of ROC (Taiwan) for ﬁnancially supporting this research under Contract No. NSC 91-2313-B-127-002.

References

Bouwer H (1996) Issue artiﬁcial recharge. Water Sci Tech 33:381– 390

Bouwer H (2002) Artiﬁcial recharge of groundwater: hydrogeology and engineering. Hydrogol J 10:121–142

Brooks RH, Corey AT (1964) Hydraulic properties of porous media affecting fluid flow. J Irrig Drain Div Am Soc Civ Eng 92(IR2):61–88

Chang NB (2005) Sustainable water resources management under uncertainty. Stoch Environ Res Risk Assess 19:97–98

Chang TC, Lin SJ (1999) Grey relation analysis of carbon dioxide emissions from industrial production. J Environ Manag 56:247–257

Chapman TG, Malone RW (2002) Comparison of models for estimation of groundwater recharge, using data from a deep weighing lysimeter. Math Comput Simul 59:3–17

Chen JF (2002) Analysis and application of a water budget model in the unsaturated zone. PhD Thesis, National Cheng Kung University, Tainan

Cook PG, Robinson NI (2002) Estimating groundwater recharge in fractured rock from environmental 3H and 36Cl, Clare Valley, South Australia. Water Resour Res 38. DOI 10.1029/ 2001WR000772

van Dam J C, Feddes RA (2000) Numerical simulation of inﬁl-tration, evaporation and shallow groundwater levels with the Richards equation. J Hydrol 233:72–85

Deng JL (1989) Introduction to the grey system theory. J Grey Syst 1:1–24 Deng JL (1996) Fundamental methods on grey system. Hua-Zhong

University of Science and Technology Press, Wu-Han El-Hames AS, Richards KS (1995) Testing the numerical diﬃculty applying

Richards’ equation to sandy and clayey soils. J Hydrol 167:381–391 EMRL (2003) Groundwater modeling system tutorials volume

IIV. Brigham Young University, UT

Eriksson E, Khunakasem V (1969) Chloride concentration in groundwater, recharge rate and rate of deposition of chloride in the Israel coastal plain. J Hydrol 7(2):178–197

Finch JW (1998) Estimating direct groundwater recharge using a simple water balance model—sensitivity to land surface parameters. J Hydrol 211:112–125

Fung CP (2003) Manufacturing process optimization for wear property of ﬁber reinforced polybutlene terephthalate com-posite with grey relation analysis. Wear 254:298–306

Gau HS, Chen JS, Peng JK, Hsieh SY (2001) The inﬂuence of hydraulic conductivity, pond water depth and ground water depth on artiﬁcial groundwater recharge (in Chinese). Taiwan Water Conserv l49:46–54

van Genuchten MT (1980) A close-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J 44:892–898 Giao PH, Phien-Wej N, Honjo Y (1999) FEM quasi-3D modeling of responses to artiﬁcial recharge in the Bangkok multiaquifer system. Environ Model Softw 14:141–151

Hodnett MG, Tomasella J (2002) Marked diﬀerences between van Genuchten soil water-retention parameters for temperate and tropical soils: a new water retention pedo-transfer functions developed for tropical soils. Geoderma 128:155–180

Hubert J, Meyer PD, Nachabe M, Touma J, van Genuchten MT, Lenhard RJ (1996) Parameter equivalence for the Brooks– Croey and van Genuchten soil characteristics: preserving the eﬀective capillary drive. Water Resour Res 32(5):1251–1258 Kao C, Bouarfa S, Zimmer D (2001) Steady state analysis of

unsaturated ﬂow above a shallow water-table aquifer drained by ditches. J Hydrol 250:122–133

Liang RH (1999) Application of grey relation analysis to hydroelectric generation scheduling. Electr Power Energy Syst 21:357–364 Lin HC, Richards DR, Yeh GT, Cheng JR, Cheng HP, Jones NL

(1997) FEMWATER: a three-dimensional ﬁnite element com-puter model for simulating density-dependent ﬂow and trans-port in variably saturated media, Technical Retrans-port CHL-97-12, U.S. Army Corps of Engineer, 151pp

Luckner L, van Genuchten MTh, Nielsen DR (1989) A consistent set of parametric models for the two phase ﬂow of immiscible ﬂuids in the subsurface. Water Resour Res 25:2187–2193

Cho-Shui River Pei-Kong Ri ver 0 5Km N 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Index value

A

B C

Fig. 11 Distribution of index values over Yun-Lin

(15)

Parlange JY, Hogarth WL, Barry DA, Parlange MB, Haverkamp R, Ross PJ, Steenhuis TS, DiCarlo DA, Katul G (1999) Ana-lytical approximation to the solutions of Richards’ equation with applications to inﬁltration, ponding, and time compression approximation. Adv Water Resour 23:189–194

Ragab R, Finch J, Harding R (1997) Estimation of groundwater recharge to chalk and sandstone aquifer using simple soil models. J Hydrol 190:19–41

Rastogi AK, Pandey SN (1998) Modeling of artificial recharge b of different shapes and effect on underlying system. J Hydrol Eng 3:62–68

Salama RB, Tapley I, Ishii T, Hawkes G (1994) Identiﬁcation of areas of recharge and disharge using Landsat-TM satellite imagery and aerial photography mapping techniques. J Hydrol 162:119–141

Saraf AK, Choudhury PR (1998) Integrated remote sensing and GIS for groundwater exploration and identiﬁcation of artiﬁcial recharge sites. Int J Remote Sensing 19:1825–1841

Sukhija BS, Reddy DV, Nandakumar MV, Rama (1997) A method for evaluation of artiﬁcial recharge through percolation tanks using environmental chloride. Ground Water 35:161–165 Todd DK (1980) Groundwater hydrology. Wiley, New York,

535 pp

Warrick AW (1991) Numerical approximations of darcian ﬂow through unsaturated soil. Water Resour Res 27:1215–1222 WRA (1998) Planning of groundwater artiﬁcial recharge on

Cho-Shui-River alluvial. Technical Report, Water Resources Agen-cy, Taiwan

WRA (1999) Corpus book for the reports of the entire ground-water monitoring plan-the ﬁrst stage from 1992 to 1998. Technical Report, Water Resources Agency, Taiwan

WRA (2002) Hydrological year book of Taiwan, Water Resources Agency, Taiwan

Yeh YL, Chen TC (2004) Application of grey correlation analysis for evaluating the artiﬁcial lake site in Pingtung Plain, Taiwan. Can J Environ Eng 31:56–64