Interactive comment on “Modeling hydrological processes influenced by soil, rock and vegetation in a small karst basin of southwest China” by Z.-C. Zhang et al

We would like to thank you and two reviewers for providing useful comments and suggestions on our submission to Hydrology and Earth System Sciences. We are pleased to report that we have now completed all the required revisions and are in position of submitting a new and improved version where we took into account all the comments provided. In this revision report we will provide details on how we dealt with each of the comments and suggestions. We will be very happy to address any other issues that you may judge necessary in addition to the corrections we already completed. Please


Introduction
Karst terrain covers about 15% of the world's land area or about 2.2 million km 2 , and is home for around 1 billion people (17% of the world's population) (Yuan and Chai, 1988).About 25% of the world's population is supplied largely or entirely by karst waters, including deep carbonate aquifers (Ford and Williams, 1989).One of the largest, continuous karst areas in the world is located in Yunnan-Guizhou Plateau of southwest China .It covers 540 000 km 2 and has a population of 100 million.In the karst region of southwest China, soils, developed on carbonate rocks, are generally 30-50 cm thick.Limestone fragments are mixed in soils and can act as a controlling factor for erosional rates and patterns in the landscape.
Since the beginning of China's economic reforms in 1978, China has experienced dramatic changes in its environment.Because its land is dominated by karst with thin soil, southwest karst region has experienced more significant environmental problems than other parts of the nation (Xiao et al., 2003).One of the most critical problems in karst areas is "rock desertification", which is the transformation of a karst area formerly covered by soil and vegetation into a rocky landscape or lithologic desert almost devoid of soil and vegetation (Yuan, 1997).In most cases, rock desertification is caused by human activities, especially clearing of natural vegetation.About 41.2% of the land experienced soil erosion and loss of water due to deforestation and agriculture.Karst systems typically have poor water storage capacity, which can be further diminished by deforestration (Huntoon, 1992).
From the hydrological aspect, in the vertical profile of a karst zone could be divided into three zones as shown in Fig. 1.An upper soil and epikarst zone (the near-surface weathered zone of the exposed carbonates), an infiltration zone and a deep conduit system.Because the uppermost layer soil is usually thin, epikarst zones in the karst basin take a key role for water storage and movement, as well as for transpiration of vegetation root systems.Because of effects of CO 2 supply and water flow on soluble carbonate rocks, the network of epikarst fissures through which percolation water passes is widened by dissolution near the surface.Therefore, near the surface epikarst has a large permeability, offering a fast water infiltration.As the extent and frequency of widening diminishes gradually with depth epikarst permeability diminishes with depth (Ford and Williams, 2007).Consequently, after recharge percolating rainwater is retained near the base of the epikarst leading to the formation of an epikarstic aquifer (Williams, 2008).The temporarily stored water flows via the vertical shaft system into 563 conduit systems of the deep aquifer.Groundwater in the deep aquifer flows through interconnected, solutionally enlarged conduits and cave systems and eventually flows out basin as a spring (Fig. 1).
The simulations of the karstic spring hydrographs were achieved by using two main modelling methods: (i) physically distributed models, and (ii) lumped models, with two different approaches either as black-or grey-box or reservoirs.Distributed models are usually based on numerical groundwater models, e.g.MODFLOW and FEFLOW, in the same way as for porous or fractured aquifers (Scanlon et al., 2003).They are improved to be able to represent the flow conditions in the saturated zone of karst aquifers in an underground conduit-and -fracture network through a porous matrix, like MODFLOW-DCM (Sun et al., 2005).However, application of these models in the karst area of southwest China is very difficult because most basins are located in the mountainous areas where observation data of groundwater table are not available.Moreover, the groundwater based models do not describe spatial distribution of soil moisture content, hillslope surface water or evapotranspiration from soil or epikarst.Therefore, these kinds of models are not suitable for estimation of eco-hydrological relationships.
Distributed hydrological models, like TOPMODEL (Beven and Kirkby, 1979) and DHSVM (Wigmosta et al., 1994), describe hydrological processes, especially surface, subsurface and groundwater flows.They are often used to evaluate impacts of climatic change and land use and land cover (LULC) on hydrology and water resources because they are easily developed to integrate with climate model and ecological processes.
In this study, we improved the distributed hydrology-soil-vegetation model (DHSVM) for describing hydrological processes in both porous and fissure mediums of the karst basin.In the improved DHSVM, a mixed runoff routing method integrating interactions among Darcy flow, fissure flow and channel storage routing was developed.The improved model was applied in a small watershed of Chenqi within the Puding Karst Ecohydrological Observation Station, Guizhou province of China (Fig. 2).The model was calibrated and validated based on the observation data of streamflow, vegetation interception and soil moisture contents.Modeling results reveal hydrological processes influenced by vegetation, soil and fracture features of the karst basin.

DHSVM
Distributed Hydrology-Soils-Vegetation Model (DHSVM) is a physically based distributed parameter model that provides a DEM based representation of a watershed (Wigmosta et al., 1994).The spatial scale of the representation is based on the grid size of the DEM.The model consists of a two-layer canopy (overstorey vegetation, understorey vegetation) representation of evapotranspiration, a two layer energy balance model for ground snow pack, a multilayer unsaturated soil model, a saturated subsurface flow model, and a three dimensional overland flow representation.The DEM provides the topographic controls on the incoming short-wave radiation, precipitation, air temperature, and downslope water movement.Grid cells are assigned vegetation and soil characteristics and are hydrologically linked through surface and subsurface flow routing.

Multilayer unsaturated zone model
The karst basin is divided, in the vertical direction, into three zones (Fig. 1): upper soil and epikarst, infiltration zone and deep aquifer.Since the infiltration flow via the vertical shaft system into the conduit system of the aquifer is fast, we neglected regulation of the infiltration zone.
Based on the original DHSVM structure, a multi-layer model is used for accounting soil moisture and epikarstic flow dynamics: where d 1 ,d 2 are soil or/and epikarst thickness of the upper and lower rooting zones, respectively, d 3 is the left epikarst thickness perched by saturated water, θ n (n=1, 2, 3) is the average soil moisture content of nth zones, P 0 is the volume of infiltrated rainfall, V sat is the volume of water supplied by a rising water table, V r is the volume of return flow (generated when a rising water table reaches the ground surface), E T o and E T u are evapotranspiration from overstorey vegetation (o) and understorey vegetation (u), respectively, E s is soil evaporation, P n (θ) (n=1,2,3) is downward volumes of water discharged from nth zones over the time step, Q t sin and Q t s is subsurface flow rate to pass in and out the epikarst zone, respectively.

Vertical infiltration and percolation
An equivalent hydraulic conductivity K v (θ) as follows could be used to compute infiltration and percolation P n (θ) based on Darcy's law assuming a unit hydraulic gradient: For the uppermost soil, the Brooks-Corey (Brooks and Corey, 1964) equation is used to calculate hydraulic conductivity: where m is the pore size distribution index, φ is the soil porosity, and θ r is the residual soil moisture content.For simplicity, the saturated moisture content θ s is taken equal to φ.
The epikarstic zone primarily consists of fractured medium.The mathematical formulation of the average velocity of groundwater flow in fractures is given by the "cubic law" (Kiraly, 1969;Snow, 1969).The transmissivity and hydraulic conductivity of fractural medium with smooth fractural surface can be expressed as follows: where b i ,j is the fracture aperture at grid i , j , µ is dynamic viscosity of water, ρ is fluid density, g is acceleration due to gravity.The equivalent hydraulic conductivity of a rock mass with one parallel set of fractures with rough surface is expressed by: where s is fracture spacing (m), C is a coefficient which is related with the roughness of the fracture surfaces and fracture aperture.
Infiltration and percolation P n (θ) is estimated by Eq. ( 4) using hydraulic conductivity K v (θ) estimated by Eq. ( 5) for the uppermost soil zone and Eq. ( 8) for the epikarst fracture zone.For the uppermost soil zone with partially rock exposure, an equivalent infiltration rate can be estimated by following equation: where f r is proportion of rock exposure.Hydraulic conductivity and transmissivity in the epikarst zone depend on fracture aperture and spacing.In an idealized fracture, with a 1 mm aperture and the same spacing, the hydraulic conductivity will be 10 −3 m/s, similar to that of loose clean sand (Singhal and Gupta, 1999).Limestone fracture aperture may be as large as 1×10 −2 m 567 in the exposed surface, decreasing to 1×10 −3 ∼1×10 −4 m in compact and impermeable bed rock.Epikarst porosity ranges between 1 (Smart and Friederich, 1986) and 2-10% (Gouisset, 1981) and 5-10% (Williams, 1985).

Saturated water flow in epikarst aquifer
In the saturated epikarst zone, subsurface water flowing cell by cell is calculated by the following equation with different values of horizontal transmissivities.
where w i ,j,k is the grid (i ,j ) width at k flow direction, T (t) i ,j is hydraulic transmissivity for grid (i ,j ), f i ,j,k is flow fraction for cell i ,j in direction k.
As the extent and frequency of fracture widening diminishes gradually with depth (Ford and Williams, 2007), the fracture aperture and spacing can be regarded as exponential decreasing towards the base of the epikarst, resulting in exponential decrease with depth of hydraulic conductivity (K s =K s i ,j e −αz ) Thus where K s i ,j is the lateral component of saturated hydraulic conductivity for cell i ,j at epikarst surface and is estimated by Eq. ( 8) with horizontal fractural aperture and spacing, z i ,j is the distance from the epikarst surface to the water table (positive downward), α i ,j is a parameter related to the decay of saturated conductivity with depth, and D i ,j is the total epikarst depth.For a soil matrix, lateral flow direction in the subsurface depends on hydraulic gradient which is usually assumed as topographic slope.However, in karst basins, strong anisotropy of karst fracture dominates subsurface flow directions as well.Hence, we developed the following formula to calculate flow fraction f i ,j,k for cell i ,j in direction k, considering topographic slope or groundwater slope β i ,j,k and aquifer anisotropy p i ,j,k : where p i ,j,k (= ) is a flow fraction for cell i ,j in direction k resulting from karst fractural anisotropy.Since the exact the location, geometry and hydraulic properties of each fracture are usually unknown, we applied statistical attributes of fractures, using a stochastic models, to describe spatial heterogeneity of rock fractures.In absence of field description of fractural anisotropy, the weight p i ,j,k can be randomly generated on the basis of uniform distribution of p-value between 0-1 in the eight directions.

Flow routing of underground conduit systems
The conduit systems within aquifers dominated by recharge through internal runoff often exhibit branch work patterns as small tributaries which, draining individual closed depressions, merge to form larger drainage channels which ultimately lead to a spring (White, 1999).Simulation of a karst system composed of dendritic paths (e.g., Milanovic, 1981;White, 1988;White and Deike, 1989) may require a great deal of sitespecific information for the karst channels and flow conditions (e.g., elevation, slope, fill material, roughness, cross sectional area, Reynolds number, Froude number, diameter, etc.) (Field and Nash, 1997;Field, 1997).
In this study, the dendritic paths are assumed to be similar to a surface channel network.The water may flow as open channel underground streams or as pipe flow with varying degrees of pressure head.Gradients within the matrix and fracture flow portions of the aquifer are orientated toward the conduits (Fig. 3).
Flow in the conduit systems is routed using a series of cascading linear channel reaches.Individual hydraulic parameters describe each reach.As the reach passes through grid cells, lateral inflow into the channel reach consists of the matrix and fracture flow.Flow is routed between channel reaches as a linear routing algorithm where 569 each reach is treated as a reservoir of constant width with outflow linearly related to storage.
where k is coefficient of underground channel storage, R r hydraulic radius, S o is hydraulic gradient, ∆L is channel segment, n is roughness, and ∆t is time step.

Study area and site description
The small catchment of Chenqi is located in the Puding basin of Guizhou with an area of 1.5 km 2 (Fig. 2).The study site has a subtropical wet monsoon climate with a mean annual temperature of 20.1 • C. The highest average monthly temperature is in July, and the lowest is in January.Annual precipitation is 1140 mm, with a distinct wet summerseason and a dry winter season.Average monthly humidity ranges from 74% to 78%.The elevation of the study area varies between 1320 and 1520 m above the sea level (Fig. 2).With the typical cone karst and cockpit karst geomorphology, the cone peaks of Chenqi catchment are generally 200 m above the adjacent doline depressions, and its catchment surface relief and slope are very steep.Geological properties include dolostone, thick and thin limestone, marlite and Quaternary soil (Fig. 2).
The vegetation species include grasses, karst montane deciduous broad-leaved shrubs and evergreens and deciduous broad-leaved mixed forest (Fig. 4a).According to field investigation by Pen et al. (2008), in the forest area, tree height varies from 2 to 5 m and its canopy fraction is about 80-90%; and rock outcrop rate is 20-40%.In the shrub area, tree height varies from 1 to 1.5 m and its canopy fraction is about 20-50%; and rock outcrop rate is 30-50%.In the grass area, the canopy fraction is about 80%; and rock outcrop rate is 35%.Rate of rainfall interception by forest is 13.2 by average, ranging between 3 and 25% for rainfall events with the total of 6.4-90.4mm.For shrub area, the rate of interception is about 9.7 on average.Distribution of soil properties, e.g.soil types, composition and density, are determined through field investigation and laboratory experiments from 49 soil samples.Generally, soil of the upper 20 cm is classified as two types of loam in the low and depression areas and sand loam in the middle and upper hill areas, with bulk density of 0.66 to 0.90 g/cm 3 and porosity of 0.2 to 0.4 (Fig. 4b).The lower layer is primarily the brown clay with most calcium content.Its density is 1-1.28 g/cm 3 and porosity 0.42-0.48.The average of soil thickness is 30 cm over the calcareous rock.Spatial distribution of soil thickness is strongly related to terrain, rock fragment and vegetative cover.Spatial distribution of the soil thickness is derived from the digital elevation model (DEM).The thickness ranges from 0 for exposed rock areas to 2.0 m for the lower cultivated areas.
In-situ experiments for the determination of saturated hydraulic conductivity (K v ) of soil were conducted at 15 measuring points using Guelph-Infiltration (GUELPH).The upper soil K v varies from 4.5×10 −5 m/s for loam to 5.67×10 −5 m/s for sand loam.
Shallow soil with thickness less than 15 cm overlying fractured bedrock significantly increases K v , which can be as large as 1.47×10 −4 m/s.Soil below 0.5 cm depth is usually compacted and K v is 1.0×10 −5 m/s (Chen et al., 2008).Field investigations show that epikarst thickness in Guizhou province, China,is about 2 m (Jang, 1998).Average thickness of the epikarst zone in the study catchment is estimated based on overland flow discharges collected by six troughs on four mountainous hillslops (Pen et al., 2008).Statistical analyses of the collected overland flow demonstrate that surface runoff only occurs when the cumulative amount of a successive rainfall events is larger than 80 mm during the wet season (Pen et al., 2008).According to the saturation mechanism, overflow only occurs when the whole epikarst 571 is saturated.Therefore, the amount of 80 mm can be regarded as a critical value for the whole epikarst zone saturation.The estimated average thickness of the epikarst zone is about 1.6 m with an effective porosity of 0.05.Spatial distribution of the epikarst thickness, ranging from 0.8 to 2.5 m shown in Fig. 4c is generated based on terrain for assumption that epikarst thickness is larger in flat areas and smaller on hillslopes (Tan, 2005).
Based on field observations of rock fracture distribution at 22 sites, the mean fracture aperture was estimated to be around 13.6×10 −3 m with standard variance of 15.7×10 −3 m and the mean fracture spacing was determined to be 0.142 m with standard variance of 0.0986 m.The exposed fractural apertures larger than 3×10 −3 m are mostly filled with soil.For apertures less than 3×10 −3 m, mean fracture aperture is 1.55×10 −3 m with standard variance of 0.705×10 −3 m.

Observation data and watershed digitalization
A set of meteorological and eco-hydrological observation stations were established since June 2007 within the Chenqi catchment.Three meteorological stations are located at hills of Dongjiashan (DJS), Houshaopo (HSP) and Yinpo (YinP) with forest, grass and shrub covers, respectively.Precipitation, air temperature and pressure, relative humidity, wind speed, and radiation are recorded in a time interval of five minutes.Four groups of the soil moisture probes at the depths of 15-30 cm below the ground surface were installed at DJS, HSP, YinPu and YanP, automatically recording soil moisture content every five minutes.50 rainfall collection tubes, with a diameter of 29 cm, were placed under the forest and shrub areas to collect penetrating rainfall.An automatic water level observation station was installed at the catchment outlet to record water level at every 15 min interval.
Meteorological observations demonstrate that forest canopy shades the soil and consequently prevents direct radiation absorption, resulting in lower temperature (T ) and higher humidity (Fig. 5).Daily irradiance at the ground surface within 10-20 cm high grass was 118.2 W/m 2 , much greater than 11.1 W/m 2 under 3-12 m high forest.The daily mean relative humidity was 86.3% and 95.5% for the grass and forest cover, respectively.The daily mean air temperature at 2 m height is 17.5 for the grass and forest cover, respectively.
The gridded DEM from model computation was derived from the 1:10 000 digital topography map from the Puding Karst Ecosystem Observation Station, Guizhou Province, China.The resulting rectangular subset 10 m×10 m resolution DEM containing 16 950 pixels (113 rows by 150 columns) was analyzed to delineate the basin boundary.In the DHSVM, a soil class is assigned to each model pixel.Weighted averages of percent sand, silt, and clay were calculated for each soil class.
Using statistical results of fractural features from field investigation, we mapped fractures onto a two-dimensional domain (x-y plane) with the third dimension assumed to be decreasing with depth to groundwater.Distributions of the fracture aperture and spacing on the ground surface in the vertical direction are randomly generated and follow a normal distribution function centered around 1.55×10 −3 m with a variance of 0.705×10 −3 m for aperture, and 0.142 m with a variance of 0.0986 m for fracture spacing.The weight p of fracture anisotropy is randomly generated on the basis of uniform distribution of p value between 0 and 1 in the eight directions.Thus, spatial distribution of K s value can be estimated by Eq. ( 8).
The underground channel drainage network was initially derived from the DEM on the basis of the assumption that large underground channels are located in the lower area of the basin.The channel network was further edited through field surveys.GIS analysis was used to determine the channel orders, and to assign mean values of channel width, depth, and Manning roughness for each order.Three orders of the underground channel network were generated, ranging 1∼3 m for channel depths and 0.3∼1.5 m for widths.573 5 Determination of model parameters

Procedures for parameter determination
Parameters of hydrological models are usually determined through model calibration according to basin outlet flow discharges.Calibration of physically based distributed models such as the DHSVM requires particular attention to the number of parameters used (e.g.Beven, 1989).Equivalence of different parameter combinations (parameter equifinality) is usually expected, and this causes non-uniqueness of the model calibration.The problem of parameter equifinality can, to some extent, be overcome by adding additional sources of data (e.g.parameterization results of other basins with similar land surface conditions, field experiment results) and by considering different hydrologic processes that can be separated using the observation data (e.g.soil moisture content and flow discharges).
In this study, the model parameters were determined as follows: (1) Parameters related with vegetation were determined on the basis of Land Data Assimilation System (LDAS) (http://ldas.gsfc.nasa.gov/LDAS8th/MAPPED.VEG/LDASmapveg.shtml) and the field observation of vegetation interception (Fig. 6).Except for hydraulic conductivity, parameters related with soil properties, such as porosity, field capacity and wilting point, were determined by field experiments and laboratory analysis.(2) Field measured soil hydraulic conductivities are used as initial values for model calibration.The final values of soil hydraulic conductivities and other soil related parameters (e.g.α and m in Eq. 6) are further calibrated making use of SMC data measured at the four sites (Fig. 7).( 3) Hydraulic conductivities of epikarst were initially estimated according to observed flow discharges using Darcy's Law.These hydraulic conductivities and conduit roughness were further calibrated using the basin outflow hydrograph (Fig. 9).The calibrated hydraulic conductivities represent catchment average values.Their spatial distribution was described according to spatial statistical results of fractural features.
In each calibration step, the four large flood discharges from 28 July 2007 to 19 October 2007 were used for the model initiation and calibration, and the period from 20 April 2008 to 26 June 2008 was held back for the verification.The computation step is 1 h.

Estimation of hydraulic conductivities of epikarst
Hydraulic conductivity of the fracture medium can be roughly estimated according to watershed flow discharge (Moore, 1997) if the flow discharge is primarily from fractures.The discharge is adequately described by Darcy's Law.The general expression is where T t is hydraulic transmissivity, I t is the average slope of the study basin, A is drainage area, and a is path length, the average distance from a drainage divide to a stream (Moore, 1992).
From field observations and map measurements, the average distance from a drainage divide to a flowing stream in the study watershed is estimated to be about 1225 m.The average hydraulic gradient is about 0.268 for most mountainous areas with average slope angle of 15 0 .If the entire basin is assumed to be contributing, Eq. ( 16) can be used to calculate horizontal transmissivities during the low flow and the high flow periods using the minimum flow discharge and the maximum flood discharge of 0.01 and 2.28 m 3 /s, respectively.
The computed transmissivities are 3.05×10 −5 for the lowest flow and 6.95×10 −3 m 2 /s for the highest flow.The horizontal K s is estimated between 1.90×10 −5 and 4.34×10 −3 m/s for an epikarst thickness of 1.6 m.Fracture apertures with smooth fracture surface (C=1) can be thus computed by Eq. ( 8).A K s value of 1.90×10 −5 is equivalent to the lowest flow via small fracture apertures of 0.286 mm, and 4.34×10 −3 m/s is the highest flow via the large fracture apertures of 1.75 mm.
The above results indicate that drought flow from the epikarst aquifer with a lower groundwater table is primarily controlled by a low permeability matrix while the storm 575 response in a high groundwater table is significantly controlled by a high permeability fracture system.Thus vertical decrease of K s can be expressed according to vertical variation of fracture apertures.If fracture aperture exponentially decreases from 1.75×10 −3 m on the epikarst surface to 0.286×10 −3 m at the base rock of 1.6 m depth, vertical variation of fracture aperture is assumed as b = b 0 e −1.13z , where b 0 is fracture aperture on the epikarst surface.The final spatial mean value of K s is determined by model calibration for adjusting the C in Eq. (8).

Simulation of soil moisture contents
The dynamics of soil moisture content (SMC) are influenced by soil properties, upper and lower boundary conditions: land surface including climate and vegetation, and underlying epikarst.Climate and vegetation canopy influences can be calculated using the two-layer canopy evapotranspiration package of DHSVM.The underlying epikarst boundary can be assumed as a large infiltration rate without water ponded on the soilepikarst interface because epikarst fracture K v is usually much larger than that of the upper soil, and thus the underlying permeability does not influence downward water percolation from the upper soil zone (P 1 (θ) in Eq. 1).Therefore, we focus on calibration of the upper soil hydraulic conductivities which strongly influence infiltration and water flow, and thus soil moisture content.If the horizontal hydraulic conductivity K s is assumed to be same as the vertical hydraulic conductivity K v , the hydraulic conductivities of surface soil K v in Eq. ( 5) is calibrated as 1.5×10 −4 and 1.6×10 −4 m/s for loam and sand loam, respectively.Exponential decrease coefficient α is 3 and the pore size distribution index m for Eq. ( 5) is 0.19 and 0.20 for loam and sand loam, respectively.Figure 7 shows that simulated and observed SMCs at DJS, YinP, YanP and HSP have similar change patterns, represented by the normalized SMC.However, the magnitude of simulated SMC is different from of the observed SMC because simulated SMC represents the average value within the grid of 10 m×10 m, while observed SMC represents a point value from the sensor.Variations of the normalized SMC in the four sites reflect different behaviors of soil dynamics (Eqs.5 and 6) response to different surface shadings in a same climate forcing.Soil moisture variations indicate that vegetation cover inhibits soil moisture loss of the shallow soil layer through evaporation because the canopy shades the soil and consequently prevents direct radiation absorption (Fig. 5).SMC variations for DJS with forest are relatively smooth than those on other hillslopes.Simulated soil evaporation for the forest area in DJS and the grass area in HSP during 20 August to 19 October 2007 (Fig. 8) shows that soil evaporation in the forest area is much less than in the grass area, resulting in less variations of soil moisture content during the drought periods.

Simulation of streamflow discharges
Generally, the influence of small storage of the thin soil on the hydrograph is not significant compared with epikarst and conduit system in the study area.Epikarst hydraulic conductivity and conduit roughness are regarded as two key parameters controlling the hydrograph.They can be further calibrated based on measured underground flow discharge.For underground flow discharge simulation, three objective functions are selected for evaluating the success of the calibration: the relative volume error R ve between the observed and simulated flows, the root mean squared error RMSE, relating how well the calculated and observed hydrographs compare in both volume and shape, and the coefficient of determination D, relating how well the calculated hydrograph compares in shape to that observed and depends only on timing, not on volume.
The conductivity profile suffice to match the outflow hydrograph with reasonable accuracy (Fig. 9).If the relationship b = b 0 e −1.13z is adopted for estimation of horizontal K s variation along the vertical direction and an anisotropy factor K v /K s is assumed as 1.The calibrated roughness coefficient of the fractural surfaces C is 0.5 and thus estimated mean value of horizontal hydraulic conductivities at epikarst surface K s is 2.22×10 −3 m/s.Roughness which controls hydrograph oscillation is 0.025.

Effects of epikarst fractural width on hydrological processes
The hydraulic conductivity and discharge are a function of, and directly proportional to, the size of the aperture (i.e. an increase the aperture results in an increase in the hydraulic conductivity and an increase in the discharge).For scenarios of the fracture apertures with mean values of 0.25 mm, 0.5 mm, 1 mm and 1.5 mm and standard variance of 1/2 mean values, simulated flow discharges shown as Fig. 10 demonstrate that 0.42 m 3 /s of peak flow discharge for the aperture width of 0.25 mm increases to 2.29 m 3 /s for the width of 1.5 mm.The peak discharge for the aperture larger than 1.0 mm occurs 2 h earlier than those of the smaller apertures.Figure 10 shows that flow recession becomes slower as the widths decrease.
Model simulation results show that fractural aperture controls the flood volume through influencing infiltrated rainfall division between unsaturated and saturated zones.A high density of fractures results in large infiltration and sharp hydrograph response and thus little infiltrated water remains in soil and gets eventually lost through evapotranspiration.An example shown in Fig. 10 is the simulated streamflow for a scenario where the epikarst zone is replaced with upper soil layer, compared with observed and simulated streamflow with epikarst during 29-31 July 2007.Basin average SMC of the upper and lower layers of the root zone overlying epikarst is 0.135 and 0.368, respectively, less than 0.158 and 0.390 when epikarst zone is replaced with soil.Consequently, basin actual evapotranspiration increases from 0.047 to 0.070 mm/h.

Conclusions
Because soil is very thin in the karst region of southwest China, hydrology in the karst basin is closely related to basin characteristics of fissured media in addition to topography, soil and vegetation.Modeling of water flow in the karst basin should be based on runoff routing integrating Darcy's flow, fissure flow and channel storage routing.In this study, a distributed hydrological model focused on epikarst fracture hydrological dynamics is developed on the basis of the distributed hydrology-soil-vegetation model (DHSVM).The improved model was applied in a small watershed of Chenqi within the Puding Karst Ecohydrological Observation Station.
As the developed model becomes complex, it increases additional parameters related with fractural characteristics.Since we can not measure the fractural characteristics at every site, a statistical method based on field investigation at typical sites is a better way to generate spatial distribution of the fracture features.The fractural parameters influenced flow discharges were estimated from the flow discharge data and was further calibrated based on the observation data of streamflow.
The results show that this improved model successfully captures the sharp variation in the hydrograph.It can also simulate soil moisture content and evapotranspiration associated with precipitation, karst fracture and vegetation cover.This study has implication for the analyses of ecohydrological processes and is likely to help understand the response of these unique systems to major changes in land cover and land use.

Fig. 9 Fig. 9 Fig. 9 .
Fig.9 Observed and simulated streamflow for model calibration (a) and validation (b) Fig.10 Comparison of simulated streamflow for different epikarst apertures
model computation results demonstrate: R ve =−0.018,RMSE=0.05m 3 /s, and D=0.847 for the calibration period from 28 July 2007 to 19 October 2007; R ve =0.077, RMSE=0.01 m 3 /s, and D=0.836 for the validation period from 20 April 2008 to 26 June 2008.These results illustrate that adjustments to the lateral hydraulic 577 Observed and simulated results of normalized SMC.
8 Fig. 7. Fig.8 Simulated results of normalized SMC and soil evaporation (Es) in the forest and grass areas Fig. 8. Simulated results of normalized SMC and soil evaporation (E s ) in the forest and grass areas.