MIKE 11 modeling system for an arid watershed

Journal of Hydrology (2007) 347, 448– 459

available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/jhydrol

Investigation of groundwater response to overland flow and topography using a coupled MIKE SHE/MIKE 11 modeling system for an arid watershed Hai-Long Liu a b c

a,b,*

, Xi Chen a, An-Ming Bao a, Ling Wang

c

Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences, Urumchi 830011, PR China Graduate School of the Chinese Academy of Sciences, Beijing 100039, PR China Normal College of Shihezi University, Shihezi 832000, PR China

Received 29 October 2006; received in revised form 7 September 2007; accepted 18 September 2007

KEYWORDS Groundwater recharge; Distributed hydrological model; Wavelet–fractal; Flooding in arid areas

The distribution of groundwater and dynamic fluctuations in groundwater levels have direct impacts on the eco-environment of arid areas. Investigations of groundwater recharge in arid areas are typically limited by a lack of adequate meteorological and hydrogeological records. This study focuses on groundwater recharge in a seasonally flooded arid area within the Tarim Basin, China, with the aim of analyzing the influence of groundwater and topography on the response characteristics of overland flow. We conducted a simulation using a coupled MIKE SHE/MIKE 11 model over 112 days of the flood season and calculated the average water balance. Based on the properties of the multiscale recognition of wavelets and the self-comparability of fractals in analyzing the detailed characteristics of groundwater diffusion and fluctuations in groundwater levels, a hybrid fractal–wavelet method was used to explain the recharge response associated with overland flow, distance from midstream, topography, and flooding depth. The results of the model simulations are generally consistent with observed data, indicating that the hybrid fractal–wavelet method is able to recognize the detailed characteristics of the groundwater response. Furthermore, the groundwater levels show a significant relationship with the orientation of the riverway prior to flooding. When flooding occurred, the groundwater levels showed a rapid response to changes in the depth of long-term overland flow. A total of 71.31% of the study area showed a strong correlation between groundwater level and the distance from midstream. The results demonstrate that the relationship exists for fluctuations in groundwater levels of more than 1.4 m. Variations in the height of the water table were significantly influenced by topographic elevation; in contrast, slope and aspect had little effect. In conclusion, the above results indicate that the pro-

Summary

* Corresponding author. Address: Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences, Urumchi 830011, PR China. Tel.: +86 991 7885373; fax: +86 991 7885378. E-mail address: [email protected] (H.-L. Liu). 0022-1694/$ - see front matter ª 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2007.09.053

Investigation of groundwater response to overland flow and topography using a coupled MIKE SHE/MIKE

449

posed methodology is applicable for the management of water resources in arid regions. The modeling and hybrid fractal–wavelet method study allowed quantification of the processes affecting groundwater levels and provided an insight into their implications in exploring groundwater level management. ª 2007 Elsevier B.V. All rights reserved.

Introduction In arid desert areas, insufficient rainfall and the vulnerability of the ecosystem mean that the growth of desert vegetation is closely related to the level of the water table (Fan et al., 2004; Tim et al., 2005). In flood-prone regions, riverbeds rise continuously due to sedimentation and flooding occurs in distinct flood seasons. For these reasons, plants in such areas are mainly distributed close to rivers and along the flood belts at the margins of alluvial fans (Wang et al., 2002). As groundwater plays an important role as the link between flooding and plant distribution, it is important to study the groundwater regime in order to conserve the eco-environment of arid areas. The Tarim basin, which includes the Taklimakan Desert, is home to nearly 10 million residents; however, the average groundwater levels in the Tarim River, the longest inland river in China, has fallen from 7 m in the 1970s to well over 12 m in the 1990s (Zhao et al., 2006; Gong et al., 2006). Extensive upstream irrigation means that the downstream stretches of the river have run dry since 1972 (Hao et al., 2006). Forests of poplar trees located along the river, formerly known as the ‘‘Green Corridor,’’ are vanishing day by day, and the falling water table is leading to desertification. Ecological restoration and conservation is therefore urgently required. Many groundwater studies have been carried out in arid and semi-arid areas around the world. Farah et al. (1997) assessed the groundwater resources in a semi-arid area by comparing and classifying lithological, hydrogeological, and hydrochemical characteristics. Urbano et al. (2000) studied groundwater–lake interactions in semi-arid environments, and Bouhlassa and Aiachi (2002) and Le Gal La Salle et al. (2001) used radiocarbon methods to investigate groundwater recharge under semi-arid conditions. de Marsily (2003) analyzed the interaction between groundwater and surface water and the behavior of temporary ponds in arid climates. Oren et al. (2004) investigated groundwater processes in an arid area and drilled several wells to compare natural groundwater with contaminated groundwater. The authors calculated the flow geometry by integrating several types of data (Navarro et al., 2005) and developed an effective approach for understanding groundwater dynamics in a semi-arid basin. Berendrecht et al. (2003) used a numerical simulation to investigate the movement and distribution of groundwater, while distributed hydrologic models have been used to understand the macro-scale movement mechanism of groundwater (Bogena et al., 2005). Time-series models (Box and Jenkins, 1976) and transfer function-noise (TFN) models (Hipel and Mcleod, 1994; Berendrecht et al., 2003) have both been extensively applied in analyzing groundwater dynamics.

In contrast to the above, few studies have investigated the temporal characteristics of groundwater in regions characterized by seasonal flooding. This work is essential to restore and conserve the ecology of arid areas such as those in western China. Through flood control, water consumption can be optimized, and ecological water demand can be supplied to the fullest extent. Furthermore, while a number of studies have investigated the dynamic characteristics of groundwater, there is currently insufficient detailed information on groundwater response. Wavelet and fractal theory is a new approach for analyzing fluctuating signals. Few such studies have been undertaken in the field of groundwater flow. Sahimi (2000) put forward a fractal–wavelet neural-network approach for analyzing characterization of field-scale porous media, such as oil, gas and geothermal reservoirs. Ruiz-Medina et al. (2003) solved the inverse estimation problem for random fields based on fractional-order regularization and wavelet approximation, The authors thought this method also can be used for the analysis of groundwater. The wavelet function is able to magnify local signals and recognize texture, and fractal theory is able to depict the general laws behind irregular phenomena; accordingly, we expect that these methods can be successfully applied in the field of groundwater flow. This paper aims to analyze the temporal characteristics of groundwater in a region characterized by seasonal flooding. The main objectives of the study are as follows: (1) couple a spatially distributed hydrological model for catchment hydrology and groundwater, implemented in the MIKE-SHE hydrological modeling software of DHI Water and Environment, with a MIKE 11 hydraulic model to simulate dynamic changes in groundwater within the study area for both flood and dry seasons; (2) analyze the characteristics of groundwater changes and their relations to topography, overland flow, distance to the river, etc., during both flood and fry seasons, using the wavelet transform and fractal approaches; (3) apply the proposed methodology to the Yingsu sub-watershed, a representative seasonally flooded tributary of the Tarim River.

Methodology MIKE SHE/MIKE 11 coupling MIKE SHE, a physically based distributed modeling system, is derived from the SHE hydrological model (Refsgaard and Storm, 1995). MIKE SHE was developed to model water movement, including overland flow, rivers and lakes, saturated and unsaturated flow, and evapotranspiration (Refsgaard, 1997). It can be used to describe the main physical processes of the hydrological cycle. In this model, overland

450

H.-L. Liu et al.

flow is modeled by solving the diffusive wave approximation in two horizontal directions, and vertical flow in unsaturated zones is modeled by a one-dimensional Richards equation. Flow in saturated zones is described by the three-dimensional Boussinesq equation. The methods proposed by Kristensen and Jensen (1975) are used to calculate evapotranspiration, including canopy interception, evaporation from the canopy, plant transpiration, and soil evaporation. MIKE 11 is a river modeling system based on the complete dynamic wave formulation of the Saint Venant equations. Through the coupling of the MIKE-SHE implementation for the watershed with the MIKE11 model of the river network in that watershed, the bi-directional interaction between the watershed hydrological processes and the river hydrodynamics can be accounted for. This allows the groundwater regimes and the interactions with the river bed flooding and the vegetations along the flooded areas in the arid watershed under study to be investigated. The problems involved in coupling MIKE SHE and MIKE 11 have already been solved (Sørensen et al., 1996; Refsgaard and Sørensen, 1997). In the coupled modeling system, simulations calculated using MIKE SHE and MIKE 11 are simultaneous, in which data exchange between the two models is realized through a shared storage space (Thompson et al., 2004; Sahoo et al., 2006). In the present paper, the water levels of the river and alluvial plain are calculated using MIKE 11, and the results are transferred to MIKE SHE. Depth and areal extent graphs of flooding can be compiled by comparing calculated water levels with topographic information stored in MIKE SHE. The residual water flow of the hydrological cycle is also calculated. Water exchange between the two models takes place via evapotranspiration, infiltration, overland runoff, and exchange between the stream and aquifer. Using source/sink terms of the Saint Venant continuous equations in MIKE 11, flow from MIKE SHE is exchanged with MIKE 11 (Hao and Wang, 2003). Coupling MIKE SHE and MIKE 11 (see Fig. 1) is essential in accurately describing the dynamic process of interaction between stream water and groundwater. This is because MIKE SHE is not appropriate for modeling exchange between

the stream and aquifer when the stream width is much smaller than the grid size calculated in MIKE SHE. Moreover, the complex branch-system-bearing loop and flood unit requires an expressional form of hydrodynamics that is as powerful as that of MIKE 11.

Hybrid fractal–wavelet (HFW) method In many physical and applied sciences, the method used to characterize irregular functions is always an issue, e.g., financial time-series, geologic shapes, medical time-series, interfaces that develop from non-equilibrium growth processes, and turbulent velocity signals (Arneodo et al., 1996). The hybrid fractal–wavelet (HFW) method (Podsiadlo and Stachowiak, 2003) is capable of handling such a problem. The detailed algorithm of the HFW is illustrated as follows. Assuming the signal fðtÞ 2 L2 ðRÞ (Cui, 1997; Qiu et al., 2003), the continuous wavelet transform is   Z 1 tb dt ð1Þ W f ða; bÞ 6 f; wa;b P jaj2 fðtÞw a R where t is a time variable, a is a scale parameter, and b is a translation parameter. Here, wa;b ðtÞ is a wavelet function that must fit: Z jwðxÞj2 dx < 1 ð2Þ Cw ¼ jxj R where x is frequency and:   1 tb wa;b ¼ pffiffiffiffiffiffi w a jaj

ð3Þ

Derived from the wavelet function wðtÞ, the wavelet wa;b ðtÞ serves as an observation window in the transformation. Consequently, wðtÞ also must fit: Z þ1 jwðtÞjdt < 1 ð4Þ 1

For a discrete wavelet (Xiong et al., 2001), the scale parameter a and translation parameter b need to be discretized. In Eq. (3), parameters a and b are given as a ¼ aj 0 , b ¼ kaj b , respectively, where j; k 2 Z; 0 < a < 1, and 0 0 0 b0 2 R. Therefore, the corresponding wavelet function wj;k is transformed as ! 1 t  kaj0 b0 ¼ aj=2 wðaj ð5Þ wj;k ðtÞ ¼ qffiffiffiffiffi w 0 0 t  kb0 Þ aj0 aj 0

The signal fðtÞ can be decomposed into multi-scale components interlacing each other by transforming the wavelet scale from small to large. Accordingly, every fine detail can be observed by adjusting corresponding sampling steps in the temporal or frequency domain. The typical characteristic of a fractal is its self-similarity at any scale. Assuming that fðtÞ is a real fractal function of t, fðt þ t0 Þ is the signal observed at the different scale k near t0 , and fðt0 þ ktÞ is the observed result (The Instrumental analysis section of department of chemistry of Peking University, 1997), we have: Figure 1

Coupling structure of MIKE 11 and MIKE SHE.

fðt þ t0 Þ ¼ kHðt0 Þ fðt0 þ ktÞ

ð6Þ

Investigation of groundwater response to overland flow and topography using a coupled MIKE SHE/MIKE where H is the Ho ¨lder or part-scale exponent. If each point of the function has the same fractal characteristic near the domain, then Eq. (4) can be extended to the entire course and is expressed as fðtÞ ¼ aH fðatÞ

ð7Þ

The dimensional number is used to describe the fractal as a quantitative parameter, and the capacity dimension Dc is derived from Kolmogorov (1958), based on covering. Supposing that the figure is a finite set in the n-dimensional Euclidean space R00 , the set is covered with a d-dimensional sphere with a radius of e, and NðeÞ is named as the minimum number of spheres, the capacity dimension Dc can be defined as ln NðeÞ Dc ¼ lim e!0 lnð1=eÞ

ð8Þ

The relationship between the capacity dimension Dc and the Ho ¨lder exponent (H) is Dc ¼ 2  H

ð9Þ

Decomposing the signal fðtÞ with the orthonormalized basis wavelet wðtÞ (Acharya et al., 2005), we have: X X ðmÞ fðtÞ ¼ d k wmk ðtÞ ð10Þ m

ðmÞ dk

¼

Z

k

fðtÞwmk ðtÞdt

ð11Þ

The following relationship exists between the Ho ¨lder exponent (H) and the wavelet coefficient for the fractal signal fðtÞ with self-comparability: ðmÞ

dk

m

ð0Þ

¼ b 2 d k ;

b ¼ 22Hþ1

ð12Þ

Based on the above discussion, we propose a HFW method. First, the water-flow field signal fðtÞ of every grid is transformed using Eq. (1) to obtain the discrete wavelet coefficient. Next, the fractal characteristics of the water-flow field are described by analyzing variations in the capacity dimension Dc ðtÞ over time t according to the relationship between the wavelet coefficients and the Ho ¨lder exponent, as in Eq. (12). This HFW method takes advantage of the strengths of the two components, i.e., the inherent ability of wavelets to decompose overland flow and groundwater data into individual scales and the inherent ability of fractals to describe overland flow and groundwater data over many scales in a scale-invariant manner (Arneodo et al., 1996). Groundwater flow is complicated because it is influenced by the characteristics of rock and soil, or water supply, and it is highly variable; however, when the scale of analysis is restricted to certain limits, the governing rules become apparent. Accordingly, we consider that the HFW method is appropriate in the analysis of groundwater flow.

Pn

j¼1 ðO

 PÞ2

j¼1 ðO

 OÞ2

E ¼ 1  Pn

451 ð13Þ

where O and O are the observed water table and the mean observed water table, respectively, and P is the predicted water table. If the predictions equal the mean observed water table, then E = 0. When all the predictions equal their corresponding observations, then E = 1.

Case study Description of the study area The study area, the Yingsu sub-watershed in the middle reaches of the Tarim River, is located within the arid Tarim Basin and covers a surface area of about 91.16 km2 between latitudes 413 0 2200 –4059 0 3400 N and longitudes 8629 0 5100 – 8640 0 1300 E (Fig. 2). The climate in this area is extremely dry, with only rare rainfall. The average annual precipitation is 33.1 mm, and the average annual temperatures is 10.7 C, with maximum temperatures of around 43 C and minimum temperatures of around 31 C; annual sunlight hours are approximately 3000. The area consists of alluvial plain and a transition zone to wind-accumulated dunes. The soil is predominantly silt loam, fine sand, and sandy loam. Groundwater recharge mainly occurs by lateral infiltration from stream flow and flooding during the flood season, whereas during dry periods groundwater discharge includes evaporation and plant transpiration. The landforms of the area are characterized by low relief between elevations of 866.4 and 927.8 m, within which flooding readily occurs under large flows because of the narrow width of the riverway.

Available data We obtained data on daily temperature, precipitation, and potential evapotranspiration for the period June 1 to September 20, 2002, from the Xinjiang Autonomous Region Climate Center. The size of the grid cells for the model of the Yingsu sub-watershed was 50 · 50 m. The data that describe soil layers and the hydraulic characteristics of the soil were acquired from the Design Institute of the Yellow River Water Conservancy, while streamflow and groundwater level data were obtained from the Tarim River Basin Administrative Bureau. Land-use types and Leaf Area Index (LAI) data were calculated based on Normalized Difference Vegetation Index (NDVI) from Landsat thematic mapper data (TM) and Moderate Resolution Imaging Spectroradiometer data (MODIS) for the corresponding period. Stream profile properties, such as cross section area and storage width, were obtained from Digital Elevation Model data (DEM). TM data were also used to regulate the river location.

Model uncertainty estimation

Calibration parameters

To statistically evaluate model performance, a standard performance metric, the coefficient of efficiency (Nash and Sutcliffe, 1970; McMichael et al., 2006), is generally used to evaluate model performance following each run:

The main parameter sets used in the modeling system are listed in Table 1. The actual evapotranspiration was calculated using the Kristensen–Jensen model (Kristensen and Jensen, 1975):

452

H.-L. Liu et al.

Figure 2 Table 1

Relief map of the studied watershed.

Calibrated parameter values

Parameter

Calibrated value

C1 C2 C3 Kv (m s1) Kh (m s1) Ks (m s1) n (s m1/3)

0.3 0.2 20 1.2 · 105 5 · 103 1.2 · 108 0.017

E at ¼ f1 ðLAIÞ  f2 ðhÞ  RDF  E p

on the soil moisture content in the root zone, and RDF is a root distribution function, and Ep is the potential evapotranspiration rate. f1 ðLAIÞ ¼ C2 þ C1 LAI

ð15Þ

where C1 and C2 are empirical parameters. 

hFC  h f2 ðhÞ ¼ 1  hFC  hW

ð14Þ

where Eat is the actual evapotranspiration, f1(LAI) is a function based on the leaf area index, f2 ðhÞ is a function based

CE 3

p

ð16Þ

where hFC is the volumetric moisture content at field capacity, hW is the volumetric moisture content at the wilting point, h is the actual volumetric moisture content and C3 is an empirical parameter. RDF ¼

Z

z2

z1

eðln R0 AROOTzÞ dz

Z 0

L

eðln R0 AROOTzÞ dz

ð17Þ

Investigation of groundwater response to overland flow and topography using a coupled MIKE SHE/MIKE

453

where R0 is the root extraction at the soil surface, AROOT is a parameter that describes the root mass distribution, z, z1 and z2 are the depths below ground surface, and L is the maximum root depth. The parameters C1 and C2 control the distribution of actual evapotranspiration between transpiration and soil evaporation, while C3 influences the value of the moisture content function (f2). We calibrated the parameters by comparing simulated evapotranspiration with that estimated using the Penman formula for the period from June 1 to September 20, 2002. The values of vertical saturated hydraulic conductivity (Kv), horizontal saturated hydraulic conductivity (Kh), saturated hydraulic conductivity (Ks), and Manning coefficient (n) of overland flow used in the calibrated model are given in Table 1. The initial values were derived from a preliminary design report for the Yingsu ecological sluice gate published by the Design Institute of the Yellow River Conservancy in July 2003 (Xu et al., 2003). The values listed in Table 1 are the final calibration parameters values.

Results and discussion Model calibration, testing, and estimate of uncertainties The time scale and amount of data used by calibration and testing depend on the physical process and period of change (Zhang et al., 2007; Dobrzan ´ski et al., 2007; Tsai and Li, 2008). The final values of the calibration parameters in Table 1 are consistent with the values proposed in previous studies conducted within arid and semi-arid areas (Christine et al., 2006; Zhao et al., 2005). Fig. 3 shows a comparison of the simulated and measured daily stream flow at the outlet of the watershed for a period of 28 days (June 1–28). Fig. 4 shows a comparison of simulated and measured groundwater depth over the same period. The visual comparison indicates that the model performance is satisfactory. Parameter sets (see Table 1 for the main initial parameters) were used to conduct the Yingsu simulations for a calibration period of 28 days (June 1–28, 2002). Values of E ranged from 0.80 to 0.89 (Table 2) for the groundwater depths at the A, B, C, D, E, F, G, H, I and J observation sites (see Fig. 2). The threshold value of E applied in our application of MIKE SHE follows that used by Andersen et al. (2001) and Christine et al. (2006), who found good model runs for E P 0.80. The results indicate that the model is acceptable

Figure 3

Observed and simulated stream flow.

Figure 4 Observed and simulated depth of groundwater at four locations within the Yingsu watershed.

for conducting a simulation of the Yingsu sub-watershed. Overland flow depth has not been measured in this study for it is affected by surface roughness and changes quickly. This makes accurate observation difficult. So it cannot be tested. We calculated prediction errors for daily data over the test period (June 29 to July 13, 2002). Approximately 75% of the observed water table values fell within the 5% and 95% uncertainty bounds based on regression of simulated versus observed values. Compared with studies of semi-arid shrublands in California, the accuracy of the model predictions is considered to be acceptable (Christine et al., 2006).

454 Table 2

H.-L. Liu et al. Estimate of uncertainties in selected observation locations

Observation site

Calibration period

Ecalibration

Test period

Etest

Date of flooding occurring

Date of field observations

A B C D E F G H I J

June 1–June 28

0.86 0.84 0.83 0.89 0.87 0.84 0.80 0.82 0.83 0.80

June 29–July 13

0.83 0.82 0.81 0.85 0.85 0.83 0.81 0.80 0.84 0.81

June 12–June 15

June 12–June 15

In the following analysis, the MIKE SHE/MIKE 11 model is used to simulate overland flow and groundwater distribution for the period from June 1 to September 20, 2006.

Characteristics of changes in groundwater patterns Flow in the Yingsu stream, which is a branch of the Tarim River, is recorded when flow in the Tarim river exceeds a certain level (Xuan and Wang, 2004). In 2002, flow was recorded in the Yingsu stream from June 1 to September 20 (i.e., the flood season). The characteristics of variations in groundwater during this period are described below.

Diffusion prior to flooding The simulation results indicate that no flooding occurs before June 12. This is consistent with field observations. Fig. 5 shows temporal fluctuations in groundwater levels on June 4, 7, and 12 compared with those on June 1, 4, and 7, respectively. The groundwater level increases along the direction of streamflow from west to east. The rate of increase in groundwater level is much faster in streamside areas than in remote areas prior to flooding. This pattern occurs because the groundwater is mainly recharged by lateral infiltration from stream flow during this period. Moreover, the rise in groundwater level is less significant (Zhang and Schilling, 2006), and the fluctuation range is

Figure 5 Changes in groundwater level prior to flooding. (a) Changes in groundwater level (m) from June 1–4. (b) Changes in groundwater level (m) from June 4–7. (C) Changes in groundwater level (m) from June 7–12.

Investigation of groundwater response to overland flow and topography using a coupled MIKE SHE/MIKE 1–1.5 m in most areas, which is closely related to the inflow of the stream during this period. To analyze the amount of groundwater recharge caused by infiltration at different orientations, we consider the fluctuation gradient, which is defined as the ratio of the fluctuation value of the water table to the distance from midstream (Land and Paull, 2001). This ratio indicates the influence of the distance from midstream on groundwater recharge at different sites. First, we calculated the varying quantity of phreatic fluctuation for each grid point over the period June 1–12 according to the corresponding graph of groundwater distribution. We then calculated the shortest distance from each grid point to the riverway (distance from midstream). This enabled us to generate a graph of the distribution of the fluctuation gradient. As shown in Fig. 6, there are directional differences in spatial variations in the mean fluctuation gradient. The results indicate that groundwater levels alternately increased and decreased in each direction, and the sensitivity of the values to the distance from midstream alternated between strong and weak; this is in accordance with the zonal distribution rule for groundwater. However, the fluctuation gradient shows an increasing trend from north to south, while there is little variation from east to west. These trends suggest that the rate of groundwater recharge is higher in the southern area than in the north; there is no such trend from east to west. Comparison of the mean values of the fluctuation gradient in the south–north direction (0.0067) with that in the east–west direction (0.0071) reveals that spatial variations in soil transmission capacity are larger in the east–west direction than in the south– north direction. The difference in fluctuation gradient in south–north direction compared with east–west direction is significant at the 1% level. Significance testing proves that spatial variations do exist. Relationship between fluctuations in groundwater level and depth of overland flow When overland flow occurs, groundwater is recharged at both streamsides and through vertical infiltration (Wallach et al., 1997). The vertical infiltration rate is closely related to the timing of overland flow and its depth, as well as the properties of the soil profile. We overlaid all the two-dimensional (2D) maps (221 · 165 pixels, 99 maps) of groundwater level and those of overland

455

flow depth in a temporal sequence for the period from June 13 to September 20, producing two corresponding series of arrays. Subsequently, we decomposed the two groups of 2D signals using a compactly supported and double-orthogonal wavelet function to filter out detailed components of the fluctuation signal (Acharya et al., 2005). We also obtained two groups of decomposed signal series (1 row · 3,610,035 columns). Finally, we analyzed the temporal response relationship between a variety of groundwater levels and the depth of overland flow. The results of this analysis are shown in Table 3. Layer in Table 3 is the decomposed times of signals. The results indicate that the depth of overland flow is strongly related to fluctuations in the groundwater level. The more that the layers of signals are decomposed, the stronger the correlations (Table 3). The next layer means that the signals would be decomposed on this time. Part of high-frequency signals (i.e. short-period signals) would be filtered out (Cui, 1997). The decomposed result reveals that correlation between the depth of overland flow and degree of fluctuation in groundwater level is stronger for longer analysis periods. Therefore, there is a long-period time–response relationship between these two factors. The hybrid fractal–wavelet approach is used to analyze the space–response relationship between groundwater fluctuation and depth of overland flow (Acharya et al., 2005; Arneodo et al., 1996). Details of this manipulation are described in Section ‘‘Hybrid fractal–wavelet method’’. First, mean values of daily overland flow depth and groundwater level are calculated for the duration of the study period. These values are then further transformed by wavelet transformation and fractal dimension. According to the comparison of spatial variations in fractal dimensions between overland flow depth and groundwater level (Fig. 7), the fluctuation range of fractal dimensions is 1.62–2.63 for overland flow depth and 1.64–2.69 for groundwater level. It is apparent that overland flow depth records less spatial variation than groundwater levels. The variance of two-dimensional density maps is useful to estimate the confidence interval of a map and information on local characteristics (Cheng and Yeager, 2007). The variance from nearest neighbors(20 · 20 pixels) has been used to estimate the fluctuation range. The pixel numbers of

Table 3 Relationship between variations in groundwater level and depth of overland flow

Figure 6 Spatial variations in the fluctuation gradient for groundwater level.

Layer

Correlation coefficient

0 1 2 3 4 5 6 7 8 9 10 11

0.1625 0.5137 0.5133 0.5132 0.4848 0.4785 0.5816 0.6952 0.7997 0.8801 0.9332 0.9645

456

H.-L. Liu et al.

Figure 7

Spatial variations in the fractal dimensions of overland flow depth (a) and groundwater level (b).

overland flow depth with value more than 2.0 are 2240, and those of groundwater levels are 11,887. They occupy 6.14% and 32.60% of the whole pixel numbers, respectively. The mean difference is significant at the 1% level by t-test. The reasons for this are as follows: (1) overland flow changed markedly before June 20, and this required the filling of pocket storage across the entire sub-watershed, whereas after June 20 the level changed slowly; (2) in contrast, the process of groundwater recharge continued throughout the entire simulation period because the water table is very low in the arid area. Both fractal dimension maps are classified into four classes (high, moderate, low, and no change) based on the KMeans classification approach and a minimum distance technique (Smet and Guzma ´n, 2004). The overlay analysis approach was used to analyze the corresponding response region of fluctuation intensity (Shi et al., 2004). The results indicate that the response areas for high, moderate, low and no change occupy 14.47%, 27.75%, 29.09%, and 28.69% of the study area, respectively. The areas with a high response relationship mainly occur in the upper reaches and areas with little relief in the lower reaches. Areas of no response mainly occur at sites with high topography (Fig. 8). In addition, correlation analysis of mean values in east–west and south–north directions demonstrate that the correlation coefficients of overland flow depth and groundwater level in south–north and east–west directions are 0.7475 and 0.8282, respectively. These response characteristics

are consistent with the variations in soil transmission capacity in different directions described above. Relationship between fluctuations in groundwater level and distance from midstream We now consider the influence of distance from midstream on groundwater distribution. First, we calculated variations in the depth to the phreatic surface for every grid point from June 1 to September 20. The shortest distance to the midstream for every grid point was estimated using the buffer analysis method (Xiang, 1996). Second, regions of interest (ROI) (Bradley and Stentiford, 2003) were demarcated based on the degree of variation in groundwater level. The distance map was then clipped according to the ROI and the shortest-distance map, and we obtained the map showing the distance from midstream, which corresponds to the map showing change in depth to the phreatic surface. Finally, we estimated mean distances and areal percentage for each ROI. Table 4 shows the relationship between fluctuations in groundwater level and distances from midstream in terms of percentages and mean values. The ranges of fluctuations are mainly within 1.4 m, and the area within this scope occupies approximately 89.4% of the total area. There is no synchronized response between these two factors; however, the extent of fluctuation increases with increasing average distance when the fluctuation in the water table is greater than 1.4 m. The reason for this finding is that

Table 4 Relationship between fluctuations in groundwater level and distance from midstream

Figure 8 Intensity of groundwater level response to overland flow depth.

Groundwater level (m)

Average distance (m)

Percentage (%)

0–0.2 0.2–0.4 0.4–0.6 0.6–0.8 0.8–1.0 1.0–1.2 1.2–1.4 1.4–1.6 1.6–1.8 >1.8

1799.39 2071.28 2427.62 2062.49 2379.56 2231.29 1729.42 2520.49 2949.78 3360.87

0.41 11.63 28.84 13.68 10.56 10.10 14.18 5.82 4.54 0.24

Investigation of groundwater response to overland flow and topography using a coupled MIKE SHE/MIKE

Figure 9

457

Comparison of the classification of sub-watershed elevation (a) and change in groundwater level (b).

the distance is not a major impact factor for those areas with high initial groundwater levels, which are areas close to the river channel where only minor fluctuations are recorded. In contrast, for areas far from the channel, the initial groundwater level is lower and greater fluctuations in groundwater levels are recorded. Kebede et al. (2005) undertook a similar analysis for the Blue Nile River using geochemical and environmental isotope data. They found that the distance from midstream is the dominant factor. Relationship between fluctuations in groundwater level and topography Lin et al. (2005) found that topography plays an important role in controlling variations in soil morphology and groundwater distribution. Our comparison of the classification of elevation and fluctuations in groundwater level from June 1 to September 20 indicate a significant correlation (r = 0.6533, significant at 1% level), especially along the channel (Fig. 9). To estimate the impact of topographic factors on fluctuations in groundwater level, various values of groundwater level, elevation, slope, and aspect were decomposed by a near-symmetrical compactly supported double-orthogonal wavelet (Beta et al., 2003). All four factors were decomposed into detailed horizontal and vertical components that correspond to detailed characteristics in east–west and south–north directions, respectively. We then examined the relationship between the detailed components. The results indicate that groundwater recharge is strongly influenced by elevation, and only weakly influenced by slope and aspect (Table 5). Groundwater levels increase with increasing elevation, which is consistent with the distribution map of Michael et al. (2000). This trend reflects the fact that the influence of various flow-driving mechanisms, such as pressure distribution, changes with depth.

Table 5 Correlations between various amplitudes of groundwater levels and detailed components of topographic factors Correlation coefficient

Elevation

Slope

Aspect

East–west direction South–north direction

0.1954 0.8634

0.1928 0.5742

0.1980 0.1343

Furthermore, the impact of elevation is more significant in a south–north direction than in an east–west direction (Table 5). This is because the riverway flows toward east and topography varies only gradually in this direction. Initial groundwater levels vary markedly in a north–south direction. By comparison of the modeling maps, the groundwater is recharged intensively in that direction during the whole course. On the basis of the above analysis, we conclude that topography has an important impact on groundwater recharge and distribution.

Conclusions The coupled MIKE SHE/MIKE 11 model has been applied extensively (French and Clifford, 2000), and the performance of the model has been assessed using the coefficient of efficiency E (Nash and Sutcliffe, 1970; McMichael et al., 2006). This coefficient is greater than 0.80 for groundwater depths in all 10 observation sites of the present study and approximately 75% of the observed water table values were within the 5% and 95% uncertainty bounds of the regression of simulated versus observed groundwater depths. This result demonstrates that the model is acceptable for the simulation of overland flow and groundwater in the study area. Before flooding occurs, fluctuations in groundwater range from 1 to 1.5 m in most areas. A quantitative analysis of fluctuation gradients indicates that spatial variations in soil transmission capacity are greater in a stream-parallel direction than in a vertical direction. During flooding, there is a long-term spatial response relationship of the response between overland flow and groundwater fluctuation. Distance is not a dominant factor in areas close to the river channel with high initial groundwater levels. In contrast, for areas far from the river, the distance from midstream is the dominant factor (Fan et al., 2004). The distance from midstream shows a positive response to the amount of groundwater recharge at sites far from the river (further than 2.5 km); there is no such linear relationship for sites close to the river. In general, groundwater recharge is affected by both riverway direction and orientation, as well as topographic factors such as elevation. Based on the results of previous work (Acharya et al., 2005; Arneodo et al., 1996), we successfully used the hybrid

458 fractal–wavelet approach to analyze the temporal characteristics of groundwater. This approach can also be used in the future to analyze flow fields and the diffusion of contamination. In addition, the fluctuation gradient is useful for assessing the intensity of fluctuations. This study describes an efficient method for analyzing groundwater dynamics and their response relationship to environmental factors. This is important for considering groundwater pollution and ecological risk in arid areas.

Acknowledgements Funding for this research was provided in part by a National Science Foundation of China (No. 40571030) and the Tarim River Basin Management Information Systems Project in China, supported by the World Bank. The authors are grateful to Drs. Gordon Huang and Patrick Willems for their insightful comments on this manuscript. We wish to thank Ms. Yongqin Wang, an engineer in the Tarim River Basin Administrative Bureau, for her help in organizing the groundwater data. The authors also appreciate other members of the Tarim River Basin Management Project for their assistance.

References Andersen, J., Refsgaard, J.C., Jensen, K.H., 2001. Distributed hydrological modelling of the Senegal River Catchment-model construction and validation. Journal of Hydrology 247, 200–214. Acharya, U.R., Bhat, P.S., Kannathal, N., Rao, A., Lim, C.M., 2005. Analysis of cardiac health using fractal dimension and wavelet transformation. ITBM-RBM 26, 133–139. Arneodo, A., Aubenton, Y.D., Bacry, E., Manneville, S., Muzy, J.F., Roux, S.G., 1996. Wavelet based fractal analysis of DNA sequences. Physica D: Nonlinear Phenomena 96, 291–320. Berendrecht, W.L., Heemink, A.W., Geer, F.C.V., Gehrels, J.C., 2003. Decoupling of modeling and measuring interval in groundwater time series analysis based on response characteristics. Journal of Hydrology 278, 1–16. Beta, C., Schneider, K., Farge, M., Bockhorn, H., 2003. Numerical study of mixing of passive and reactive scalars in two-dimensional turbulent flows using orthogonal wavelet filtering. Chemical Engineering Science 58, 1463–1477. Bogena, H., Kunkel, R., Scho ¨bel, T., Schrey, H.P., Wendland, F., 2005. Distributed modeling of groundwater recharge at the macroscale. Ecological Modelling 187, 15–26. Bouhlassa, S., Aiachi, A., 2002. Groundwater dating with radiocarbon: application to an aquifer under semi-arid conditions in the south of Morocco (Guelmime). Applied Radiation and Isotopes 56, 637–647. Box, G.E.P., Jenkins, G.M., 1976. Time series analysis: forecasting and control. Holden Day, San Francisco. Bradley, B.P., Stentiford, F.W.M., 2003. Visual attention for region of interest coding in JPEG 2000. Journal of Visual Communication and Image Representation 14, 232–250. Cheng, A., Yeager, M., 2007. Bootstrap resampling for voxel-wise variance analysis of three-dimensional density maps derived by image analysis of two-dimensional crystals. Journal of Structural Biology 158, 19–32. Christine, E.M., Allen, S.H., Hugo, A.L., 2006. Distributed hydrological modeling in California semi-arid shrublands: MIKE SHE model calibration and uncertainty estimation. Journal of Hydrology 317, 307–324.

H.-L. Liu et al. Cui, J.T., 1997. Theories on Wavelet Analysis. Xi’An jiaotong Universities, Xi’an, pp. 47–71. de Marsily, G., 2003. Importance of the maintenance of temporary ponds in arid climates for the recharge of groundwater. Comptes Rendus Geosciences 335, 933–934. Dobrzan ´ski, L.A., Maniara, R., Sokolowski, J., Kasprzak, W., Krupinski, M., Brytan, Z., 2007. Applications of the artificial intelligence methods for modeling of the ACAlSi7Cu alloy crystallization process. Journal of Materials Processing Technology 192, 582–587. Fan, Z.L., Ma, Y.J., Zhang, H., 2004. Research of eco-water table and rational depth of groundwater of Tarim River Drainage Basin. Arid Land Geography 27, 8–13. Farah, E.A., Abdullatif, O.M., Kheir, O.M., Barazi, N., 1997. Groundwater resources in a semi-arid area: a case study from central Sudan. Journal of African Earth Sciences 25, 453–466. French, J.R., Clifford, N.J., 2000. Hydrodynamic modelling as a basis for explaining estuarine environmental dynamics: some computational and methodological issues. Hydrological Processes 14, 2089–2108. Gong, L., Wang, H.W., Bao, P.Y., Lv, G.H., Pan, X.L., 2006. Landscape change and ecological effect of typical oasis in the upper reaches of Tarim River. Journal of Desert Research 26, 421–425. Hao, X.M., Chen, Y.N., Li, W.H., 2006. The driving forces of environmental change during the last 50 years in the Tarim River Basin. Acta Geographica Sinica 61, 262–272. Hao, F.H., Wang, L., 2003. Distributed hydrological modeling. Yellow river conservancy press, Zhenzhou, China, pp. 138–144. Hipel, K.W., Mcleod, A.I., 1994. Time series modelling of water resources and environmental systems. Developments in Water Science. Elsevier, New York, pp. 45. Kebede, S., Travi, Y., Alemayehu, T., Ayenew, T., 2005. Groundwater recharge, circulation and geochemical evolution in the source region of the Blue Nile River, Ethiopia. Applied Geochemistry 20, 1658–1676. Kolmogorov, A.N., 1958. A new metric invariant of transitive dynamic system and automorphysms in Lebesgue spaces. Dokl Akad Nauk SSSR 119, 861–864. Kristensen, K.J., Jensen, S.E., 1975. A model of estimating actual evapotranspiration from potential evapotranspiration. Nordic Hydrology 6, 170–188. Land, L.A., Paull, C.K., 2001. Thermal gradients as a tool for estimating groundwater advective rates in a coastal estuary: White Oak River, North Carolina, USA. Journal of Hydrology 248, 198–215. Le Gal La Salle, C., Marlin, C., Leduc, C., Taupin, J.D., Massault, M., Favreau, G., 2001. Renewal rate estimation of groundwater based on radioactive tracers (3H, 14C) in an unconfined aquifer in a semi-arid area, Iullemeden Basin, Niger. Journal of Hydrology 254, 145–156. Lin, Y.S., Chen, Y.G., Chen, Z.S., Hsieh, M.L., 2005. Soil morphological variations on the Taoyuan Terrace, Northwestern Taiwan: roles of topography and groundwater. Geomorphology 69, 138–151. McMichael, C.E., Hope, A.S., Loaiciga, H., 2006. Distributed hydrological modelling in California semi-arid shrublands: MIKE SHE model calibration and uncertainty estimation. Journal of Hydrology 317, 307–324. Michael, K., Bachu, S., Machel, H.G., 2000. Groundwater flow in response to ground surface topography, erosional rebound, and hydrocarbon generation in Cretaceous strata in the Alberta Basin, Canada. Journal of Geochemical Exploration 69, 657–661. Nash, I.E., Sutcliffe, I.V., 1970. River flow forecasting through conceptual models. Journal of Hydrology 10, 282–290. Navarro, I., Soliz, J.G., Mahlknecht, J., 2005. Groundwater flow regime under natural conditions as inferred from past evidence and contemporary field observations in a semi-arid basin:

Investigation of groundwater response to overland flow and topography using a coupled MIKE SHE/MIKE Cuenca de la Independencia, Guanajuato, Me ´xico. Journal of Arid Environments 63, 756–771. Oren, O., Yechieli, Y., Bo ¨hlke, J.K., Dody, A., 2004. Contamination of groundwater under cultivated fields in an arid environment, central Arava Valley, Israel. Journal of Hydrology 290, 312–328. Podsiadlo, P., Stachowiak, G.W., 2003. Fractal–wavelet based classification of tribological surfaces. Wear 254, 1189–1198. Qiu, J.D., Liang, R.P., Zou, X.Y., 2003. The resolution of overlapping peaks by utilizing wavelet fractal. Journal of Instrumental Analysis 22, 1–5. Refsgaard, J.C., 1997. Parameterization, calibration and validation of distributed hydrological models. Journal of Hydrology 198, 69–97. Refsgaard, J.C., Sørensen, H.R., 1997. Water management of the Gabcikovo Scheme for balancing the interest of hydropower and environment. In: Refsgaard, J.C., Karalis, E.A. (Eds.), Operational Water Management. Balkema, Rotterdam, The Netherlands, pp. 365–372. Refsgaard, J.C., Storm, B., 1995. Computer Models of Watershed Hydrology. Water Resources Publications, Englewood, USA, pp. 809–846. Ruiz-Medina, M.D., Angulo, J.M., Anh, V.V., 2003. Fractional-order regularization and wavelet approximation to the inverse estimation problem for random fields. Journal of Multivariate Analysis 85, 192–216. Sahimi, M., 2000. Fractal–wavelet neural-network approach to characterization and upscaling of fractured reservoirs. Computers & Geosciences 26, 877–905. Sahoo, G.B., Ray, C., Carlo, E.H.D., 2006. Calibration and validation of a physically distributed hydrological model, MIKE SHE, to predict streamflow at high frequency in a flashy mountainous Hawaii stream. Journal of Hydrology 327, 94–109. Shi, W.Z., Cheung, C.K., Tong, X.H., 2004. Modelling error propagation in vector-based overlay analysis. Journal of Photogrammetry and Remote Sensing 59, 47–59. Smet, Y.D., Guzma ´n, L.M., 2004. Towards multicriteria clustering: an extension of the k-means algorithm. European Journal of Operational Research 158, 390–398. Sørensen, H.R., Klucovska, J., Topolska, J., Clausen, T., Refsgaard, J.C., 1996. An engineering case study-modelling the influences of Gabcikovo hydropower plant on the hydrology and ecology in the Slovakian part of the river branch system of Zitny Ostrov. In: Refsgaard, J.C. (Ed.), Distributed Hydrological Modelling. Kluwer, Dordrecht, The Netherlands, pp. 233–253. The instrumental analysis section of department of chemistry of Peking University, 1997. Instrumental Analysis Tutorial. Peking University Press, Beijing, pp. 304–306.

459

Thompson, J.R., Sørenson, H.R., Gavin, H., Refsgaard, A., 2004. Application of the coupled MIKE SHE/MIKE 11 modelling system to a lowland wet grassland in southeast England. Journal of Hydrology 293, 151–179. Tim, E., Tom, H., Ian, N., 2005. An ecological optimality approach for predicting deep drainage from tree belts of alley farms in water–limited environments. Agricultural Water Management 75, 92–116. Tsai, T.I., Li, D.C., 2008. Approximate modeling for high order nonlinear functions using small sample sets. Expert Systems with Applications 34, 564–569. Urbano, L.D., Person, M., Hanor, J., 2000. Groundwater–lake interactions in semi-arid environments. Journal of Geochemical Exploration 69, 423–427. Wallach, R., Grigorin, G., Rivlin, J., 1997. The errors in surface runoff prediction by neglecting the relationship between infiltration rate and overland flow depth. Journal of Hydrology 200, 243–259. Wang, Y.G., Hu, C.H., Zhou, W.H., 2002. Study on river patterns of the Tarim River. Journal of Sediment Research 6, 19–25. Xiang, W.L., 1996. GIS-based riparian buffer analysis: injecting geographic information into landscape planning. Landscape and Urban Planning 34, 1–10. Xiong, F., Jiang, X.Q., Gao, Y.S., Li, Z., 2001. Evaluation of engineering surfaces using a combined fractal modeling and wavelet analysis method. International Journal of Machine Tools and Manufacture 41, 2187–2193. Xu, R., Wu, C.Z., Lu, X.J., and Wang, Y.C., 2003. The preliminary design reports on ecological gate and weir project of Tarim River. Zhenzhou, Design Institute of the Yellow River Conservancy, pp. 49–69. Xuan, B.D., Wang, L., 2004. Design practice of ecological sluices on middle and lower Tarim River. Water and Electricity in Northwest 2, 10–11. Zhang, Y.K., Schilling, K.E., 2006. Effects of land cover on water table, soil moisture, evapotranspiration, and groundwater recharge: a field observation and analysis. Journal of Hydrology 319, 328–338. Zhang, X.G., Guo, Y.C., Chan, C.K., Lin, W.Y., 2007. Numerical simulations on fire spread and smoke movement in an underground car park. Building and Environment 42, 3466–3475. Zhao, C.Y., Wang, Y.H., Li, B.G., 2005. Simulation of the effects of groundwater level on vegetation change by combining FEFLOW software. Ecological Modelling 187, 341–351. Zhao, Z.Y., Wang, R.H., Zhang, H.Z., Shun, H.B., 2006. Degradation mechanism of desert ecosystem in lower reaches of Tarim River. Journal of Desert Research 26, 220–225.