water

Article
Watershed Characterization and Hydrograph Recession
Analysis: A Comparative Look at a Karst vs. Non-Karst
Watershed and Implications for Groundwater Resources
in Gaolan River Basin, Southern China
Hamza Jakada 1 , Zhihua Chen 1, *, Mingming Luo 1 , Hong Zhou 2 , Zejun Wang 2 and
Mukhtar Habib 3
1 School of Environmental Studies, China University of Geosciences, Wuhan 430074, China;
jakadahamza@gmail.com (H.J.); luomingming@cug.edu.cn (M.L.)
2 Geological Survey of China University of Geosciences, Wuhan 430074, China; zhouhong@cug.edu.cn (H.Z.);
wangzejun@cug.edu.cn (Z.W.)
3 Department of Minerals and Petroleum Engineering Kaduna Polytechnic, P.M.B 2021 Kaduna State, Nigeria;
mukhtarhabib102@gmail.com
* Correspondence: zhchen@cug.edu.cn; Tel.: +86-139-0716-6202




Received: 24 February 2019; Accepted: 2 April 2019; Published: 10 April 2019


Abstract: Karst watersheds are often treated as non-karst watersheds that can lead to several
hazards. Hence, how do karst watersheds differ from non-karst watersheds and what are the
effects of karstification on groundwater availability and quality? In this study, we contrast between
a karst and non-karst watershed by elucidating their geomorphological peculiarities and potential
impact on spatio-temporal availability and quality of groundwater. GIS morphometric mapping and
hydrograph recession analysis are applied to map the watershed features and estimate hydrograph
recession coefficient to define the groundwater drainage characteristics as well as the influence of
karst drainage attributes (KDA). Furthermore, we characterize streamflow components based on
the hydrograph recession limbs (segments) and infer their contributing geomorphological factors.
Results show that the karst watershed has higher recession coefficients for successive recession
limbs. Consequently, it drains larger volumes of groundwater primarily due to the KDAs, which
transmit interflow and groundwater flow more rapidly through large cavities to springs as well
as stream channels. The KDAs generate what we term karst drainage flow (KDF), defined by the
second recession limb which has high recession coefficient as the first limb (overland flow) and
strongly contrasts with the non-karst watershed from visual and ANOVA analysis. The effect is that
karst aquifer yield over time is significantly lower and highly exposed to pollution compared to
the non-karst aquifer. Consequently, sustainable water management practices should be adopted to
ensure the availability and safety of groundwater reserves.

Keywords: karst; watershed; hydrograph recession analysis; GIS; groundwater

1. Introduction
Groundwater resources in karst regions are under increased pressure from rapid population
growth and increased agricultural as well as industrial activities [1]. In China alone, over 100 million
people rely on karst water resources [2]. The multiple porosity matrix in carbonate rocks (limestone,
dolostone, gypsum, quartzite, and halite) gives karst aquifers the capacity to house large volumes
of water, but also, this extreme heterogeneity (secondary and tertiary porosities) from fractures and
underground conduits can result in enhanced flow thereby reducing storage and also increase aquifer

Water 2019, 11, 743; doi:10.3390/w11040743 www.mdpi.com/journal/water

Water 2019, 11, 743 2 of 20

sensitivity to pollution [3]. Therefore, due to these special characteristics of karst systems, accurate
assessment of quality and quantity of karst water is imperative. Unfortunately, karst groundwater
systems are often treated as that of non-karst areas, which has led to groundwater pollution risks as
well as other environmental and geological hazards [4,5].
Generally, watersheds are bound by high relief which serves as drainage for all precipitation
that occurs within it to a single outlet. Therein, precipitation is the most potent catalyst that drives
geomorphological work. Typical watersheds comprise of some common components such as a natural
stream network, defined catchment divide, general orientation downstream as well as composite
geology. These components and their morphometric characteristics define the general catchment
response to rainfall events. The primary groundwater drainage controls are topography and geology,
with other essential factors being basin physiographic characteristics and stream network length and
density. In this study, we focus mainly on the geomorphological peculiarities since karst geological
features have direct impacts on groundwater infiltration and drainage, which can lead to large yield
due to higher hydraulic conductivity. Conversely, non-karst terrains such as crystalline formations,
are known for absences of primary porosity and lower hydraulic potential, but lend itself to the
development of thick saprolite overburden that can hold and convey large quantities of water over
time [6]. Thus, in basins where the underlying geology is crystalline (hard rock) with exceptionally
high regolith and fluvial deposit, this sediment layer can retain large volumes of delayed flow and as
such can become a more critical baseflow-sustaining reservoir than the underlying solid bedrock [7].
In this study, we compare between two watersheds, a karst (Miaogou) basin characterized by
sedimentary carbonate rocks and a non-karst (Gaojiaping) basin of mainly crystalline rocks. The karst
watershed is our primary focus as the objective is to determine how karst geomorphological features
impact the interaction between surface and groundwater hydrology in contrast to the non-karst
watersheds. Furthermore, we assess the consequences on groundwater availability and quality, as
these issues have not been adequately addressed in previous literature. Using streamflow data collected
at the outlets of both Miaogou and Gaojiaping watersheds, we analyzed the storm hydrographs of
both watersheds with particular focus on the recession curve to determine the recession coefficient of
successive limb segments on the curve. Streamflow recession curve holds important information
concerning the watershed geology, storage and aquifer properties [8–10]. Geomorphological
characteristics influencing the hydrograph can be inferred from hydrograph recession curve analysis
due to its correlation with watershed geomorphology and hydrogeology [11–17]. For example,
Kullman [15], Luo et al. [16], and Fiorillo et al. [17] applied hydrograph recession analysis to
characterize properties of karstic aquifers. Also, Malı’k [18] applied the recession coefficient to
determine regional karstification degree as well as groundwater sensitivity to pollution. Hence,
recession curve analysis can improve our understanding of low flow periods and aquifer contribution
to streamflow which is significant for ensuring sustainable water resources management. Additionally,
it is noteworthy that, recently, residents of Miaogou (karst) watershed have been experiencing a
shortage of groundwater especially during dry periods of the hydrologic year. This has raised concerns
on the availability and quality of groundwater resources in the area, which makes this particular study
of importance.

Study Area
Gaolan River basin is located in Xingshan County, western Hubei Province, Southern China. It is
characterized by a subtropical monsoon climate with abundant rainfall that averages 900–1200 mm
annually. These climate parameters have a direct influence on the karst geomorphology since
karstification process is water dependent. The two case study watersheds, Miaogou and Gaojiaping,
are less than 9 km apart from each other in Gaolan River basin. Miaogou karst watershed is underlain
by a well-developed karst terrain with a surface elevation that ranges between 430 to 1770 m.a.s.l.
The predominant carbonate rock in the area is limestone of Cambrian age, which appears like
rounded hills or isolated towers (typical of tropical karst). Non-carbonate rocks present in the area are
Water 2019, 11, 743 3 of 20

Water 2019, 11,

mudstone andx FOR PEER
shale ofREVIEW
Silurian and Lower Cambrian ages, respectively (Figure 1). The karst areas 3 ofare
20
Water 2019, 11, x FOR PEER REVIEW 3 of 20
characterized by large depressions, sinkholes, caves, and several springs, some of which are shown
are characterized
in Figure by large depressions,
2. The landcover is predominantlysinkholes, caves,ofand
a mixture several springs,
Agriculture land and some of which
Evergreen are
Forest.
are characterized
shown in Figure 2. by
Thelarge depressions,
landcover is sinkholes, caves,
predominantly a and several
mixture of springs, some
Agriculture land ofandwhich are
Evergreen
Gaojiaping (non-karst)
shown in Figure 2. Thewatershed
landcover has about 75% of aits
is predominantly entire of
mixture catchment
Agriculture underlain
land andbyEvergreen
an igneous
Forest.
rock Gaojiaping
while 15% and (non-karst)
10% is watershed and
metamorphic has about 75%rocks,
carbonate of its respectively.
entire catchment The underlainigneous
exposed by an
Forest. Gaojiaping (non-karst) watershed has about 75% of its entire catchment underlain by an
igneous
and rock whilestrata
metamorphic 15% andmainly
10% is metamorphic and carbonate rocks, respectively. The exposed
igneous rock while 15%are and 10% isProterozoic
metamorphic where
and the igneous
carbonate comprises
rocks, of granite
respectively. The and diorite.
exposed
igneous
The and
metamorphicmetamorphic
rocks arestrata are
comprised mainly
mainlyProterozoic
of slant where
gneiss, the igneous
graphite comprises
gneiss, gneiss, of
andgranite and
granulite.
igneous and metamorphic strata are mainly Proterozoic where the igneous comprises of granite and
diorite.
This The metamorphic
composition collectively rocks arethe
forms comprised mainly of slant gneiss,
ancient crystalline graphite gneiss, landcover
gneiss, and
diorite. The metamorphic rocks are comprised mainly of basement
slant gneiss,in the area where
graphite gneiss,the
gneiss, and is
granulite.
mainly This
evergreen
granulite.
composition
forest.
This composition
collectively forms the ancient crystalline basement in the area
collectively forms the ancient crystalline basement in the area where the
where the
landcover is mainly evergreen
landcover is mainly evergreen forest.forest.

Figure 1. Major geological units in Gaolan River basin, Hubei Province, South China (Source: China
Figure 1.
Figure Major geological
1. Major geological units
units in
in Gaolan
Gaolan River
River basin,
basin, Hubei
Hubei Province,
Province, South
South China
China (Source:
(Source: China
China
Geological Survey).
Geological Survey).
Geological Survey).

Figure2.2.Some
Figure Some karst featureslocated
karst features locatedinin Miaogou
Miaogou watershed
watershed and and environs
environs (A) Xiaosanyou
(A) Xiaosanyou Cave;
Cave; (B)
(B)doline
dolinewith
withdotted
dottedlines
lines enunciating the perimeter; (C) large sinkhole with dotted lines enunciating
enunciating the perimeter; (C) large sinkhole with dotted lines enunciating
the
theperimeter;
Figure 2. Some (D)
perimeter; (D) vracture
karst with
featureswith
vracture dotted
located lines
in lines
dotted Miaogou enunciating
watershed
enunciating the
andorientation
environs
the orientation of Xiaosanyou
(A) fracture;
of fracture; (E)Cave;
epikarst;
(E) epikarst; and(B)
and (F)
doline bailongquan spring.
with dottedspring.
(F) bailongquan lines enunciating the perimeter; (C) large sinkhole with dotted lines enunciating
the perimeter; (D) vracture with dotted lines enunciating the orientation of fracture; (E) epikarst; and
2. Data and Methods
(F) bailongquan spring.

2. Data and Methods

Water 2019, 11, 743 4 of 20

2. Data and Methods

2.1. GIS Watershed Characterization

Water 2019, 11, x FOR PEER REVIEW 4 of 20
Several studies have applied GIS and Remote Sensing data and methods for different terrains, and
they
2.1. have proved Characterization
GIS Watershed to be efficient for the characterization of drainage parameters [19–21]. The digital
elevation
Several studies haveand
model (DEM) its derivative
applied datasetSensing
GIS and Remote (slope, data
aspect
andetc.) constitute
methods the building
for different blocks for
terrains,
defining
and theythe watershed
have proved to boundary, flow
be efficient for direction, and stream
the characterization network,parameters
of drainage which are [19–21].
generally
Theprepared
asdigital
inputelevation
for hydrologic models
model (DEM) and[22–29]. Extrapolating
its derivative the stream
dataset (slope, order
aspect etc.) is conducted
constitute based on the
the building
hierarchical order method
blocks for defining as proposed
the watershed boundary, byflow
Strahler [30]and
direction, from flownetwork,
stream accumulation
which areand flow direction
generally
preparedderived
data-sets as inputbased
for hydrologic models [22–29].pour
on the eight-direction Extrapolating the stream
point algorithm order is conducted
[29,31–33] as a modulebased
in ArcGIS
on the analysis
surface hierarchical order DEM
utilizes method as proposed
raster by Strahlerflow
for determining [30] direction
from flowof accumulation
grid cells outandofflow
their eight
direction data-sets
neighbouring cells.derived based on the
Extrapolating the eight-direction
stream order pour point algorithm
is conducted based[29,31–33] as a moduleorder as
on the hierarchical
in ArcGIS surface analysis utilizes DEM raster for determining flow direction of grid cells out of their
proposed by Strahler [30]. This is achieved from the flow accumulation and flow direction data-set
eight neighbouring cells. Extrapolating the stream order is conducted based on the hierarchical order
as derived from the DEM. Since several studies have been carried out describing these methods we
as proposed by Strahler [30]. This is achieved from the flow accumulation and flow direction data-
will
set not dwell from
as derived on itthe
here.
DEM.For a more
Since detailed
several studies description,
have been carried the reader may refer
out describing thesetomethods
Reference [29].
Since several studies have been carried out describing these methods, we will not
we will not dwell on it here. For a more detailed description, the reader may refer to Reference [29]. dwell on it here.
For a more
Since detailed
several studiesdescription, the reader
have been carried may refer
out describing to methods,
these Dixon and weUddameri [29].
will not dwell onFigure
it here. 3 shows
m a×more
30For 30 mdetailed
resolution SRTM DEM
description, datamay
the reader for both study
refer to Dixon areas.
and Uddameri [29]. Figure 3 shows 30
m × 30 m resolution SRTM DEM data for both study areas.

Figure3.3. SRTM
Figure SRTM DEM
DEM data (30 × 30for
data (30x30m) m)Miaogou
for Miaogou (A)Gaojiaping
(A) and and Gaojiaping (B) watersheds
(B) watersheds (source: (source:
https://earthexplorer.usgs.gov/).
https://earthexplorer.usgs.gov/).

2.2.
2.2.Hydrograph
Hydrograph Recession CurveAnalysis
Recession Curve Analysis
AsAsstated
stated earlier,
earlier, the
therecession
recession curve
curverefers to the
refers partpart
to the of the
ofhydrograph that comes
the hydrograph that after
comes theafter the
crest(peak).
crest (peak).The
The recession
recessioncurve
curvecan bebe
can characterized
characterizedby aby
uniform recession
a uniform limb orlimb
recession a setor
of turbulent
a set of turbulent
segmentsas
segments asthe
the curve
curve regresses
regressesover
over time. TheThe
time. segments or inflection
segments pointspoints
or inflection are generally associated
are generally associated
with the non-linearity in rock and soil media storage structure which that contributes to streamflow,
with the non-linearity in rock and soil media storage structure which that contributes to streamflow,
generally in the form of interflow and baseflow. The segments on the recession curve have been
generally in the form of interflow and baseflow. The segments on the recession curve have been
analyzed in the past using several methods, including the exponential method to determine the
analyzed
dischargein the pastthat
processes using
makeseveral methods,
up the total including
streamflow. the exponential
The exponential functionmethod to determine the
was first expressed
discharge processes that make up the total
by Boussinesq [34] and Maillet [35], expressed as; streamflow. The exponential function was first expressed
by Boussinesq [34] and Maillet [35], expressed as;
𝑄 =𝑄 𝑒 (1)
0 −αt
Qt = Q e recession and base-flow analyses for low
Since then, it has been widely used for hydrograph (1)
flow forecasting, drainage basin studies, pollution sensitivity studies [12,14–18,36,37]. This model can
Since athen,
describe it has closed,
relatively been widely used for
independent andhydrograph recession
large reservoirs whoseand base-flow
main analyses
water source is for low
flow forecasting,
precipitation anddrainage basin
major losses arestudies,
throughpollution
streamflow sensitivity
dischargestudies [12,14–18,36,37].
[16]. Miao This model can
and Miao [38] showed
how it could
describe be used
a relatively to describe
closed, the drainage
independent of thereservoirs
and large aquifer reserves over time.
whose main waterThe exponential
source is precipitation
equation can be expanded and further expressed for multiple time-steps as:
and major losses are through streamflow discharge [16]. Miao and Miao [38] showed how it could
Water 2019, 11, 743 5 of 20

be used to describe the drainage of the aquifer reserves over time. The exponential equation can be
expanded and further expressed for multiple time-steps as:
Water 2019, 11, x FOR PEER REVIEW 5 of 20
logQ0 − logQt
α= (2)
0.4343t
𝑙𝑜𝑔𝑄 − 𝑙𝑜𝑔𝑄
where Q0 is the initial flow; Qt is the flow 𝛼at = any moment,
0.4343𝑡
(2)
and α is the recession coefficient (coefficient
− 1
of discharge) with dimension [T ]. It is estimated by plotting the recession limb in a semi-logarithmic
where
plot where Q0 is the
Equation (1)initial
appears flow;
as Q is the flow at
a tstraight-line anywith
plot moment, andslope
constant α is the recession coefficient
independent of the initial
Q0(coefficient
value [16].ofBased
discharge)
on thewith dimension scale,
logarithmic [T ]. Itthe
−1 is estimated
recessionby plotting
curves canthe recessionbe
generally limb in a semi-
classified into
logarithmic plot where equation 1 appears as a straight-line plot with constant slope
several segments, as shown in Figure 4, based on inflection points which are indicative of a transition independent of
the initial Q0 value
(micro-regimes) [16].drainage
in the Based onstructure.
the logarithmic scale, the recession
Consequently, we cancurves can generally
infer major sourcesbeand
classified
storage
into several segments, as shown in Figure 4, based on inflection points which
systems such as overland and subsurface flows as well as estimate specific volumes of each segment are indicative of a
transition (micro-regimes) in the drainage structure. Consequently, we can infer major sources and
(area under the curve) (Figure 4). For heterogeneous aquifers, on a small scale, especially in areas with
storage systems such as overland and subsurface flows as well as estimate specific volumes of each
highly karstified media, the recession curve can be classified into three (3) or more segments [16,17].
segment (area under the curve) (Figure 4). For heterogeneous aquifers, on a small scale, especially in
Aquifers with a linear storage function have been reported to have no significant inflections and thus
areas with highly karstified media, the recession curve can be classified into three (3) or more
can be hardly segmented. Therefore, it can be understood that the non-linear recession curve cannot
segments [16,17]. Aquifers with a linear storage function have been reported to have no significant
be inflections
solved onlyand by thus
Equation
can be(1).hardly
However, by assuming
segmented. each itrecession
Therefore, segment as that
can be understood one the
being “constant”
non-linear
and the second
recession “decaying,”
curve then itonly
cannot be solved canby beEquation
solved using the formby
(1). However, expressed
assumingineach
Equations
recession(3)segment
and (4).
as one being “constant” and the second “decaying,” then it can be solved using the form expressed
logQi−1 − logQi
in Equations (3) and (4). αi = (3)
0.4343(ti − ti−1 )
 𝑙𝑜𝑔𝑄 − 𝑙𝑜𝑔𝑄
𝛼=Q e−α1 t (0 ≤ t ≤ t ) (3)
 0.4343(𝑡
 1 −𝑡 ) 1

Qt = Q 2 e − α2 t ( t 1 ≤ t ≤ t 2 ) (4)
 𝑄 𝑒−α t (0 ≤ 𝑡 ≤ 𝑡 )
 Q e ((𝑡t2 ≤ t≤t )
 3
𝑄 = 𝑄3 𝑒 ≤ 𝑡 ≤ 𝑡 3)

(4)
𝑄 𝑒 (𝑡 ≤ 𝑡 ≤ 𝑡 )

Figure
Figure 4. 4. Schematicpresentation
Schematic presentationof
ofthe
thehydrograph
hydrograph recession
recessionlimb
limbwith
withfour
fourrecession
recession segments and
segments and
corresponding
corresponding volumesofofdischarged
volumes dischargedwater
water(adapted
(adapted from
from Kresic
Kresic N.
N. and
andBonacci
BonacciO.O.[39]).
[39]).

From
From equation(4),
Equation (4),considering
consideringQ
Q is known, α
is known, α (recession
(recessioncoefficient)
coefficient)can
canthen
thenbebe
solved forfor
solved each
each
segment. Additionally, each segment can be quantified and uniquely defined by integrating
segment. Additionally, each segment can be quantified and uniquely defined by integrating the flow the flow
at initial time step t0, which represents the lower limit and the next time step being the upper limit,
while maintaining the recession coefficients for each time step, respectively, expressed as;
Water 2019, 11, x FOR PEER REVIEW 6 of 20
Water 2019, 11, 743 6 of 20

⎧ 𝑉 = the(𝑄
at initial time step t0 , which represents 𝑒 limit
lower − 𝑄 and𝑒 the )𝑑𝑡next time step being the upper limit,
⎪
while maintaining the recession ⎪ coefficients for each time step, respectively, expressed as;
𝑉 = (𝑄 𝑒 −𝑄 𝑒 )𝑑𝑡 (5)
⎨ R t1 − t − t
⎪ V1 = 0 ( Q1 e
α 1 − Q2 e α 2 )dt
⎪

𝑉 = R(𝑄
t2 𝑒 −𝑑𝑡 ) −
⎩ V2 = 0 ( Q2 e 2 − Q3 e 3 )dt
α t α t (5)

 R∞
V3 = ( Q3 e−α3 t dt)

In addition, the total percentage volume0 of each component (Vi, i=1, 2, 3) to total water storage
capacity (V0) is calculated
In addition, the totalas;
percentage volume of each component (V , i=1, 2, 3) to total water storage
i
capacity (V 0 ) is calculated as;
3
𝑉 V= 𝑉 Vi
0 = ∑ (6) (6)
i =1

V
Ki =𝑉 i × 100% (7)
𝐾 = ×0 100%
V (7)
𝑉
2.3. Streamflow Data
2.3. Streamflow Data data were collected using an automatic ultrasonic water level gauge (CJ 800)
The streamflow
stationed at the outlet
The streamflow dataofwere
bothcollected
Miaogou andanGaojiaping
using watersheds.
automatic ultrasonic waterBylevel
visually
gaugeanalysing
(CJ 800) the
hydrograph
stationed at the trend (Figure
outlet 5), especially
of both Miaogou and the Gaojiaping
recession limb, it would
watersheds. Bybe noted analysing
visually that for thethekarst
watershed,
hydrographthe recession
trend (Figureis 5),
smoother andthe
especially lessrecession
turbulentlimb,
than itthat of thebenon-karst
would noted thatwatershed. This can
for the karst
be attributed
watershed, thetorecession
the very islarge subsurface
smoother karst
and less conduitsthan
turbulent thatthat
easily drain
of the percolated
non-karst precipitation
watershed. This out
of
canthe
bewatershed.
attributed toThe
therecession
very largelimb of the non-karst
subsurface watershed
karst conduits is more
that easily drainturbulent and
percolated thus, indicative
precipitation
outa of
of morethe segmented
watershed. drainage
The recession limb
due to someof the non-karst
active watershed is impediments.
geomorphological more turbulentInand thus,
other words,
indicative of a more segmented drainage due to some active geomorphological impediments.
the karst watershed seems to have a more linear reservoir system. Also, the data shows that the karst In other
words, the generally
watershed karst watershed
has lessseems to have
volume a more linear
of discharge duringreservoir system.
dry periods Also, the data shows
(November-March, whichthat
marks
the dry
the karstwinter
watershed
seasongenerally has less volume
in the catchment of discharge
area). During during groundwater
this period, dry periods (November-March,
levels are generally at
whichlowest
their marksmaking
the dryabstraction
winter season in difficult,
more the catchment area).inDuring
especially Maiogou.this Thus,
period,fifteen
groundwater levelswere
hydrographs
are generally at their lowest making abstraction more difficult,
selected from 2015 and 2016 for the recession curve analysis for this study. especially in Maiogou. Thus, fifteen
hydrographs were selected from 2015 and 2016 for the recession curve analysis for this study.

Streamflowdata
Figure5.5.Streamflow
Figure dataforfor Miaogou
Miaogou and
and Gaojiaping
Gaojiaping Watersheds
Watersheds for 2015
for 2015 (A) and
(A) and 20162016
(B). (B).

3. Results and Discussion

Water 2019, 11, 743 7 of 20

3. Results and Discussion

Water 2019, 11, x FOR PEER REVIEW 7 of 20

3.1. Topographic Analysis

3.1. Topographic Analysis
Topography is defined by the general relief of an area, slope as well as the drainage divide which
Topography is defined by the general relief of an area, slope as well as the drainage divide which
delineates the surface water basin. The hypsographic curve of Miaogou watershed shown in Figure 6
delineates the surface water basin. The hypsographic curve of Miaogou watershed shown in Figure
illustrates
6 illustratescomplex
a a complex surface
surfacerelief
reliefwhere
wherethe theminimum
minimum altitude
altitude is 436m while
is 436m whilethe themaximum
maximumisis1780 1780 m
inman 2
in area of only
an area of only4545
km km. 2The
. Thefrequency
frequency polygons
polygons arearedivided
dividedintointo
1010 topographic
topographic classes
classes for both
for both
watersheds.
watersheds. In the karst catchment, there are considerable amounts of low to flat-lying areas, which are
In the karst catchment, there are considerable amounts of low to flat-lying areas, which
actually karstkarst
are actually depressions
depressionsthat that
extend between
extend betweenkarst cones
karst (domes).
cones (domes).About
About 8.5%8.5%of of
thethe
total area
total of the
area
watershed ranges ranges
of the watershed between 436–770
between m while
436–770 m 63%
whilerange between
63% range 840–1100
between m andm31.5%
840–1100 rangerange
and 31.5% between
1175–1780 m. This will
between 1175–1780 havewill
m. This a substantial impact on
have a substantial time on
impact of concentration for runoff
time of concentration for as wellas
runoff aswell
erosion
and sediment
as erosion andtransport.
sedimentFurthermore, since infiltration
transport. Furthermore, sinceisinfiltration
inversely related to slope,
is inversely it, therefore,
related to slope, means
it,
therefore,
that the flatmeans
slope that
areasthe flatcontribute
will slope areas will to
more contribute
groundwatermore recharge,
to groundwater
especiallyrecharge,
with especially
the extensive
with the extensive
distribution of sinksdistribution
in and around of sinks
theseinlowandareas.
around these low the
Conversely, areas. Conversely,
relief in Gaojiaping the relief in
watershed
Gaojiaping watershed is more evenly distributed with gentler slopes. No consecutive
is more evenly distributed with gentler slopes. No consecutive areas between two contour intervals areas between
two contour
were found tointervals
be higherwerethan found
10% toofbethehigher
entirethan 10% of the
watershed. entire
This watershed.
is very differentThis is very different
to Miaogou watershed
to Miaogou watershed where four classes (7–10) of contour intervals alone cover
where four classes (7–10) of contour intervals alone cover 15.02, 19.4, 17.6, and 12.45% of the watershed 15.02, 19.4, 17.6, and
12.45% of the watershed area, respectively (Figure 6). This highlights the
area, respectively (Figure 6). This highlights the impact of carbonate dissolution in creating a rugged impact of carbonate
dissolution in creating a rugged topography compared to the more homogenous crystalline
topography compared to the more homogenous crystalline watershed. The phenomenon of highly
watershed. The phenomenon of highly heterogeneous relief system in Miaogou makes for a more
heterogeneous relief system in Miaogou makes for a more dynamic inter-fluvial process within the
dynamic inter-fluvial process within the watershed thereby promoting more carbonate dissolution,
watershed thereby promoting more carbonate dissolution, higher drainage density (as is the case), and
higher drainage density (as is the case), and easy diversion of runoff processes into caves, sinks, and
easy diversion of runoff processes into caves, sinks, and dolines within the catchment.
dolines within the catchment.

Figure HypsographicCurves
6. Hypsographic
Figure 6. CurvesofofMiaogou
Miaogou(A)(A)
andand Gaojiaping
Gaojiaping (B) (B) Watershed.
Watershed.
3.2. Drainage Network Analysis
3.2. Drainage Network Analysis
Surface hydraulic transport mechanisms are driven by stream channels and their hierarchy which
Surface hydraulic transport mechanisms are driven by stream channels and their hierarchy
defines drainage density. Miaogou and Gaojiaping stream networks were extracted, and the main
which defines drainage density. Miaogou and Gaojiaping stream networks were extracted, and the
channel lengthlength
main channel as illustrated in Figure
as illustrated 7. Results
in Figure fromfrom
7. Results hydrographic
hydrographic computations
computations (Table 1) 1)
(Table show
the karst watershed area to be 45 km 2 and main stream channel to be about 15 km long while the
show the karst watershed area to be 45 km and main stream channel to be about 15 km long while
2

total length
the total of drainage
length of drainageis 30.86 km.
is 30.86 km.Gravelius
Gravelius index
index waswascalculated
calculatedtotobebe 1.93,
1.93, which
which is indicative
is indicative
ofofan
anelongated
elongated drainage systemwith
drainage system withthe
thepotential
potentialofofhaving
having low
low runoff
runoff time
time of concentration
of concentration and and
consequently gentler hydrograph peak. The drainage density is about 0.68 km 2 , which is same as the
consequently gentler hydrograph peak. The drainage density is about 0.68 km , which is same as the
2

stream
streamdensity
density is
is 0.68,
0.68, while the average
while the averagelength
lengthofofsurface
surfacerunoff
runoffis is 0.36.
0.36.
Water 2019, 11, 743 8 of 20
Water 2019, 11, x FOR PEER REVIEW 8 of 20

Figure Streamnetwork
7. Stream
Figure 7. networkorder
order (Strahler
(Strahler Order,
Order, [30])
[30]) of Miaogou
of Miaogou (A) Gaojiaping
(A) and and Gaojiaping (B) watershed.
(B) watershed.

On
On the other hand,
the other hand,Gaojiaping
Gaojiaping watershed
watershed area area is computed
is computed to beto54bekm km2 (Table
542 (Table 1) while
1) while its its
primary stream channel is 11.91 km long and shorter than that of Miaogou. The
primary stream channel is 11.91 km long and shorter than that of Miaogou. The total length of its total length of its
drainage is 25.06 km and also lesser than Miaogou, implying that overall, Gaojiaping
drainage 25.06 km and also lesser than Miaogou, implying that overall, Gaojiaping watershed has watershed has
shorter
shorter stream channelsfor
stream channels forthe
theconveyance
conveyance of of
its its total
total runoff
runoff as compared
as compared to Miaogou.
to Miaogou. FigureFigure
8 8
shows Gaojiaping
shows Gaojiapingtotohave a fan-shaped
have perimeter
a fan-shaped with awith
perimeter Gravelius index ofindex
a Gravelius 1.65, also lesseralso
of 1.65, thanlesser
that than
of Miaogou
that of Miaogouwatershed which
watershed implies
which a higher
implies runoff
a higher time time
runoff of concentration
of concentrationand consequently
and consequently
generate higher flow in shorter period. In addition, the drainage in Gaojiaping is characterized
generate higher flow in shorter period. In addition, the drainage in Gaojiaping is characterized by a by a
dendritic pattern.
dendritic pattern.

3.3. Geomorphic
3.3. Geomorphic Analysis
Analysis
Generally, geomorphic
Generally, geomorphicdynamics
dynamics is is
defined
definedby bythree factors
three as proposed
factors as proposed by Strahler [40], they
by Strahler [40], they
include gravitational stresses, molecular stresses and chemical processes. The impact of these forces
include gravitational stresses, molecular stresses and chemical processes. The impact of these forces is
is generally driven by fluvial process, which is governed principally by relief and underlying
generally driven by fluvial process, which is governed principally by relief and underlying geology.
geology. In Miaogou for example, which is mainly carbonate rocks, interaction with water with
In Miaogou for example, which is mainly carbonate rocks, interaction with water with presence of
presence of CO2 results in dissolution of the bed rock. This, however is not the case for metamorphic
CO 2 results
rocks in dissolution
in Gaojiaping, of the large
which require bed rock.
amount This, however
of heat (solar)istonot theitscase
drive for metamorphic
weathering process (rockrocks in
Gaojiaping, which require large amount of heat (solar) to drive its weathering
exfoliation). Thus, karst in humid-tropics which have abundant rainfall are characterized by karst process (rock exfoliation).
Thus,
towerskarst in humid-tropics
or karst cones (domes)which have abundant
that appear like conical rainfall
haystackarerounded
characterized by karst
hills [41] towers or karst
with elongated
cones
ground (domes) that appear
depressions extendinglikebetween
conical haystack
them which rounded
createshills [41] areas
low-flat with elongated ground
as in our study areadepressions
[42].
extending between them
Visual interpretation which8(A)
of Figure creates low-flat
reveals some areas
of theseas features,
in our studywherearea [42].19Visual
about interpretation
karst cones were of
delineated
Figure in the watershed.
8A reveals some of these The features,
highest cone
where peak is 1514
about 19 m above
karst mean
cones wereseadelineated
level whileinthe lowest
the watershed.
is 888m
The above
highest mean
cone peakseaislevel.
1514 Additionally,
m above mean five
seamajor
levelkarst
whiledepressions
the lowestwere alsoabove
is 888m delineated
meanandsea level.
measured to be
Additionally, fivearound
major 2.42,
karst0.48, 1.13 andwere
depressions 0.18 also
km delineated
2 respectively.and In measured
addition, 18 tosinkholes
be aroundwere2.42, 0.48,
identified
1.13 during
and 0.18 km2ground survey In
respectively. and where found
addition, to be mainly
18 sinkholes weredistributed
identifiedinduring
the upper half of
ground the and
survey
watershed. Also, it should be noted that Figure 8 has no vertical exaggeration,
where found to be mainly distributed in the upper half of the watershed. Also, it should be noted thus, the high
concavity of the karstified watershed further highlights the impact of geomorphic work on the overall
that Figure 8 has no vertical exaggeration, thus, the high concavity of the karstified watershed further
karstification phenomenon. The valleys along the depressions illustrate the weathering of folding
highlights the impact of geomorphic work on the overall karstification phenomenon. The valleys along
cavities. Gaojiaping on the hand has low concavity, the slopes are generally homogeneous perhaps
the depressions illustrate the weathering of folding cavities. Gaojiaping on the hand has low concavity,
due to higher resistance to geomorphic work by its crystalline formations (Figure 8B).
the slopes are generally homogeneous perhaps due to higher resistance to geomorphic work by its
crystalline formations (Figure 8B).
Water 2019, 11, 743 9 of 20
Water 2019, 11, x FOR PEER REVIEW 9 of 20

Figure8.8.3D
Figure 3DGeomorphometric
Geomorphometric illustration
illustration (A)
(A) Miaogou watershed characterized
Miaogou watershed characterizedby
byan
anelongated
elongated
shape with the presence of spring, sinkholes, caves, karst depressions and karst cones (B) Gaojiaping
shape with the presence of spring, sinkholes, caves, karst depressions and karst cones (B) Gaojiaping
watershed
watershedcharacterized
characterizedby
byaabroader
broaderfan-shaped
fan-shaped area
area with
with no
no karst
karst features.
features.
Table 1. Summary of parameters for both Miaogou (karst) and Gaojiaping (Non-karst) watersheds.
Table 1. Summary of parameters for both Miaogou (karst) and Gaojiaping (Non-karst) watersheds.
Domain Parameter Units Karst Non-Karst
Domain Parameter Units Karst Non-Karst
Dominant Geology - Sedimentary Metamorphic
Dominant Geology - Sedimentary Metamorphic
Sinkholes - 18 0
Sinkholes
Caves - - 18
4 00
Geomorphic Caves
Springs - - 24 00
Basin Shape
Springs - - High-Concavity
2 Low-Concavity
0
Geomorphic Cone Cones - 19 0
Low-
KarstBasin Shape
Depressions - - High-Concavity
4 0
Concavity
Area (A) km2 45 54
Cone Cones - 19 0
Perimeter (P) km2 46 43
Max.Karst Depressions
Altitude (ALTmax) m - 4
1770 0
1968
Area(ALTmin)
Min. Altitude (A) mkm² 45
430 54
762
Topographic
Avg. Altitude (ALTavg)
Perimeter (P) mkm² 1099
46 1365
43
Total Relief (RTotal)
Max. Altitude (ALTmax) m m 1340
1770 1206
1968
Watershed average slope % 31 40.35
Min. Altitude (ALTmin) m 430 762
Topographic Most Frequent Altitude m 898 1304.70
Avg. Altitude (ALTavg) m 1099 1365
Gravelius’s Shape Index (Cg) Un 1.93 1.65
Total Relief (RTotal) m 1340 1206
Drainage Network Avg. Slope % 22.56 19.34
Watershed
Main average(ML)
Channel Length slope km % 31
15 40.35
11.91
Most Frequent Altitude
Stream Order 1 Length km m 898
14.23 1304.70
16.80
Stream Order
Gravelius’s 2 Length
Shape Index (Cg) kmUn 5.63
1.93 6.35
1.65
Hydrographic Stream Order 3 Length
Drainage Network Avg. Slope km % 8.47
22.56 1.91
19.34
Stream Order 4 Length km 2.53 0
Main Channel Length (ML) km 15 11.91
Total Length of Drainage
Stream Order 1 Length kmkm 30.86
14.23 25.06
16.80
Network (Lt)
Stream Order
Drainage Density2 (Dd)
Length km/km km2 5.63
0.68 6.35
0.46
Stream
Stream Density
Order 3 Length km 0.68
8.47 0.28
1.91
Hydrographic Avg. Stream
Length of Surface Runoff kmkm 0.36 0.54
Order 4 Length 2.53 0
Total Length of Drainage
km 30.86 25.06
3.4. Hydrograph Recession Coefficient
Network(α) (Lt)
Estimation
Drainage Density (Dd) km/km² 0.68 0.46
From the streamflow data, fifteen (15) storms enumerated in Tables 2 and 3 were selected for
Stream Density 0.68 0.28
hydrograph recession analysis. The recession curve corresponds to the concave part of the hydrograph
Avg. Length of Surface Runoff km 0.36 0.54
following a peak i.e., when the crest of the hydrograph is at its highest and is not influenced by recharge
processes [43]. Using Equation (3), recession coefficients of successive limbs are calculated. In all, four
3.4. Hydrograph Recession Coefficient (α) Estimation
major limbs were segmented for the fifteen storms even though not all storms had four distinct stages.
From the
The recession streamflow
time data, hydrographs
of the selected fifteen (15) storms enumerated
ranged in Tables
from a minimum of 702 and
h (2.93 days)
were to
selected forh
over 289
hydrograph recession analysis. The recession curve corresponds to the concave part of the
Water 2019, 11, 743 10 of 20

(12 days). Here, it is noteworthy, Bonacci [44] stated that one of the challenges in recession analysis is
the time-step selection, where measurement methods in the past often restrict the treatment to a time
unit of one day, however, he stresses that this unit can be too long and thus suggested future research
to focus on hourly time-step. For this study, we also found this hourly time step as well as the overall
duration of about 12 days suitable for the purpose of our study, which is to determine the impact of
geomorphological features on the drainage function of karst watersheds against non-karst watersheds.

Table 2. Table of recession coefficients and recession stages at Miaogou Karst Watershed.

Peak Flow Stage I Stage II Stage III Stage IV

No
(m3 /s) α(1/h) T(h) Q1 α(1/h) T(h) Q2 α(1/h) T(h) Q3 α(1/h) T(h) Q4
1 2.71 0.0760 26 0.38 0.0133 91 0.11 0.0036 ≥216 0.05
2 9.71 0.0760 31 0.92 0.0113 95 0.31 0.0044 196 0.13 0.0006 >132 0.12
3 19.86 0.1008 36 0.53 0.0149 86 0.51 0.0053 >150 0.07
4 1.28 0.0972 19 0.20 0.0104 100 0.07 0.0043 >162 0.04
5 12.65 0.0748 36 0.86 0.0154 94 0.20 0.0050 >127 0.11
6 51.75 0.0535 53 3.04 0.0219 75 0.59 0.0041 195 0.26 0.0005 >67 0.26
7 18.08 0.0775 36 1.11 0.0131 97 0.31 0.0038 >122 0.20
8 3.35 0.0401 36 0.79 0.0085 88 0.37 0.0045 >136 0.20
9 9.31 0.0751 26 1.32 0.0162 90 0.31 0.0020 >70 0.27
10 14.93 0.0709 25 2.54 0.0094 109 0.91 0.0051 >150 0.42
11 12.43 0.0537 47 1.00 0.0074 115 0.43 0.0026 260 0.22 0.0008 >220 0.18
12 9.51 0.0839 35 0.50 0.0107 81 0.21 0.0021 279 0.12 0.0009 >352 0.09
13 7.25 0.0481 47 0.76 0.0059 120 0.37 0.0032 197 0.20 0.0006 >104 0.19
14 7.98 0.0539 28 1.76 0.0106 97 0.63 0.0047 172 0.28 0.0009 >128 0.25
15 2.64 0.0434 30 0.72 0.0126 60 0.34 0.0033 212 0.17 0.0004 ≥1858 0.08
Average 0.0683 34 0.79 0.0121 93 0.35 0.0039 216 0.18 0.0007 0.17
Std Dev
0.0187 9.3 0.79 0.004 15.1 0.22 0.0011 38.9 0.10 0.0002 0.07
(σ)

Table 2 shows the results for recession coefficients of Miaogou for four (4) major recession stages
(limbs) of the storms. The average value for the first segment is 0.0683 for an average period of 34 h
while the second stage has an average value of 0.0121 for an average of 93 h. For the third and fourth
stages, the recession coefficient drastically depletes to third and fourth order magnitudes at 0.0039 and
0.0007 respectively. The highest coefficients for both first and second stage segments are 0.1008, 0.0149,
and 0.0972, 0.0104 for storms 3 and 4, respectively. For Gaojiaping, the average value for the first
segment is 0.0405 for an average period of 31 h while the second stage has an average value of 0.0076
for an average of 96 h. For the third and fourth stages, the recession coefficient drastically depletes to
third and fourth order magnitudes at 0.0026 and 0.0004, respectively. The highest coefficients for both
first and second stage segments are 0.0630, 0.0100, and 0.0615, 0.0150 for storms 8 and 2, respectively.
Figure 9 shows the average recession coefficients of the four stages of the recession curves of the
two catchments.
The contrast between the two watersheds is highest during the first and second stages of the
recession limb which is characterized by the fast flow drained immediately after the peak flow. While
the slow flow represented by the third and fourth stages are delayed. This phenomena is principally in
conformity “matrix-restrained flow” and “conduct-influenced flow” (fast and slow flows respectively)
regimes elucidated by Kovács et al. [45]. Fiorillo [43] asserts that the “matrix-restrained flow” is
assumed to discharge instatenously to the outlet and therefore, requires well-developed conduit
pervading the entire karst system while the “conduct-influenced flow” is slower and thus, modelled
as an equal porous medium (non-karst). The longest period analyzed for Miaogou was greater than
1858 h which had maintained a recession coefficient of about 0.0004, while that of Gaojiaping is slightly
shorter at about 1729 h in Gaojiaping where the recession coefficient was steady about at 0.0002 due
the denser rock matrix.
 Gaojiaping is slightly shorter at about 1729 hours in Gaojiaping where the recession coefficient was
steady about at 0.0002 due the denser rock matrix.

Table 3. Table of recession coefficients and recession stages for Gaojiaping Watershed.
Water
No2019, 11, 743
Peak Stage Ⅰ Stage Ⅱ Stage Ⅲ Stage Ⅳ 11 of 20
flow
α(1/h) T(h) Q1 α(1/h) T(h) Q2 α(1/h) T(h) Q3 α(1/h) T(h) Q4
（m³/s）
1 11.21 3. Table
Table of recession
0.0444 27 coefficients
3.38 0.0095and recession
93 1.40stages for Gaojiaping
0.0065 ≥128 Watershed.
0.61
2 41.39 0.0615 32 5.78 0.0150 70 2.02 0.0050 ＞84 1.33
3No 4.89Flow 0.0195Stage27
Peak I 2.89 Stage II98
0.0094 1.15 Stage III
0.0024 ＞87 0.93 Stage IV
(m3 /s) α(1/h) T(h) Q1 α(1/h) T(h) Q2 α(1/h) T(h) Q3 α(1/h) T(h)＞ Q4
4 2.37 0.0213 22 1.48 0.0044 119 0.88 0.0012 194 0.70 0.0002 0.67
185
1 11.21 0.0444 27 3.38 0.0095 93 1.40 0.0065 ≥128 0.61
＞
52 41.39
7.27 0.0615 24
0.0537 32 5.78 0.0075
2.00 0.0150 70 74 2.02
1.15 0.0050
0.0013 >84220 1.33
0.86 0.0006 0.72
3 4.89 0.0195 27 2.89 0.0094 98 1.15 0.0024 >87 0.93 294
≥154
64 2.37
1.79 0.0213 33
0.0129 22 1.48 0.0030
1.17 0.0044 119
109 0.88
0.84 0.0012
0.0010 194271 0.70
0.64 0.0002
0.0002 >185 0.67
0.47
5 7.27 0.0537 24 2.00 0.0075 74 1.15 0.0013 220 0.86 0.0006 >2945 0.72
76 1.57
1.79 0.0072
0.0129 47 33 1.12
1.17 0.0032 104
0.0030 109 0.80
0.84 0.0014 271
0.0010 210 0.64
0.60 0.0002 ＞67 0.47
0.0005 ≥1545 0.58
87 9.10
1.57 0.0630
0.0072 22 47 2.28
1.12 0.0100 110
0.0032 104 0.76
0.80 0.0025 210
0.0014 ≥219 0.60
0.44 0.0005 >67 0.58
98 15.33
9.10 0.0496
0.0630 29 22 3.64
2.28 0.0108 108
0.0100 110 1.13
0.76 0.0040 ≥≥125
0.0025 0.69
219 0.44
109 15.33
16.83 0.0496 33
0.0373 29 3.64 0.0070
4.91 0.0108 108
104 1.13
2.37 0.0036 ≥125
0.0040 ≥97 0.69
1.67
1110 16.83
7.51 0.0373 37
0.0370 33 4.91 0.0037
1.91 0.0070 104 96 2.37
1.34 0.0011 ≥≥263
0.0036 97 1.67
1.00
11 7.51 0.0370 37 1.91 0.0037 96 1.34 0.0011 ≥263＞ 1.00
1212 11.53
11.53 0.0536
0.0536 28 28 2.57
2.57 0.0057
0.0057 98 98 1.47
1.47 0.0019 >120
0.0019 1.17
120 1.17
1313 37.18
37.18 0.0604 35
0.0604 35 4.49 0.0059
4.49 0.0059 89 89 2.66
2.66 0.0030 ≥≥194
0.0030 194 1.48
1.48
14 2.75 0.0270 35 1.07 0.0052 86 0.68 0.0007 285 0.56 0.0002 ≥1729≥172 0.40
1415 2.75
36.16 0.0270
0.0603 35 30 1.07
5.92 0.0052
0.0143 75 86 0.68
2.03 0.0007 208
0.0038 285 0.92
0.56 0.0005
0.0002 >96 0.40
9 0.88
Average
15 36.16 0.0405 30
0.0603 31 2.97 0.0143
5.92 0.0076 96 75 1.38
2.03 0.0026
0.0038 235208 0.91
0.92 0.0004
0.0005 ＞96 0.620.88
Average
Std Dev 0.0405 31 2.97 0.0076 96 1.38 0.0026 235 0.91 0.0004 0.62
Std (σ)
Dev 0.019 6.5 1.66 0.0038 14.5 0.62 0.0017 37.4 0.36 0.0002 0.17
0.019 6.5 1.66 0.0038 14.5 0.62 0.0017 37.4 0.36 0.0002 0.17
(σ)

Figure 9. Average recession coefficients of the four stages of the recession curves.
Figure 9. Average recession coefficients of the four stages of the recession curves.
3.5. ANOVA on Recession Coefficient
The analysis of variance (ANOVA) technique was used to determine the significant statistical
difference between the recession coefficients of both Miaogou and Gaojiaping watersheds. Results from
Tables 2 and 3 for each recession stage of both watersheds were individually analysed as a single group
with an alpha level of 0.1. It was observed that comparing each recession segment would statistically
be better as opposed to lumping all the recession stages together. This is done with a view to ensuring
singular groupings for the ANOVA can identify which stage or stages have the most statistical variance
and therefrom infer which watershed attributes have the major impact on streamflow.
The first stage of the recession curve, with a total of 15 records for each watershed showed
statistically significant difference in the recession coefficient of the two watersheds as determined by
one-way ANOVA at (F (1, 28) = 16.232, P = 0.000389). Hence, the null hypothesis is rejected based on
the alpha value set at 0.1, and the P value being significantly lower at 3.89 × 10−4 . Similarly, the second
stage of the recession curve was analysed where the one-way ANOVA was found to have statistically
significant difference at (F (1, 28) = 9.884, P = 0.00392). The P value of 3.92 × 10−3 was found to be
 The first stage of the recession curve, with a total of 15 records for each watershed showed
statistically significant difference in the recession coefficient of the two watersheds as determined by
one-way ANOVA at (F (1, 28) = 16.232, P = 0.000389). Hence, the null hypothesis is rejected based on
the alpha value set at 0.1, and the P value being significantly lower at 3.89 × 10 − 4. Similarly, the
Water 2019, 11, 743
second stage of the recession curve was analysed where the one-way ANOVA was found to12have of 20

statistically significant difference at (F (1, 28) = 9.884, P = 0.00392). The P value of 3.92 × 10 − 3 was
found to be lower
significantly significantly
than the lower than the
alpha value. alpha
Finally, forvalue. Finally,
the third for the
stage, there was third stage, theresignificant
the statistically was the
statistically significant difference between the recession coefficient of
difference between the recession coefficient of the two watersheds as determined by one-way ANOVA the two watersheds as
determined
at (F (1, 28) = by5.859,
one-way ANOVAItatis(Fnoteworthy
P = 0.0222). (1, 28) = 5.859,
thatPfrom
= 0.0222). It is noteworthy
the third that fromdifference
stage, the statistical the third
stage, the statistical difference begins to reduce in magnitude as the degree of
begins to reduce in magnitude as the degree of variance between the alpha and P value begins to close variance between the
alpha
up. andtherefore
It can P value begins to closefrom
be concluded up. visual
It can interpretation
therefore be concluded
of Figure 9from visual
that for the interpretation
third and fourth of
Figure 9 that for the third and fourth stages,
stages, the statistical variance is not extremely wide. the statistical variance is not extremely wide.
Figure 10
Figure 10 also
also shows
showsthe thevariability
variabilitydischarge
dischargeratio
ratioln(Q
ln(Qmax/Q/Qmean) )of the karst watershed to be
max mean of the karst watershed to be
much higher
much higher by by an
an order
order of
of nearly
nearly two
two magnitudes.
magnitudes. In In addition,
addition, the
the karst
karst discharge
discharge variability
variability shows
shows
an exponential growth from 2.26 to 4.14 in about 220 hours whereas the non-karst
an exponential growth from 2.26 to 4.14 in about 220 h whereas the non-karst discharge, at about discharge, at about
220 h
220 hours
begins begins
to decay to decay
rather rather
rapidly. Therapidly. Thedischarge
variability variability discharge
ratio has often ratio
beenhasusedoften been usedthe
to determine to
determine the
“flashiness” of a“flashiness”
spring, which of adistinguishes
spring, which distinguishes
diffuse (slow) and diffuse
conduit (slow)
(fast)and
flowsconduit
[46]. (fast) flows
[46].

Figure 10. Variability discharge ratio ln(Qmax /Qmean ) showing higher discharge rates with respect to
Figure 10. Variability discharge ratio ln(Qmax/Qmean) showing higher discharge rates with respect to
time in karst watershed.
time in karst watershed.
3.6. Streamflow Component Separation
3.6. Streamflow Component Separation
Three typical storms characterized by low, medium, and high volumes of stream discharges were
Three
selected for typical storms characterized
the streamflow by low, medium,
components separation. This is toand
have high volumes
a balance of stream
between discharges
dry period and
wereperiods
wet selectedduring
for thethe
streamflow
hydrologiccomponents
year. Usingseparation.
EquationsThis is toas
(5)–(7), have
wella balance betweencoefficients
as the recession dry period
and wet periods
generated in Tables during
2 and the hydrologic the
3, respectively, year. Using Equations
streamflow recession(5–7), as well
segments andastheir
the respective
recession
volumetric flows, computed. The initial and steepest part of the recession should characterize a
different hydraulic flow type, as it generally reflects the fall of the hydraulic head in the conduit
network. This is because the hydraulic behavior depends on two main hydraulic heads which
characterize karst aquifers, and connected to the peculiar structure of the medium [47]. The result,
shown in Tables 4 and 5 were quantified based on the function that higher recession coefficient implies
higher water drainage capacity per unit time which is primarily attributed to topographical and
geomorphological characteristics. The geology is particularly critical here since it defines the drainage
porosity which in effect is more impactful after the first overland discharge. Also, the function of
the aquifer specific yield is known to be given by the ratio of the volume of water drained to the
rock volume. Hence, for Miaogou, we defined four (4) classes of flow regime that contribute to
streamflow; (1) overland flow, (2) karst drainage flow, (3) medium fracture flow and (4) micro fracture
flow. For Gaojiaping the classes are (1) overland flow, (2) macro fracture flow, (3) medium fracture
flow, and (4) micro fracture flow. The streamflow separation results show that overland flow and karst
drainage flow have the largest percentages for all three storms that were evaluated. This phenomenon
Water 2019, 11, 743 13 of 20

correlates well with the drainage characteristics of karst aquifers having multiple porosity structures
due to karstification [48]. Conversely, in Gaojiaping, micro fracture flow, which is generally through
regolith deposits and soil sediments contribute the most significant percentage of water ratio. Thus
we can conclude that the aquifer matrix of the non-karst basin is more homogenous than that of the
karst basin. When the volumes of each stream component are examined and contrasted between the
case areas, the disparity becomes apparent. For example, Figure 11 shows the disparity in combined
overland and karst drainage flows when compared to the overland, and large fracture flows in between
the case study areas. In the first stage of the recession, there is a difference of about 40 m3 while the
second stage has a difference of about 30 m3 . This implies that the karst watershed tends to discharge
larger volumes of water per unit time over the surface as well as subsurface at nearly the same rate.
Bonacci [44] elucidated this phenomenon with cases from Switzerland and France [49] and Russia [50]
as well as other parts of the world [12], which show that the water retention capacity in karst areas
may not be great. This corroborates the findings in this study, even though, the high surficial flow
observed in Miaogou, is not primarily attributed to any particular karst feature but more to the basin
geometry and morphology, which are ultimately defined by the underlying carbonate system. Hence,
there is an indirect correlation but it is merely circumstantial since other geological formations can
produce basins geometries that discharge runoff quickly.

Table 4. Table of streamflow component separation for Miaogou watershed.

Miaogou α, duration 0.0535 (0,54] 0.0219 (55,129] 0.0041 (130,324] 0.0005 (325,∞]
Stages Total I II III IV
Flow Type V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 %
Overland Flow 262.99 39.96 262.99 68.24
Karst Drainage Flow 124.06 18.85 101.94 26.45 22.12 50.30
Medium Fracture Flow 34.54 5.25 13.15 3.41 12.02 27.35 9.37 28.16
Micro Fracture Flow 236.53 35.94 7.31 1.90 9.83 22.35 23.89 71.84 195.50 100.00
Total 658.13 100.00 385.40 100.00 43.97 100.00 33.26 100.00 195.50 100.00
Miaogou α, duration 0.0839 (0,36] 0.0107 (37,117] 0.0021 (118,396] 0.0009 (397,∞]
Stages Total I II III IV
Flow Type V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 %
Overland Flow 30.79 33.06 30.79 77.52
Karst Drainage Flow 9.18 9.86 5.31 13.37 3.87 34.95
Medium Fracture 8.54 9.17 1.55 3.91 2.79 25.22 4.20 24.47
Micro Fracture Flow 44.62 47.91 2.07 5.20 4.41 39.83 12.96 75.53 25.18 100.00
Total 93.13 100.00 39.72 100.00 11.08 100.00 17.15 100.00 25.18 100.00
Miaogou α, duration 0.0434 (0,31] 0.0126 (31,91] 0.0033 (91,303] 0.0004 (303,∞]
Stages Total I II III IV
Flow Type V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 %
Overland Flow 5.08 3.79 5.08 37.31
Karst Drainage Flow 7.23 5.38 4.49 32.98 2.74 28.85
Medium Fracture Flow 12.33 9.18 2.53 18.54 3.86 40.66 5.94 38.03
Micro Fissure Flow 109.69 81.66 1.52 11.17 2.89 30.48 9.68 61.97 95.60 100.00
Total 134.34 100.00 13.63 100.00 9.49 100.00 15.62 100.00 95.60 100.00

Table 5. Table of streamflow Component Separation for Gaojiaping watershed.

Gaojiaping α, duration 0.0537 (0,25] 0.0075 (26,99] 0.0013 (100,319] 0.0006 (320,∞]
Stages Total I II III IV
Flow Type V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 %
Overland Flow 16.14 2.74 16.14 44.92
Macro Fracture Flow 18.99 3.23 8.74 24.32 10.25 25.05
Medium Fracture Flow 16.28 2.77 2.37 6.59 5.73 13.99 8.18 10.75
Micro Fracture Flow 536.51 91.26 8.68 24.17 24.95 60.96 67.97 89.25 434.90 100.00
Total 587.91 100.00 35.92 100.00 40.93 100.00 76.15 100.00 434.90 100.00
Water 2019, 11, 743 14 of 20

Table 5. Cont.

Gaojiaping α, duration 0.0129 (0,34] 0.003 (35,142] 0.001 (143,413] 0.0002 (414,∞]
Stages Total I II III IV
Flow Type V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 %
Overland Flow 2.89 0.53 2.89 20.79
Macro Fracture Flow 1.57 0.29 1.13 8.13 0.44 1.56
Medium Fracture Flow 29.34 5.35 3.46 24.88 9.02 31.88 16.85 23.66
Micro Fracture Flow 514.53 93.84 6.43 46.20 18.83 66.55 54.37 76.34 434.90 100.00
Total 548.33 100.00 13.91 100.00 28.30 100.00 71.22 100.00 434.90 100.00
Gaojiaping α, duration 0.0270 (0,25] 0.0052 (26,110] 0.0007 (111,395] 0.0002 (396,∞]
Stages Total I II III IV
Flow Type V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 % V,104 m3 %
Overland Flow 4.99 0.42 4.99 27.84
Macro Fracture Flow 12.05 1.02 4.82 26.88 7.23 20.52
Medium Fracture Flow 17.29 1.46 1.73 9.64 5.35 15.18 10.21 12.36
Micro Fracture Flow 1149.10 97.10 6.40 35.65 22.67 64.31 72.39 87.64 1047.65 100.00
Water 2019,Total
11, x FOR PEER1183.43
REVIEW 100.00 17.94 100.00 35.25 100.00 82.59 100.00 1047.65 15 of 20
100.00

Figure
Figure 11. Contrast
11. Contrast in Overland
in Overland + Karst
+ Karst Drainage
Drainage Flow Flow Volume
Volume vs. Overland
vs. Overland + Macro
+ Macro Fracture
Fracture Flow
Flow Volume.
Volume.

Similarly, Figure 12 shows the disparity in micro fracture flow for both watersheds where the
Similarly, Figure 12 shows the disparity in micro fracture flow for both watersheds where the
non-karst catchment has a larger volumetric discharge rate. This property of non-karst watershed is
non-karst catchment has a larger volumetric discharge rate. This property of non-karst watershed is
even more pronounced during the wetter phases beginning from the first stage through to the third
even more pronounced during the wetter phases beginning from the first stage through to the third
stage which is drier due to the absence of precipitation. From Figure 11 the drainage volumes for all
stage which is drier due to the absence of precipitation. From Figure 11 the drainage volumes for all
for flow types for each of our case study areas and stacked against each other. It can be concluded that
for flow types for each of our case study areas and stacked against each other. It can be concluded
the karst watershed has a faster and more robust drainage system due to the impacts of topography
that the karst watershed has a faster and more robust drainage system due to the impacts of
and especially geomorphology, as exemplified by the karst drainage attributes (KDA) such as caves,
topography and especially geomorphology, as exemplified by the karst drainage attributes (KDA)
sinks, and closed depressions. The non-karst watershed, on the other hand, has the potential to hold
such as caves, sinks, and closed depressions. The non-karst watershed, on the other hand, has the
more groundwater reserves despite being a crystalline hard rock formation.
potential to hold more groundwater reserves despite being a crystalline hard rock formation.
stage which is drier due to the absence of precipitation. From Figure 11 the drainage volumes for all
for flow types for each of our case study areas and stacked against each other. It can be concluded
that the karst watershed has a faster and more robust drainage system due to the impacts of
topography and especially geomorphology, as exemplified by the karst drainage attributes (KDA)
Water
such 2019, 11, 743 sinks, and closed depressions. The non-karst watershed, on the other hand, has
as caves, 15 of 20
the
potential to hold more groundwater reserves despite being a crystalline hard rock formation.

Figure 12.
Figure Contrast in
12. Contrast in Micro
Micro Fracture
Fracture Flow
Flow Volume
Volume for
for all
all four
four recession
recession stages.
stages.
3.7. Implications on Groundwater Availability and Quality in Karst Areas
3.7. Implications on Groundwater Availability and Quality in Karst Areas
High recession coefficient watershed for successive stages of the recession hydrograph has
High recession coefficient watershed for successive stages of the recession hydrograph has
implications on the availability of groundwater especially during dry periods with no precipitation [17].
implications on the availability of groundwater especially during dry periods with no precipitation
Fiorillo et al. [17] also reports that when the recession coefficient decreases or increases, the aquifer
[17]. Fiorillo et al. [17] also reports that when the recession coefficient decreases or increases, the
is emptying more slowly or more quickly than expected, respectively. This is due to the significant
influence in hydraulic potential of the aquifers storage capacity resulting from an increase in effective
porosity in depth. In the case of Miaogou, higher recession coefficient, along with growing population
and industrial activities in the area puts the groundwater system under more pressure. Recent reports
from Huangliang Township state that, during the dry winter season, groundwater reserves have been
found to be short.
Table 6 shows the average percentile of each streamflow component from the three storms that
were selected in this study. In Miaogou, overland flow accounts for 68.12%, karst drainage flow (KDF)
accounts for 25.47%, medium fracture flow is 3.93% while micro fracture flow is only 2.48%. However,
in Gaojiaping, only 35.44% of its average streamflow is from overland flow compared to Miaogou’s
68.12%. Furthermore, the delayed flow, that is, micro fracture flow is 31.74% making it more than
twelve times the size of Miaogou karst watershed. This highlights low groundwater yield in Miaogou
due to drainage excess by KDAs. These two characteristics agree with morphometric analyses in
above where it was established that Miaogou has a higher concavity that dominates about 50% of its
cumulated surface while its average stream network slope is higher than that of Gaojiaping watershed.
Thus, we see a higher overland flow percentage as well as a higher macro fracture percentage while
the converse is seen in the medium and micro fracture flows, respectively. This disparity between the
karst and non-karst watershed is shown in Figure 13, illustrating the inverse relationship between the
respective streamflow components and also strongly highlights the groundwater resources availability
and potential of both catchments.
Water 2019, 11, 743 16 of 20

Table 6. Flow type statistics for all three storms and their average percentages.

Miaogou
Average Flow Percentage of Average
Flow Type Storm 1 Storm 2 Storm 3
Type Flow Type (%)
Overland Flow 262.99 30.79 5.08 99.62 68.12
Karst Drainage Flow 101.94 5.31 4.49 37.25 25.47
Medium Fracture Flow 13.15 1.55 2.53 5.74 3.93
Micro Fracture Flow 7.31 2.07 1.52 3.63 2.48
Total 385.39 39.72 13.62 146.24 100.00
Gaojiaping
Overland Flow 16.14 2.89 4.99 8.01 35.44
Macro Fracture Flow 8.74 1.13 4.82 4.90 21.67
Medium Fracture Flow 2.37 3.46 1.73 2.52 11.15
Micro Fracture Flow 8.68 6.43 6.4 7.17 31.74
Water 2019, 11, Total 35.93
x FOR PEER REVIEW 13.91 17.94 22.59 100.00 17 of 20

Figure13.
Figure Theplot
13.The plotof
ofdisparity
disparityin
inflow
flowtype
type average
average percentages
percentages contribution
contribution to
to streamflow.
streamflow.

Besides the availability of groundwater reserves, another significant concern is groundwater

Besides the availability of groundwater reserves, another significant concern is groundwater
quality. Malik’s [18] revised aquifer pollution sensitivity index show that multiple recession stages
quality. Malik’s [18] revised aquifer pollution sensitivity index show that multiple recession stages
with high recession coefficients values imply high heterogeneity and thus higher sensitivity to pollution.
with high recession coefficients values imply high heterogeneity and thus higher sensitivity to
Also, Jakada et al., [5] carried out an intrinsic vulnerability assessment of Miaogou catchment area
pollution. Also, Jakada et al., [5] carried out an intrinsic vulnerability assessment of Miaogou
based on the EPIK parametric method and found that most of the catchment is highly vulnerable
catchment area based on the EPIK parametric method and found that most of the catchment is highly
to pollution. Specifically, Huangliang Township, which has recently been reporting groundwater
vulnerable to pollution. Specifically, Huangliang Township, which has recently been reporting
shortages, falls under the “very high” vulnerability index. Thus, results from this work validate this
groundwater shortages, falls under the “very high” vulnerability index. Thus, results from this work
conclusion as well as establish the fact that groundwater quality will continue to deteriorate until
validate this conclusion as well as establish the fact that groundwater quality will continue to
drastic environmental management practices are enforced.
deteriorate until drastic environmental management practices are enforced.
4. Conclusions
4. Conclusion
This is the first attempt at a comparative hydro-geomorphological and hydrograph recession
This is the first attempt at a comparative hydro-geomorphological and hydrograph recession
analysis between a karst and non-karst watershed. By comparing these two distinct watersheds, under
analysis between a karst and non-karst watershed. By comparing these two distinct watersheds,
common precipitation conditions, a major contrast was observed based on their hydrograph recession
under common precipitation conditions, a major contrast was observed based on their hydrograph
characteristics, which has direct correlation with shallow and deep aquifer drainage properties. Results
recession characteristics, which has direct correlation with shallow and deep aquifer drainage
show that the karst watershed in contrast to non-karst watershed tended to drain larger volume
properties. Results show that the karst watershed in contrast to non-karst watershed tended to drain
larger volume of precipitation over shorter period. This is due to larger amounts of interflow
influenced by large fractures, while other karst features such as sinks and closed depressions act as
siphoning conduits that transmit runoff rapidly through well-developed underground cavities to
springs which drain into the surface water system (streams and river channels). These common
features usually found in highly karstified areas, referred to in this study as karst drainage attributes
Water 2019, 11, 743 17 of 20

of precipitation over shorter period. This is due to larger amounts of interflow influenced by large
fractures, while other karst features such as sinks and closed depressions act as siphoning conduits
that transmit runoff rapidly through well-developed underground cavities to springs which drain
into the surface water system (streams and river channels). These common features usually found
in highly karstified areas, referred to in this study as karst drainage attributes (KDA), generate what
we also termed karst drainage flow (KDF). The KDF is defined by the second limb of the recession
hydrograph and also shares similar properties as overland flow due to high discharge volume as well
high recession coefficient. Kresic and Bonacci [39] also argue that the steep curve one the recession
limb represents turbulent drainage of large fractures and conduits, followed by a transitional portion
that is less turbulent and reflects the contribution of smaller fractures and rock matrix, ending with the
slowly decreasing curve, the so-called master recession curve, where the drainage of rock matrix and
small fissures is dominant. However, it must be stated that this study by no means attempts to draw
conclusions that hydrograph recession analysis fully describes the internal groundwater structure, as
the initial catchment area is karst can never be reliably defined [44]. Consequently, the statistically
significant variance in recession coefficients between the second stages of the two watersheds simply
emphasizes a strong relationship between the recession coefficient and KDAs, as well as catchment
morphometric properties, as posited by references [11–13].
Additionally, with increasing pressure on groundwater reserves globally, our results seek to
draw attention to karst aquifer water retention abilities [12,49,50] relative to groundwater quality,
not only locally in the study area but also in other parts of the world. Since higher recession
coefficient index for successive recession segments (regimes) imply strong sensitivity to groundwater
pollution, the environmental implications cannot be ignored. Escolero et al. [51] state that “the high
permeability of karst terrains result in many practical problems such as scarcity and poor predictability
of groundwater supplies; scarcity of surface water bodies; instability of the ground (due to the presence
of sinkholes); leakage of surface water reservoirs, and an unreliable waste disposal environment”.
Previous groundwater vulnerability studies established that the Miaogou area has a high intrinsic
susceptibility to pollution due to karstification, hence, sustainable water management measures must
be adopted in the area.

Author Contributions: Conceptualization, H.J. and Z.C.; Methodology, H.J. and M.L.; Software, Z.W.; Validation,
H.J., Z.C.; Formal Analysis, H.J., Z.C. and M.L.; Investigation, H.J.; Resources, H.Z.; Data Curation, M.L. and H.Z.;
Writing-Original Draft Preparation, H.J.; Writing-Review & Editing, H.J., Z.C. and M.H.; Visualization, M.H. and
Z.W.; Supervision, Z.C.; Project Administration, H.Z.; Funding Acquisition, H.Z.
Funding: This research was supported through the project of the National Key Research and Development
Program of China [2017YFC0406105], Natural Science Foundation China [41807199], Karst Dynamics Laboratory,
MLR and GZAR (Institute of Karst Geology, CAGS) [KDL201702], Natural Science Foundation of Hubei
[2018CFB170], China Geological Survey [12120113103800, DD20160304] and the Fundamental Research Funds for
the Central Universities, China University of Geosciences, Wuhan [CUG170670].
Conflicts of Interest: The Authors declare no conflict of interest. The funders had no role in the design of the
study; in the collection, analyses or interpretation of data; in the writing of the manuscript and in the decision
to publish.

References
1. Hartmann, A.; Goldscheider, N.; Wagener, T.; Lange, J.; Weiler, M. Karst water resources in a changing world:
Review of hydrological modeling approaches. Rev. Geophys. 2014, 52, 218–242. [CrossRef]
2. Ford, D.; Williams, P. Karst Hydrogeology and Geomorphology; John Wiley and Sons: Hoboken, NJ, USA, 2007;
ISBN 978-0-470-84996-5.
3. Stephen, R.; Jeannin, P.Y.; Alexander, E.C.; Davies, G.J.; Schindel, G.M. Contrasting definitions for the term
‘karst aquifer’. Hydrogeol. J. 2017, 25, 1237–1240. [CrossRef]
4. Kresic, N.; Mikszweski, A. Hydrogeological Conceptual Site Models: Data Analysis and Visualization; CRC Press:
Boca Raton, FL, USA, 2013.
Water 2019, 11, 743 18 of 20

5. Jakada, H.; Chen, Z.; Luo, Z.; Zhou, H.; Luo, M.; Ibrahim, A.; Tanko, N. Coupling Intrinsic Vulnerability
Mapping and Tracer Test for Source Vulnerability and Risk Assessment in a Karst Catchment Based on
EPIK Method: A Case Study for the Xingshan County, Southern China. Arab. J. Sci. Eng. 2018, 44, 377–389.
[CrossRef]
6. Singhal, B.B.S. Nature of Hard Rock Aquifers: Hydrogeological Uncertainties and Ambiguities.
In Groundwater Dynamics in Hardrock Aquifers; Ahmed, S., Jayakumar, R., Salih, A., Eds.; Springer: Dordrecht,
Germany, 2008. [CrossRef]
7. Price, K. Effects of watershed topography, soils, land use, and climate on baseflow hydrology in humid
regions: A review. Prog. Phys. Geogr. Earth Environ. 2011, 35, 465–492. [CrossRef]
8. Christina, T.; Gordon, G.E. A geological framework for interpreting the low-flow regimes of Cascade streams.
Willamette River Basin Oregon Water Resour. Res. 2004, 40, W04303.
9. Clark, M.P.; Rupp, D.E.; Woods, R.A.; Meerveld, H.J.T.-V.; Peters, N.E.; Freer, J.E. Consistency between
hydrological models and field observations: Linking processes at the hillslope scale to hydrological responses
at the watershed scale. Hydrol. Process. 2009, 23, 311–319. [CrossRef]
10. Sayama, T.; McDonnell, J.J.; Dhakal, A.; Sullivan, K. How much water can a watershed store? Hydrol. Process.
2011, 25, 3899–3908. [CrossRef]
11. Avdagić, I. Oticanje u KrašKim HidrološKim Sistemima (Modelling of Runoff from Karstic Catchments).
Ph.D. Thesis, Sarajevo University, Sarajevo, Yugoslavia, 1987.
12. Bonacci, O. Karst Hydrology; Springer: Herdelberg, Germany, 1987.
13. Bonacci, O.; Jelin, J. Identification of a karst hyrological system in the Dinaric karst (Yugoslavia). Hydrol. Sci. J.
1988, 33, 483–497. [CrossRef]
14. Soulios, G. Contribution à l’étude des courbes de récession des sources karstiques: Exemples du pays
Hellénique. J. Hydrol. 1991, 127, 29–42. [CrossRef]
15. Kullman, E. Krasovo–Puklinové vody. Karst-Fissure Waters; Geologický ústav Dionýza Štúra: Bratislava,
Czechoslovakia, 1990.
16. Luo, M.; Chen, Z.; Yin, D.; Jakada, H.; Huang, H.; Zhou, H.; Wang, T. Surface flood and underground flood in
Xiangxi River Karst Basin: Characteristics, models, and comparisons. J. Earth Sci. 2016, 27, 15–21. [CrossRef]
17. Fiorillo, F.; Revellino, P.; Ventafridda, G. Karst aquifer draining during dry periods. J. Cave Stud. 2012, 74,
148–156. [CrossRef]
18. Malik, P. Assessment of regional karstification degree and groundwater sensitivity to pollution using
hydrograph analysis in the Velka Fatra Mountains, Slovakia. Environ. Earth Sci. 2006, 51, 707–711. [CrossRef]
19. Pankaj, A.; Kumar, P. GIS-based morphometric analysis of five major sub-watersheds of Song River,
Dehradun District, Uttarakhand with special reference to landslide incidences. J. Soc. Remote. Sens. 2009, 37,
157–166. [CrossRef]
20. Doctor, D.H.; Young, J.A. An Evaluation of Automated GIS Tools for Delineating Karst Sinkholes and
Closed Depressions from 1-Meter LIDAR-derived Digital Elevation Data. In Proceedings of the 13th
Multidisciplinary Conference on Sinkholes and the Engineering and Environmental Impacts of Karst,
Carlsbad, NM, USA, 6–10 May 2013.
21. Prafull, S.; Ankit, G.; Madhulika, S. Hydrological inferences from the watershed analysis of water resource
management using remote sensing and GIS techniques. Egypt. J. Remote Sens. Space Sci. 2014, 17, 111–121.
22. Olivera, F.; Maidment, D. GIS for hydrologic data development for the design of highway drainage facilities.
Transp. Res. Rec. 1998, 1625, 131–138. [CrossRef]
23. Olivera, F.; Maidment, D. GIS tools for HMS modeling support. In Hydrologic and Hydraulic Modeling Support
with Geographic Information Systems; Maidment, D., Djokic, D., Eds.; ESRI Press: Redlands, CA, USA, 2000;
pp. 85–112.
24. HEC. Geospatial Hydrologic Modeling Extension HEC-GeoHMS: Version 1.0 User’s Manual; Hydrologic
Engineering Center (HEC) Technical Report No CPD-77; U.S. Army Corps of Engineers: Champaign,
IL, USA, 2000.
Water 2019, 11, 743 19 of 20

25. Oliveira, E.D.; Crestani, A. Physiographic characterization of the watershed of stream Jandaia, Jandaia do
Sul/PR. Acta Geogra´Fica 2011, 5, 169–183.
26. EPA. Better Assessment Science Integrating Point and Nonpoint Sources (BASINS): Version 3.0 User’s Manual;
Office of Water Technical Report No EPA-823-B-01-001; Environmental Protection Agency: Washington, DC,
USA, 2001.
27. Di Luzio, M.; Srinivasan, R.; Arnold, J.G.; Neitsch, S.L. Soil and water assessment tool. In ArcView GIS
Interface Manual: Version 2002; International Center for Tropical Agriculture: Cali, Colombia, 2002.
28. Luzio, M.; Srinivasan, R.; Arnold, J.G. Integration of Watershed Tools and Swat Model into Basins. JAWRA J.
Am. Resour. Assoc. 2002, 38, 1127–1141. [CrossRef]
29. Dixon, B.; Uddameri, V. GIS and Geocomputation for Water Resource Science and Engineering; John Wiley & Sons,
Ltd.: Hoboken, NJ, USA, 2016.
30. Strahler, A.N. Quantitative Geomorphology of Drainage Basins and Channel Networks. In Hand Book of
Applied Hydrology; Te Chow, V., Ed.; McGraw Hill Book Company: New York, NY, USA, 1964.
31. O’Callaghan, J.F.; Mark, D.M. The extraction of drainage networks from digital elevation data. Comput. Vis.
Graph. Process. 1984, 28, 323–344. [CrossRef]
32. Jenson, S.K.; Domingue, J.O. Extracting topographic structure from digital elevation data for geographic
information system analysis. Photogramm. Eng. Remote Sens. 1988, 54, 1593–1600.
33. Jenson, S.K. Applications of hydrologic information automatically extracted from digital elevation models.
Hydrol. Process. 1991, 5, 31–44. [CrossRef]
34. Boussinesq, J. Recherches Théoriques sur l’écoulement des Nappes d’eau Infiltrées dans le sol et sur dÉbit
de Sources. J. De Mathématiques Pures Et Appliquées 1904, 10, 5–78.
35. Maillet, E. Essais d’Hydraulique Souterraine et Fluviale; Hermann: Paris, France, 1905. (In French)
36. Kullman, E. Nové metodické prístupy k riešeniu ochrany a ochranných pásiem zdrojov podzemných vôd v
horninových prostrediach s krasovo—puklinovou priepustnost’ou [New methods in groundwater protection
and delineation of protection zones in fissure-karst rock environment; in Slovak]. Podzemn. Voda 2000, 6,
31–41.
37. Malík, P.; Vojtkova, S. Use of recession-curve analysis for estimation of karstification degree and its
application in assessing overflow/underflow conditions in closely spaced karstic springs. Environ. Earth Sci.
2012, 65, 2245–2257. [CrossRef]
38. Miao, Z.L.; Miao, Z.Z. Application of Recession Equation in Groundwater Studies. Investig. Sci. Technol.
1984, 5, 1–6. (In Chinese)
39. Kresic, N.; Bonacci, O. Spring Discharge Hydrograph. In Groundwater Hydrology of Springs; Kresic, N.,
Stevanovic, Z., Eds.; Butterworth-Heinemann: Oxford, UK, 2010.
40. Strahler, A.N. Dynamic basis of geomorphology. Geol. Soc. Am. Bull. 1952, 63, 923–938. [CrossRef]
41. Lambert, A.R. Satellite Geology and Photogeomorphology: An Instructional Manual for Data Integration; Springer:
Berlin/Heidelberg, Germany, 2011; pp. 78–79, ISBN 978-3-642-20607-8.
42. Stevanovic, Z. Karst Aquifers—Characterization and Engineering. Prof. Pract. Earth Sci. 2015. [CrossRef]
43. Fiorillo, F. Tank-reservoir drainage as a simulation of the recession limb of karst spring hydrographs. Hyogeol. J.
2011, 19, 1009–1019. [CrossRef]
44. Bonacci, O. Karst springs hydrographs as indicators of karst aquifers. Hydrol. Sci. J. 1993, 38, 51–62.
[CrossRef]
45. Kovács, A.; Perrochet, P.; Király, L.; Jeannin, P. A quntitative method for characterization of karst aquifers
based on the spring hydrograph analysis. J. Hydrol. 2005, 303, 152–164. [CrossRef]
46. Florea, L.J.; Vacher, H. Springflow Hydrographs: Eogenetic vs. Telogenetic Karst. Ground Water 2006, 44,
352–361. [CrossRef]
47. Fiorillo, F. The Recession of Spring Hydrographs, Focused on Karst Aquifers. Resour. Manag. 2014, 28,
1781–1805. [CrossRef]
48. Birk, S.; Hergarten, S. Early recession behaviour of spring hydrographs. J. Hydrol. 2010, 387, 24–32. [CrossRef]
49. Burger, A.; Pasquier, F. Prospection et captage d’eau par fprages dans la vallée de la Brevine (Jura Suisse).
In Hydrogeology of Karstic Terrains; Burger, A., Dubertret, L., Eds.; UNESCO: Paris, France, 1984; Volume 1,
pp. 145–149.
Water 2019, 11, 743 20 of 20

50. Babushkina, V.D.; Lebedyanskaya, Z.P.; Plotnikov, I.I. The distinctive features of predicting total water
discharge from deep-level mines in fissure and karst rocks. In Hydrogeology of Karstic Terrains; Burger, A.,
Dubertret, L., Eds.; UNESCO: Paris, France, 1984; Volume 1, pp. 229–232.
51. Escolero, O.A.; Marin, L.E.; Steinich, B.; Pacheco, A.J.; Cabrera, S.A.; Alcocer, J. Development of a Protection
Strategy of Karst Limestone Aquifers: The Merida Yucatan, Mexico Case Study. Resour. Manag. 2002, 16,
351–367.

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC BY) license (http://creativecommons.org/licenses/by/4.0/).